1:  #include <petsc/private/petscfeimpl.h>

  3: PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscBool faceData, PetscFEGeom **geom)
  4: {
  5:   PetscFEGeom     *g;
  6:   PetscInt        dim, Nq, N;
  7:   const PetscReal *p;
  8:   PetscErrorCode  ierr;

 11:   PetscQuadratureGetData(quad,&dim,NULL,&Nq,&p,NULL);
 12:   PetscNew(&g);
 13:   g->xi        = p;
 14:   g->numCells  = numCells;
 15:   g->numPoints = Nq;
 16:   g->dim       = dim;
 17:   g->dimEmbed  = dimEmbed;
 18:   N = numCells * Nq;
 19:   PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);
 20:   if (faceData) {
 21:     PetscCalloc4(numCells, &g->face, N * dimEmbed, &g->n, N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]));
 22:   }
 23:   PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);
 24:   *geom = g;
 25:   return(0);
 26: }

 28: PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom)
 29: {

 33:   if (!*geom) return(0);
 34:   PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);
 35:   PetscFree((*geom)->invJ);
 36:   PetscFree4((*geom)->face,(*geom)->n,(*geom)->suppInvJ[0],(*geom)->suppInvJ[1]);
 37:   PetscFree(*geom);
 38:   return(0);
 39: }

 41: PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
 42: {
 43:   PetscInt       Nq;
 44:   PetscInt       dE;

 50:   if (!(*chunkGeom)) {
 51:     PetscNew(chunkGeom);
 52:   }
 53:   Nq = geom->numPoints;
 54:   dE= geom->dimEmbed;
 55:   (*chunkGeom)->dim = geom->dim;
 56:   (*chunkGeom)->dimEmbed = geom->dimEmbed;
 57:   (*chunkGeom)->numPoints = geom->numPoints;
 58:   (*chunkGeom)->numCells = cEnd - cStart;
 59:   (*chunkGeom)->xi = geom->xi;
 60:   (*chunkGeom)->v = &geom->v[Nq*dE*cStart];
 61:   (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart];
 62:   (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL;
 63:   (*chunkGeom)->detJ = &geom->detJ[Nq*cStart];
 64:   (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL;
 65:   (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL;
 66:   (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL;
 67:   (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL;
 68:   (*chunkGeom)->isAffine = geom->isAffine;
 69:   return(0);
 70: }

 72: PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
 73: {

 77:   PetscFree(*chunkGeom);
 78:   return(0);
 79: }

 81: PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
 82: {
 83:   PetscInt i, j, N, dE;

 86:   N = geom->numPoints * geom->numCells;
 87:   dE = geom->dimEmbed;
 88:   switch (dE) {
 89:   case 3:
 90:     for (i = 0; i < N; i++) {
 91:       DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
 92:       if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
 93:     }
 94:     break;
 95:   case 2:
 96:     for (i = 0; i < N; i++) {
 97:       DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
 98:       if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
 99:     }
100:     break;
101:   case 1:
102:     for (i = 0; i < N; i++) {
103:       geom->detJ[i] = PetscAbsReal(geom->J[i]);
104:       if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
105:     }
106:     break;
107:   }
108:   if (geom->n) {
109:     for (i = 0; i < N; i++) {
110:       for (j = 0; j < dE; j++) {
111:         geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.);
112:       }
113:     }
114:   }
115:   return(0);
116: }