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: - mode - Type of geometry data to store
12: Output Parameter:
13: . geom - The `PetscFEGeom` object, which is a struct not a `PetscObject`
15: Level: beginner
17: .seealso: `PetscFEGeom`, `PetscQuadrature`, `PetscFEGeomDestroy()`, `PetscFEGeomComplete()`
18: @*/
19: PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscFEGeomMode mode, PetscFEGeom **geom)
20: {
21: PetscFEGeom *g;
22: PetscInt dim, Nq, N;
23: const PetscReal *p;
25: PetscFunctionBegin;
26: PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &Nq, &p, NULL));
27: PetscCall(PetscNew(&g));
28: g->mode = mode;
29: g->xi = p;
30: g->numCells = numCells;
31: g->numPoints = Nq;
32: g->dim = dim;
33: g->dimEmbed = dimEmbed;
34: N = numCells * Nq;
35: PetscCall(PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ));
36: if (mode == PETSC_FEGEOM_BOUNDARY || mode == PETSC_FEGEOM_COHESIVE) {
37: PetscCall(PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n));
38: PetscCall(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]));
39: }
40: PetscCall(PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ));
41: *geom = g;
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: /*@C
46: PetscFEGeomDestroy - Destroy a `PetscFEGeom` object
48: Input Parameter:
49: . geom - `PetscFEGeom` object
51: Level: beginner
53: .seealso: `PetscFEGeom`, `PetscFEGeomCreate()`
54: @*/
55: PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom)
56: {
57: PetscFunctionBegin;
58: if (!*geom) PetscFunctionReturn(PETSC_SUCCESS);
59: PetscCall(PetscFree3((*geom)->v, (*geom)->J, (*geom)->detJ));
60: PetscCall(PetscFree((*geom)->invJ));
61: PetscCall(PetscFree2((*geom)->face, (*geom)->n));
62: PetscCall(PetscFree6((*geom)->suppJ[0], (*geom)->suppJ[1], (*geom)->suppInvJ[0], (*geom)->suppInvJ[1], (*geom)->suppDetJ[0], (*geom)->suppDetJ[1]));
63: PetscCall(PetscFree(*geom));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*@C
68: PetscFEGeomGetChunk - Get a chunk of cells in the group as a `PetscFEGeom`
70: Input Parameters:
71: + geom - `PetscFEGeom` object
72: . cStart - The first cell in the chunk
73: - cEnd - The first cell not in the chunk
75: Output Parameter:
76: . chunkGeom - an array of cells of length `cEnd` - `cStart`
78: Level: intermediate
80: Note:
81: Use `PetscFEGeomRestoreChunk()` to return the result
83: .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
84: @*/
85: PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom *chunkGeom[])
86: {
87: PetscInt Nq;
88: PetscInt dE;
90: PetscFunctionBegin;
91: PetscAssertPointer(geom, 1);
92: PetscAssertPointer(chunkGeom, 4);
93: if (!*chunkGeom) PetscCall(PetscNew(chunkGeom));
94: Nq = geom->numPoints;
95: dE = geom->dimEmbed;
96: (*chunkGeom)->mode = geom->mode;
97: (*chunkGeom)->dim = geom->dim;
98: (*chunkGeom)->dimEmbed = geom->dimEmbed;
99: (*chunkGeom)->numPoints = geom->numPoints;
100: (*chunkGeom)->numCells = cEnd - cStart;
101: (*chunkGeom)->xi = geom->xi;
102: (*chunkGeom)->v = PetscSafePointerPlusOffset(geom->v, Nq * dE * cStart);
103: (*chunkGeom)->J = PetscSafePointerPlusOffset(geom->J, Nq * dE * dE * cStart);
104: (*chunkGeom)->invJ = PetscSafePointerPlusOffset(geom->invJ, Nq * dE * dE * cStart);
105: (*chunkGeom)->detJ = PetscSafePointerPlusOffset(geom->detJ, Nq * cStart);
106: (*chunkGeom)->n = PetscSafePointerPlusOffset(geom->n, Nq * dE * cStart);
107: (*chunkGeom)->face = PetscSafePointerPlusOffset(geom->face, cStart);
108: (*chunkGeom)->suppJ[0] = PetscSafePointerPlusOffset(geom->suppJ[0], Nq * dE * dE * cStart);
109: (*chunkGeom)->suppJ[1] = PetscSafePointerPlusOffset(geom->suppJ[1], Nq * dE * dE * cStart);
110: (*chunkGeom)->suppInvJ[0] = PetscSafePointerPlusOffset(geom->suppInvJ[0], Nq * dE * dE * cStart);
111: (*chunkGeom)->suppInvJ[1] = PetscSafePointerPlusOffset(geom->suppInvJ[1], Nq * dE * dE * cStart);
112: (*chunkGeom)->suppDetJ[0] = PetscSafePointerPlusOffset(geom->suppDetJ[0], Nq * cStart);
113: (*chunkGeom)->suppDetJ[1] = PetscSafePointerPlusOffset(geom->suppDetJ[1], Nq * cStart);
114: (*chunkGeom)->isAffine = geom->isAffine;
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: /*@C
119: PetscFEGeomRestoreChunk - Restore the chunk obtained with `PetscFEGeomCreateChunk()`
121: Input Parameters:
122: + geom - `PetscFEGeom` object
123: . cStart - The first cell in the chunk
124: . cEnd - The first cell not in the chunk
125: - chunkGeom - The chunk of cells
127: Level: intermediate
129: .seealso: `PetscFEGeom`, `PetscFEGeomGetChunk()`, `PetscFEGeomCreate()`
130: @*/
131: PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
132: {
133: PetscFunctionBegin;
134: PetscCall(PetscFree(*chunkGeom));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: /*@C
139: PetscFEGeomGetPoint - Get the geometry for cell `c` at point `p` as a `PetscFEGeom`
141: Input Parameters:
142: + geom - `PetscFEGeom` object
143: . c - The cell
144: . p - The point
145: - pcoords - The reference coordinates of point `p`, or `NULL`
147: Output Parameter:
148: . pgeom - The geometry of cell `c` at point `p`
150: Level: intermediate
152: Notes:
153: For affine geometries, this only copies to `pgeom` at point 0. Since we copy pointers into `pgeom`,
154: nothing needs to be done with it afterwards.
156: In the affine case, `pgeom` must have storage for the integration point coordinates in pgeom->v if `pcoords` is passed in.
158: .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
159: @*/
160: PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom)
161: {
162: const PetscInt dim = geom->dim;
163: const PetscInt dE = geom->dimEmbed;
164: const PetscInt Np = geom->numPoints;
166: PetscFunctionBeginHot;
167: pgeom->mode = geom->mode;
168: pgeom->dim = dim;
169: pgeom->dimEmbed = dE;
170: //pgeom->isAffine = geom->isAffine;
171: if (geom->isAffine) {
172: if (!p) {
173: pgeom->xi = geom->xi;
174: pgeom->J = &geom->J[c * Np * dE * dE];
175: pgeom->invJ = &geom->invJ[c * Np * dE * dE];
176: pgeom->detJ = &geom->detJ[c * Np];
177: pgeom->n = PetscSafePointerPlusOffset(geom->n, c * Np * dE);
178: }
179: if (pcoords) CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c * Np * dE], pgeom->J, pcoords, pgeom->v);
180: } else {
181: pgeom->v = &geom->v[(c * Np + p) * dE];
182: pgeom->J = &geom->J[(c * Np + p) * dE * dE];
183: pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
184: pgeom->detJ = &geom->detJ[c * Np + p];
185: pgeom->n = PetscSafePointerPlusOffset(geom->n, (c * Np + p) * dE);
186: }
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: /*@C
191: PetscFEGeomGetCellPoint - Get the cell geometry for cell `c` at point `p` as a `PetscFEGeom`
193: Input Parameters:
194: + geom - `PetscFEGeom` object
195: . c - The cell
196: - p - The point
198: Output Parameter:
199: . pgeom - The cell geometry of cell `c` at point `p`
201: Level: intermediate
203: Notes:
204: For PETSC_FEGEOM_BOUNDARY mode, this gives the geometry for supporting cell 0. For PETSC_FEGEOM_COHESIVE mode,
205: this gives the bulk geometry for that internal face.
207: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into `pgeom`,
208: nothing needs to be done with it afterwards.
210: .seealso: `PetscFEGeom`, `PetscFEGeomMode`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
211: @*/
212: PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom)
213: {
214: const PetscBool bd = geom->mode == PETSC_FEGEOM_BOUNDARY ? PETSC_TRUE : PETSC_FALSE;
215: const PetscInt dim = bd ? geom->dimEmbed : geom->dim;
216: const PetscInt dE = geom->dimEmbed;
217: const PetscInt Np = geom->numPoints;
219: PetscFunctionBeginHot;
220: pgeom->mode = geom->mode;
221: pgeom->dim = dim;
222: pgeom->dimEmbed = dE;
223: //pgeom->isAffine = geom->isAffine;
224: if (geom->isAffine) {
225: if (!p) {
226: if (bd) {
227: pgeom->J = &geom->suppJ[0][c * Np * dE * dE];
228: pgeom->invJ = &geom->suppInvJ[0][c * Np * dE * dE];
229: pgeom->detJ = &geom->suppDetJ[0][c * Np];
230: } else {
231: pgeom->J = &geom->J[c * Np * dE * dE];
232: pgeom->invJ = &geom->invJ[c * Np * dE * dE];
233: pgeom->detJ = &geom->detJ[c * Np];
234: }
235: }
236: } else {
237: if (bd) {
238: pgeom->J = &geom->suppJ[0][(c * Np + p) * dE * dE];
239: pgeom->invJ = &geom->suppInvJ[0][(c * Np + p) * dE * dE];
240: pgeom->detJ = &geom->suppDetJ[0][c * Np + p];
241: } else {
242: pgeom->J = &geom->J[(c * Np + p) * dE * dE];
243: pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
244: pgeom->detJ = &geom->detJ[c * Np + p];
245: }
246: }
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: /*@C
251: PetscFEGeomComplete - Calculate derived quantities from a base geometry specification
253: Input Parameter:
254: . geom - `PetscFEGeom` object
256: Level: intermediate
258: .seealso: `PetscFEGeom`, `PetscFEGeomCreate()`
259: @*/
260: PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
261: {
262: PetscInt i, j, N, dE;
264: PetscFunctionBeginHot;
265: N = geom->numPoints * geom->numCells;
266: dE = geom->dimEmbed;
267: switch (dE) {
268: case 3:
269: for (i = 0; i < N; i++) {
270: DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
271: if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
272: }
273: break;
274: case 2:
275: for (i = 0; i < N; i++) {
276: DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
277: if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
278: }
279: break;
280: case 1:
281: for (i = 0; i < N; i++) {
282: geom->detJ[i] = PetscAbsReal(geom->J[i]);
283: if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
284: }
285: break;
286: }
287: if (geom->n) {
288: for (i = 0; i < N; i++) {
289: for (j = 0; j < dE; j++) geom->n[dE * i + j] = geom->J[dE * dE * i + dE * j + dE - 1] * ((dE == 2) ? -1. : 1.);
290: }
291: }
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }