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