Actual source code: fegeom.c
1: #include <petsc/private/petscfeimpl.h>
3: /*@C
4: PetscFEGeomCreate - Create a `PetscFEGeom` object to manage geometry for a group of cells
6: Input Parameters:
7: + quad - A `PetscQuadrature` determining the tabulation
8: . numCells - The number of cells in the group
9: . dimEmbed - The coordinate dimension
10: - faceData - Flag to construct geometry data for the faces
12: Output Parameter:
13: . geom - The `PetscFEGeom` object
15: Level: beginner
17: .seealso: `PetscFEGeom`, `PetscQuadrature`, `PetscFEGeom`, `PetscFEGeomDestroy()`, `PetscFEGeomComplete()`
18: @*/
19: PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscBool faceData, PetscFEGeom **geom)
20: {
21: PetscFEGeom *g;
22: PetscInt dim, Nq, N;
23: const PetscReal *p;
25: PetscQuadratureGetData(quad, &dim, NULL, &Nq, &p, NULL);
26: PetscNew(&g);
27: g->xi = p;
28: g->numCells = numCells;
29: g->numPoints = Nq;
30: g->dim = dim;
31: g->dimEmbed = dimEmbed;
32: g->isCohesive = PETSC_FALSE;
33: N = numCells * Nq;
34: PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);
35: if (faceData) {
36: PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n);
37: PetscCalloc6(N * dimEmbed * dimEmbed, &(g->suppJ[0]), N * dimEmbed * dimEmbed, &(g->suppJ[1]), N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]), N, &(g->suppDetJ[0]), N, &(g->suppDetJ[1]));
38: }
39: PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);
40: *geom = g;
41: return 0;
42: }
44: /*@C
45: PetscFEGeomDestroy - Destroy a `PetscFEGeom` object
47: Input Parameter:
48: . geom - `PetscFEGeom` object
50: Level: beginner
52: .seealso: `PetscFEGeom`, `PetscFEGeomCreate()`
53: @*/
54: PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom)
55: {
56: if (!*geom) return 0;
57: PetscFree3((*geom)->v, (*geom)->J, (*geom)->detJ);
58: PetscFree((*geom)->invJ);
59: PetscFree2((*geom)->face, (*geom)->n);
60: PetscFree6((*geom)->suppJ[0], (*geom)->suppJ[1], (*geom)->suppInvJ[0], (*geom)->suppInvJ[1], (*geom)->suppDetJ[0], (*geom)->suppDetJ[1]);
61: PetscFree(*geom);
62: return 0;
63: }
65: /*@C
66: PetscFEGeomGetChunk - Get a chunk of cells in the group as a `PetscFEGeom`
68: Input Parameters:
69: + geom - `PetscFEGeom` object
70: . cStart - The first cell in the chunk
71: - cEnd - The first cell not in the chunk
73: Output Parameter:
74: . chunkGeom - The chunk of cells
76: Level: intermediate
78: Note:
79: Use `PetscFEGeomRestoreChunk()` to return the result
81: .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
82: @*/
83: PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
84: {
85: PetscInt Nq;
86: PetscInt dE;
90: if (!(*chunkGeom)) PetscNew(chunkGeom);
91: Nq = geom->numPoints;
92: dE = geom->dimEmbed;
93: (*chunkGeom)->dim = geom->dim;
94: (*chunkGeom)->dimEmbed = geom->dimEmbed;
95: (*chunkGeom)->numPoints = geom->numPoints;
96: (*chunkGeom)->numCells = cEnd - cStart;
97: (*chunkGeom)->xi = geom->xi;
98: (*chunkGeom)->v = &geom->v[Nq * dE * cStart];
99: (*chunkGeom)->J = &geom->J[Nq * dE * dE * cStart];
100: (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq * dE * dE * cStart] : NULL;
101: (*chunkGeom)->detJ = &geom->detJ[Nq * cStart];
102: (*chunkGeom)->n = geom->n ? &geom->n[Nq * dE * cStart] : NULL;
103: (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL;
104: (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq * dE * dE * cStart] : NULL;
105: (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq * dE * dE * cStart] : NULL;
106: (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq * dE * dE * cStart] : NULL;
107: (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq * dE * dE * cStart] : NULL;
108: (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq * cStart] : NULL;
109: (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq * cStart] : NULL;
110: (*chunkGeom)->isAffine = geom->isAffine;
111: return 0;
112: }
114: /*@C
115: PetscFEGeomRestoreChunk - Restore the chunk obtained with `PetscFEGeomCreateChunk()`
117: Input Parameters:
118: + geom - `PetscFEGeom` object
119: . cStart - The first cell in the chunk
120: . cEnd - The first cell not in the chunk
121: - chunkGeom - The chunk of cells
123: Level: intermediate
125: .seealso: `PetscFEGeom`, `PetscFEGeomGetChunk()`, `PetscFEGeomCreate()`
126: @*/
127: PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
128: {
129: PetscFree(*chunkGeom);
130: return 0;
131: }
133: /*@C
134: PetscFEGeomGetPoint - Get the geometry for cell c at point p as a `PetscFEGeom`
136: Input Parameters:
137: + geom - `PetscFEGeom` object
138: . c - The cell
139: . p - The point
140: - pcoords - The reference coordinates of point p, or NULL
142: Output Parameter:
143: . pgeom - The geometry of cell c at point p
145: Level: intermediate
147: Notes:
148: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom,
149: nothing needs to be done with it afterwards.
151: In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in.
153: .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
154: @*/
155: PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom)
156: {
157: const PetscInt dim = geom->dim;
158: const PetscInt dE = geom->dimEmbed;
159: const PetscInt Np = geom->numPoints;
162: pgeom->dim = dim;
163: pgeom->dimEmbed = dE;
164: //pgeom->isAffine = geom->isAffine;
165: if (geom->isAffine) {
166: if (!p) {
167: pgeom->xi = geom->xi;
168: pgeom->J = &geom->J[c * Np * dE * dE];
169: pgeom->invJ = &geom->invJ[c * Np * dE * dE];
170: pgeom->detJ = &geom->detJ[c * Np];
171: pgeom->n = geom->n ? &geom->n[c * Np * dE] : NULL;
172: }
173: if (pcoords) CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c * Np * dE], pgeom->J, pcoords, pgeom->v);
174: } else {
175: pgeom->v = &geom->v[(c * Np + p) * dE];
176: pgeom->J = &geom->J[(c * Np + p) * dE * dE];
177: pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
178: pgeom->detJ = &geom->detJ[c * Np + p];
179: pgeom->n = geom->n ? &geom->n[(c * Np + p) * dE] : NULL;
180: }
181: return 0;
182: }
184: /*@C
185: PetscFEGeomGetCellPoint - Get the cell geometry for face f at point p as a `PetscFEGeom`
187: Input Parameters:
188: + geom - `PetscFEGeom` object
189: . f - The face
190: - p - The point
192: Output Parameter:
193: . pgeom - The cell geometry of face f at point p
195: Level: intermediate
197: Note:
198: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom,
199: nothing needs to be done with it afterwards.
201: .seealso: `PetscFEGeom()`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
202: @*/
203: PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom)
204: {
205: const PetscBool bd = geom->dimEmbed > geom->dim && !geom->isCohesive ? PETSC_TRUE : PETSC_FALSE;
206: const PetscInt dim = bd ? geom->dimEmbed : geom->dim;
207: const PetscInt dE = geom->dimEmbed;
208: const PetscInt Np = geom->numPoints;
211: pgeom->dim = dim;
212: pgeom->dimEmbed = dE;
213: //pgeom->isAffine = geom->isAffine;
214: if (geom->isAffine) {
215: if (!p) {
216: if (bd) {
217: pgeom->J = &geom->suppJ[0][c * Np * dE * dE];
218: pgeom->invJ = &geom->suppInvJ[0][c * Np * dE * dE];
219: pgeom->detJ = &geom->suppDetJ[0][c * Np];
220: } else {
221: pgeom->J = &geom->J[c * Np * dE * dE];
222: pgeom->invJ = &geom->invJ[c * Np * dE * dE];
223: pgeom->detJ = &geom->detJ[c * Np];
224: }
225: }
226: } else {
227: if (bd) {
228: pgeom->J = &geom->suppJ[0][(c * Np + p) * dE * dE];
229: pgeom->invJ = &geom->suppInvJ[0][(c * Np + p) * dE * dE];
230: pgeom->detJ = &geom->suppDetJ[0][c * Np + p];
231: } else {
232: pgeom->J = &geom->J[(c * Np + p) * dE * dE];
233: pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
234: pgeom->detJ = &geom->detJ[c * Np + p];
235: }
236: }
237: return 0;
238: }
240: /*@
241: PetscFEGeomComplete - Calculate derived quantities from base geometry specification
243: Input Parameter:
244: . geom - `PetscFEGeom` object
246: Level: intermediate
248: .seealso: `PetscFEGeom`, `PetscFEGeomCreate()`
249: @*/
250: PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
251: {
252: PetscInt i, j, N, dE;
255: N = geom->numPoints * geom->numCells;
256: dE = geom->dimEmbed;
257: switch (dE) {
258: case 3:
259: for (i = 0; i < N; i++) {
260: DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
261: if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
262: }
263: break;
264: case 2:
265: for (i = 0; i < N; i++) {
266: DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
267: if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
268: }
269: break;
270: case 1:
271: for (i = 0; i < N; i++) {
272: geom->detJ[i] = PetscAbsReal(geom->J[i]);
273: if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
274: }
275: break;
276: }
277: if (geom->n) {
278: for (i = 0; i < N; i++) {
279: for (j = 0; j < dE; j++) geom->n[dE * i + j] = geom->J[dE * dE * i + dE * j + dE - 1] * ((dE == 2) ? -1. : 1.);
280: }
281: }
282: return 0;
283: }