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:   } else {
 23:     PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);
 24:   }
 25:   *geom = g;
 26:   return(0);
 27: }

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

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

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

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

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

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

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

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