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