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, 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, PetscBool faceData, 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->xi = p;
29: g->numCells = numCells;
30: g->numPoints = Nq;
31: g->dim = dim;
32: g->dimEmbed = dimEmbed;
33: g->isCohesive = PETSC_FALSE;
34: N = numCells * Nq;
35: PetscCall(PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ));
36: if (faceData) {
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)->dim = geom->dim;
97: (*chunkGeom)->dimEmbed = geom->dimEmbed;
98: (*chunkGeom)->numPoints = geom->numPoints;
99: (*chunkGeom)->numCells = cEnd - cStart;
100: (*chunkGeom)->xi = geom->xi;
101: (*chunkGeom)->v = PetscSafePointerPlusOffset(geom->v, Nq * dE * cStart);
102: (*chunkGeom)->J = PetscSafePointerPlusOffset(geom->J, Nq * dE * dE * cStart);
103: (*chunkGeom)->invJ = PetscSafePointerPlusOffset(geom->invJ, Nq * dE * dE * cStart);
104: (*chunkGeom)->detJ = PetscSafePointerPlusOffset(geom->detJ, Nq * cStart);
105: (*chunkGeom)->n = PetscSafePointerPlusOffset(geom->n, Nq * dE * cStart);
106: (*chunkGeom)->face = PetscSafePointerPlusOffset(geom->face, cStart);
107: (*chunkGeom)->suppJ[0] = PetscSafePointerPlusOffset(geom->suppJ[0], Nq * dE * dE * cStart);
108: (*chunkGeom)->suppJ[1] = PetscSafePointerPlusOffset(geom->suppJ[1], Nq * dE * dE * cStart);
109: (*chunkGeom)->suppInvJ[0] = PetscSafePointerPlusOffset(geom->suppInvJ[0], Nq * dE * dE * cStart);
110: (*chunkGeom)->suppInvJ[1] = PetscSafePointerPlusOffset(geom->suppInvJ[1], Nq * dE * dE * cStart);
111: (*chunkGeom)->suppDetJ[0] = PetscSafePointerPlusOffset(geom->suppDetJ[0], Nq * cStart);
112: (*chunkGeom)->suppDetJ[1] = PetscSafePointerPlusOffset(geom->suppDetJ[1], Nq * cStart);
113: (*chunkGeom)->isAffine = geom->isAffine;
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: /*@C
118: PetscFEGeomRestoreChunk - Restore the chunk obtained with `PetscFEGeomCreateChunk()`
120: Input Parameters:
121: + geom - `PetscFEGeom` object
122: . cStart - The first cell in the chunk
123: . cEnd - The first cell not in the chunk
124: - chunkGeom - The chunk of cells
126: Level: intermediate
128: .seealso: `PetscFEGeom`, `PetscFEGeomGetChunk()`, `PetscFEGeomCreate()`
129: @*/
130: PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
131: {
132: PetscFunctionBegin;
133: PetscCall(PetscFree(*chunkGeom));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: /*@C
138: PetscFEGeomGetPoint - Get the geometry for cell `c` at point `p` as a `PetscFEGeom`
140: Input Parameters:
141: + geom - `PetscFEGeom` object
142: . c - The cell
143: . p - The point
144: - pcoords - The reference coordinates of point `p`, or `NULL`
146: Output Parameter:
147: . pgeom - The geometry of cell `c` at point `p`
149: Level: intermediate
151: Notes:
152: For affine geometries, this only copies to `pgeom` at point 0. Since we copy pointers into `pgeom`,
153: nothing needs to be done with it afterwards.
155: In the affine case, `pgeom` must have storage for the integration point coordinates in pgeom->v if `pcoords` is passed in.
157: .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
158: @*/
159: PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom)
160: {
161: const PetscInt dim = geom->dim;
162: const PetscInt dE = geom->dimEmbed;
163: const PetscInt Np = geom->numPoints;
165: PetscFunctionBeginHot;
166: pgeom->dim = dim;
167: pgeom->dimEmbed = dE;
168: //pgeom->isAffine = geom->isAffine;
169: if (geom->isAffine) {
170: if (!p) {
171: pgeom->xi = geom->xi;
172: pgeom->J = &geom->J[c * Np * dE * dE];
173: pgeom->invJ = &geom->invJ[c * Np * dE * dE];
174: pgeom->detJ = &geom->detJ[c * Np];
175: pgeom->n = PetscSafePointerPlusOffset(geom->n, c * Np * dE);
176: }
177: if (pcoords) CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c * Np * dE], pgeom->J, pcoords, pgeom->v);
178: } else {
179: pgeom->v = &geom->v[(c * Np + p) * dE];
180: pgeom->J = &geom->J[(c * Np + p) * dE * dE];
181: pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
182: pgeom->detJ = &geom->detJ[c * Np + p];
183: pgeom->n = PetscSafePointerPlusOffset(geom->n, (c * Np + p) * dE);
184: }
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: /*@C
189: PetscFEGeomGetCellPoint - Get the cell geometry for face `c` at point `p` as a `PetscFEGeom`
191: Input Parameters:
192: + geom - `PetscFEGeom` object
193: . c - The face
194: - p - The point
196: Output Parameter:
197: . pgeom - The cell geometry of face `c` at point `p`
199: Level: intermediate
201: Note:
202: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into `pgeom`,
203: nothing needs to be done with it afterwards.
205: .seealso: `PetscFEGeom()`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
206: @*/
207: PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom)
208: {
209: const PetscBool bd = geom->dimEmbed > geom->dim && !geom->isCohesive ? PETSC_TRUE : PETSC_FALSE;
210: const PetscInt dim = bd ? geom->dimEmbed : geom->dim;
211: const PetscInt dE = geom->dimEmbed;
212: const PetscInt Np = geom->numPoints;
214: PetscFunctionBeginHot;
215: pgeom->dim = dim;
216: pgeom->dimEmbed = dE;
217: //pgeom->isAffine = geom->isAffine;
218: if (geom->isAffine) {
219: if (!p) {
220: if (bd) {
221: pgeom->J = &geom->suppJ[0][c * Np * dE * dE];
222: pgeom->invJ = &geom->suppInvJ[0][c * Np * dE * dE];
223: pgeom->detJ = &geom->suppDetJ[0][c * Np];
224: } else {
225: pgeom->J = &geom->J[c * Np * dE * dE];
226: pgeom->invJ = &geom->invJ[c * Np * dE * dE];
227: pgeom->detJ = &geom->detJ[c * Np];
228: }
229: }
230: } else {
231: if (bd) {
232: pgeom->J = &geom->suppJ[0][(c * Np + p) * dE * dE];
233: pgeom->invJ = &geom->suppInvJ[0][(c * Np + p) * dE * dE];
234: pgeom->detJ = &geom->suppDetJ[0][c * Np + p];
235: } else {
236: pgeom->J = &geom->J[(c * Np + p) * dE * dE];
237: pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
238: pgeom->detJ = &geom->detJ[c * Np + p];
239: }
240: }
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: /*@C
245: PetscFEGeomComplete - Calculate derived quantities from a base geometry specification
247: Input Parameter:
248: . geom - `PetscFEGeom` object
250: Level: intermediate
252: .seealso: `PetscFEGeom`, `PetscFEGeomCreate()`
253: @*/
254: PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
255: {
256: PetscInt i, j, N, dE;
258: PetscFunctionBeginHot;
259: N = geom->numPoints * geom->numCells;
260: dE = geom->dimEmbed;
261: switch (dE) {
262: case 3:
263: for (i = 0; i < N; i++) {
264: DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
265: if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
266: }
267: break;
268: case 2:
269: for (i = 0; i < N; i++) {
270: DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
271: if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
272: }
273: break;
274: case 1:
275: for (i = 0; i < N; i++) {
276: geom->detJ[i] = PetscAbsReal(geom->J[i]);
277: if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
278: }
279: break;
280: }
281: if (geom->n) {
282: for (i = 0; i < N; i++) {
283: for (j = 0; j < dE; j++) geom->n[dE * i + j] = geom->J[dE * dE * i + dE * j + dE - 1] * ((dE == 2) ? -1. : 1.);
284: }
285: }
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }