Actual source code: fegeom.c

petsc-master 2019-12-03
Report Typos and Errors
  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;
 24:   PetscErrorCode  ierr;

 27:   PetscQuadratureGetData(quad,&dim,NULL,&Nq,&p,NULL);
 28:   PetscNew(&g);
 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:   PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);
 36:   if (faceData) {
 37:     PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n);
 38:     PetscCalloc6(N * dimEmbed * dimEmbed, &(g->suppJ[0]),    N * dimEmbed * dimEmbed, &(g->suppJ[1]),
 39:                         N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]),
 40:                         N,                       &(g->suppDetJ[0]), N,                       &(g->suppDetJ[1]));
 41:   }
 42:   PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);
 43:   *geom = g;
 44:   return(0);
 45: }

 47: /*@C
 48:   PetscFEGeomDestroy - Destroy a PetscFEGeom object

 50:   Input Parameter:
 51: . geom - PetscFEGeom object

 53:   Level: beginner

 55: .seealso: PetscFEGeomCreate()
 56: @*/
 57: PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom)
 58: {

 62:   if (!*geom) return(0);
 63:   PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);
 64:   PetscFree((*geom)->invJ);
 65:   PetscFree2((*geom)->face,(*geom)->n);
 66:   PetscFree6((*geom)->suppJ[0],(*geom)->suppJ[1],(*geom)->suppInvJ[0],(*geom)->suppInvJ[1],(*geom)->suppDetJ[0],(*geom)->suppDetJ[1]);
 67:   PetscFree(*geom);
 68:   return(0);
 69: }

 71: /*@C
 72:   PetscFEGeomGetChunk - Get a chunk of cells in the group as a PetscFEGeom

 74:   Input Parameters:
 75: + geom   - PetscFEGeom object
 76: . cStart - The first cell in the chunk
 77: - cEnd   - The first cell not in the chunk

 79:   Output Parameter:
 80: . chunkGeom - The chunk of cells

 82:   Level: intermediate

 84: .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate()
 85: @*/
 86: PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
 87: {
 88:   PetscInt       Nq;
 89:   PetscInt       dE;

 95:   if (!(*chunkGeom)) {
 96:     PetscNew(chunkGeom);
 97:   }
 98:   Nq = geom->numPoints;
 99:   dE= geom->dimEmbed;
100:   (*chunkGeom)->dim = geom->dim;
101:   (*chunkGeom)->dimEmbed = geom->dimEmbed;
102:   (*chunkGeom)->numPoints = geom->numPoints;
103:   (*chunkGeom)->numCells = cEnd - cStart;
104:   (*chunkGeom)->xi = geom->xi;
105:   (*chunkGeom)->v = &geom->v[Nq*dE*cStart];
106:   (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart];
107:   (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL;
108:   (*chunkGeom)->detJ = &geom->detJ[Nq*cStart];
109:   (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL;
110:   (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL;
111:   (*chunkGeom)->suppJ[0]    = geom->suppJ[0]    ? &geom->suppJ[0][Nq*dE*dE*cStart]    : NULL;
112:   (*chunkGeom)->suppJ[1]    = geom->suppJ[1]    ? &geom->suppJ[1][Nq*dE*dE*cStart]    : NULL;
113:   (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL;
114:   (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL;
115:   (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq*cStart]       : NULL;
116:   (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq*cStart]       : NULL;
117:   (*chunkGeom)->isAffine = geom->isAffine;
118:   return(0);
119: }

121: /*@C
122:   PetscFEGeomRestoreChunk - Restore the chunk

124:   Input Parameters:
125: + geom      - PetscFEGeom object
126: . cStart    - The first cell in the chunk
127: . cEnd      - The first cell not in the chunk
128: - chunkGeom - The chunk of cells

130:   Level: intermediate

132: .seealso: PetscFEGeomGetChunk(), PetscFEGeomCreate()
133: @*/
134: PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
135: {

139:   PetscFree(*chunkGeom);
140:   return(0);
141: }

143: /*@
144:   PetscFEGeomComplete - Calculate derived quntites from base geometry specification

146:   Input Parameter:
147: . geom - PetscFEGeom object

149:   Level: intermediate

151: .seealso: PetscFEGeomCreate()
152: @*/
153: PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
154: {
155:   PetscInt i, j, N, dE;

158:   N = geom->numPoints * geom->numCells;
159:   dE = geom->dimEmbed;
160:   switch (dE) {
161:   case 3:
162:     for (i = 0; i < N; i++) {
163:       DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
164:       if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
165:     }
166:     break;
167:   case 2:
168:     for (i = 0; i < N; i++) {
169:       DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
170:       if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
171:     }
172:     break;
173:   case 1:
174:     for (i = 0; i < N; i++) {
175:       geom->detJ[i] = PetscAbsReal(geom->J[i]);
176:       if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
177:     }
178:     break;
179:   }
180:   if (geom->n) {
181:     for (i = 0; i < N; i++) {
182:       for (j = 0; j < dE; j++) {
183:         geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.);
184:       }
185:     }
186:   }
187:   return(0);
188: }