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: `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: `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: .seealso: `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
 79: @*/
 80: PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
 81: {
 82:   PetscInt Nq;
 83:   PetscInt dE;

 87:   if (!(*chunkGeom)) PetscNew(chunkGeom);
 88:   Nq                        = geom->numPoints;
 89:   dE                        = geom->dimEmbed;
 90:   (*chunkGeom)->dim         = geom->dim;
 91:   (*chunkGeom)->dimEmbed    = geom->dimEmbed;
 92:   (*chunkGeom)->numPoints   = geom->numPoints;
 93:   (*chunkGeom)->numCells    = cEnd - cStart;
 94:   (*chunkGeom)->xi          = geom->xi;
 95:   (*chunkGeom)->v           = &geom->v[Nq * dE * cStart];
 96:   (*chunkGeom)->J           = &geom->J[Nq * dE * dE * cStart];
 97:   (*chunkGeom)->invJ        = (geom->invJ) ? &geom->invJ[Nq * dE * dE * cStart] : NULL;
 98:   (*chunkGeom)->detJ        = &geom->detJ[Nq * cStart];
 99:   (*chunkGeom)->n           = geom->n ? &geom->n[Nq * dE * cStart] : NULL;
100:   (*chunkGeom)->face        = geom->face ? &geom->face[cStart] : NULL;
101:   (*chunkGeom)->suppJ[0]    = geom->suppJ[0] ? &geom->suppJ[0][Nq * dE * dE * cStart] : NULL;
102:   (*chunkGeom)->suppJ[1]    = geom->suppJ[1] ? &geom->suppJ[1][Nq * dE * dE * cStart] : NULL;
103:   (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq * dE * dE * cStart] : NULL;
104:   (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq * dE * dE * cStart] : NULL;
105:   (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq * cStart] : NULL;
106:   (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq * cStart] : NULL;
107:   (*chunkGeom)->isAffine    = geom->isAffine;
108:   return 0;
109: }

111: /*@C
112:   PetscFEGeomRestoreChunk - Restore the chunk

114:   Input Parameters:
115: + geom      - PetscFEGeom object
116: . cStart    - The first cell in the chunk
117: . cEnd      - The first cell not in the chunk
118: - chunkGeom - The chunk of cells

120:   Level: intermediate

122: .seealso: `PetscFEGeomGetChunk()`, `PetscFEGeomCreate()`
123: @*/
124: PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
125: {
126:   PetscFree(*chunkGeom);
127:   return 0;
128: }

130: /*@C
131:   PetscFEGeomGetPoint - Get the geometry for cell c at point p as a PetscFEGeom

133:   Input Parameters:
134: + geom    - PetscFEGeom object
135: . c       - The cell
136: . p       - The point
137: - pcoords - The reference coordinates of point p, or NULL

139:   Output Parameter:
140: . pgeom - The geometry of cell c at point p

142:   Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom,
143:   nothing needs to be done with it afterwards.

145:   In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in.

147:   Level: intermediate

149: .seealso: `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
150: @*/
151: PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom)
152: {
153:   const PetscInt dim = geom->dim;
154:   const PetscInt dE  = geom->dimEmbed;
155:   const PetscInt Np  = geom->numPoints;

158:   pgeom->dim      = dim;
159:   pgeom->dimEmbed = dE;
160:   //pgeom->isAffine = geom->isAffine;
161:   if (geom->isAffine) {
162:     if (!p) {
163:       pgeom->xi   = geom->xi;
164:       pgeom->J    = &geom->J[c * Np * dE * dE];
165:       pgeom->invJ = &geom->invJ[c * Np * dE * dE];
166:       pgeom->detJ = &geom->detJ[c * Np];
167:       pgeom->n    = geom->n ? &geom->n[c * Np * dE] : NULL;
168:     }
169:     if (pcoords) CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c * Np * dE], pgeom->J, pcoords, pgeom->v);
170:   } else {
171:     pgeom->v    = &geom->v[(c * Np + p) * dE];
172:     pgeom->J    = &geom->J[(c * Np + p) * dE * dE];
173:     pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
174:     pgeom->detJ = &geom->detJ[c * Np + p];
175:     pgeom->n    = geom->n ? &geom->n[(c * Np + p) * dE] : NULL;
176:   }
177:   return 0;
178: }

180: /*@C
181:   PetscFEGeomGetCellPoint - Get the cell geometry for face f at point p as a PetscFEGeom

183:   Input Parameters:
184: + geom    - PetscFEGeom object
185: . f       - The face
186: - p       - The point

188:   Output Parameter:
189: . pgeom - The cell geometry of face f at point p

191:   Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom,
192:   nothing needs to be done with it afterwards.

194:   Level: intermediate

196: .seealso: `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
197: @*/
198: PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom)
199: {
200:   const PetscBool bd  = geom->dimEmbed > geom->dim && !geom->isCohesive ? PETSC_TRUE : PETSC_FALSE;
201:   const PetscInt  dim = bd ? geom->dimEmbed : geom->dim;
202:   const PetscInt  dE  = geom->dimEmbed;
203:   const PetscInt  Np  = geom->numPoints;

206:   pgeom->dim      = dim;
207:   pgeom->dimEmbed = dE;
208:   //pgeom->isAffine = geom->isAffine;
209:   if (geom->isAffine) {
210:     if (!p) {
211:       if (bd) {
212:         pgeom->J    = &geom->suppJ[0][c * Np * dE * dE];
213:         pgeom->invJ = &geom->suppInvJ[0][c * Np * dE * dE];
214:         pgeom->detJ = &geom->suppDetJ[0][c * Np];
215:       } else {
216:         pgeom->J    = &geom->J[c * Np * dE * dE];
217:         pgeom->invJ = &geom->invJ[c * Np * dE * dE];
218:         pgeom->detJ = &geom->detJ[c * Np];
219:       }
220:     }
221:   } else {
222:     if (bd) {
223:       pgeom->J    = &geom->suppJ[0][(c * Np + p) * dE * dE];
224:       pgeom->invJ = &geom->suppInvJ[0][(c * Np + p) * dE * dE];
225:       pgeom->detJ = &geom->suppDetJ[0][c * Np + p];
226:     } else {
227:       pgeom->J    = &geom->J[(c * Np + p) * dE * dE];
228:       pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
229:       pgeom->detJ = &geom->detJ[c * Np + p];
230:     }
231:   }
232:   return 0;
233: }

235: /*@
236:   PetscFEGeomComplete - Calculate derived quntites from base geometry specification

238:   Input Parameter:
239: . geom - PetscFEGeom object

241:   Level: intermediate

243: .seealso: `PetscFEGeomCreate()`
244: @*/
245: PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
246: {
247:   PetscInt i, j, N, dE;

250:   N  = geom->numPoints * geom->numCells;
251:   dE = geom->dimEmbed;
252:   switch (dE) {
253:   case 3:
254:     for (i = 0; i < N; i++) {
255:       DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
256:       if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
257:     }
258:     break;
259:   case 2:
260:     for (i = 0; i < N; i++) {
261:       DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
262:       if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
263:     }
264:     break;
265:   case 1:
266:     for (i = 0; i < N; i++) {
267:       geom->detJ[i] = PetscAbsReal(geom->J[i]);
268:       if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
269:     }
270:     break;
271:   }
272:   if (geom->n) {
273:     for (i = 0; i < N; i++) {
274:       for (j = 0; j < dE; j++) geom->n[dE * i + j] = geom->J[dE * dE * i + dE * j + dE - 1] * ((dE == 2) ? -1. : 1.);
275:     }
276:   }
277:   return 0;
278: }