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`, `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 - The chunk of cells

 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 f 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 f 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 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: }