Actual source code: plexceed.c

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

      3: PetscErrorCode DMGetPoints_Internal(DM dm, DMLabel domainLabel, PetscInt labelVal, PetscInt height, IS *pointIS)
      4: {
      5:   PetscInt depth;
      6:   DMLabel  depthLabel;
      7:   IS       depthIS;

      9:   PetscFunctionBegin;
     10:   PetscCall(DMPlexGetDepth(dm, &depth));
     11:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
     12:   PetscCall(DMLabelGetStratumIS(depthLabel, depth - height, &depthIS));
     13:   if (domainLabel) {
     14:     IS domainIS;

     16:     PetscCall(DMLabelGetStratumIS(domainLabel, labelVal, &domainIS));
     17:     if (domainIS) { // domainIS is non-empty
     18:       PetscCall(ISIntersect(depthIS, domainIS, pointIS));
     19:       PetscCall(ISDestroy(&domainIS));
     20:     } else { // domainIS is NULL (empty)
     21:       *pointIS = NULL;
     22:     }
     23:     PetscCall(ISDestroy(&depthIS));
     24:   } else {
     25:     *pointIS = depthIS;
     26:   }
     27:   PetscFunctionReturn(PETSC_SUCCESS);
     28: }

     30: /*@C
     31:   DMPlexGetLocalOffsets - Allocate and populate array of local offsets for each cell closure.

     33:   Not collective

     35:   Input Parameters:
     36: + dm           - The `DMPLEX` object
     37: . domain_label - label for `DMPLEX` domain, or NULL for whole domain
     38: . label_value  - Stratum value
     39: . height       - Height of target cells in `DMPLEX` topology
     40: - dm_field     - Index of `DMPLEX` field

     42:   Output Parameters:
     43: + num_cells - Number of local cells
     44: . cell_size - Size of each cell, given by cell_size * num_comp = num_dof
     45: . num_comp  - Number of components per dof
     46: . l_size    - Size of local vector
     47: - offsets   - Allocated offsets array for cells

     49:   Level: developer

     51:   Notes:
     52:   Allocate and populate array of shape [num_cells, cell_size] defining offsets for each value (cell, node) for local vector of the `DMPLEX` field. All offsets are in the range [0, l_size - 1].

     54:   Caller is responsible for freeing the offsets array using `PetscFree()`.

     56: .seealso: [](ch_unstructured), `DMPlexGetLocalOffsetsSupport()`, `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
     57: @*/
     58: PetscErrorCode DMPlexGetLocalOffsets(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, PetscInt *num_cells, PetscInt *cell_size, PetscInt *num_comp, PetscInt *l_size, PetscInt *offsets[])
     59: {
     60:   PetscDS         ds = NULL;
     61:   PetscFE         fe;
     62:   PetscSection    section;
     63:   PetscInt        dim, ds_field = -1;
     64:   PetscInt       *restr_indices;
     65:   const PetscInt *iter_indices;
     66:   IS              iter_is;

     68:   PetscFunctionBeginUser;
     70:   PetscCall(PetscLogEventBegin(DMPLEX_GetLocalOffsets, dm, 0, 0, 0));
     71:   PetscCall(DMGetDimension(dm, &dim));
     72:   PetscCall(DMGetLocalSection(dm, &section));
     73:   PetscCall(PetscSectionGetStorageSize(section, l_size));
     74:   {
     75:     IS              field_is;
     76:     const PetscInt *fields;
     77:     PetscInt        num_fields;

     79:     PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds, NULL));
     80:     PetscCheck(field_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Domain label does not have any fields associated with it");
     81:     // Translate dm_field to ds_field
     82:     PetscCall(ISGetIndices(field_is, &fields));
     83:     PetscCall(ISGetSize(field_is, &num_fields));
     84:     for (PetscInt i = 0; i < num_fields; i++) {
     85:       if (dm_field == fields[i]) {
     86:         ds_field = i;
     87:         break;
     88:       }
     89:     }
     90:     PetscCall(ISRestoreIndices(field_is, &fields));
     91:   }
     92:   PetscCheck(ds_field != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field);

     94:   PetscCall(DMGetPoints_Internal(dm, domain_label, label_value, height, &iter_is));
     95:   if (iter_is) {
     96:     PetscCall(ISGetLocalSize(iter_is, num_cells));
     97:     PetscCall(ISGetIndices(iter_is, &iter_indices));
     98:   } else {
     99:     *num_cells   = 0;
    100:     iter_indices = NULL;
    101:   }

    103:   {
    104:     PetscDualSpace dual_space;
    105:     PetscInt       num_dual_basis_vectors;

    107:     PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe));
    108:     PetscCall(PetscFEGetHeightSubspace(fe, height, &fe));
    109:     PetscCheck(fe, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Height %" PetscInt_FMT " is invalid for DG discretizations", height);
    110:     PetscCall(PetscFEGetDualSpace(fe, &dual_space));
    111:     PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors));
    112:     PetscCall(PetscDualSpaceGetNumComponents(dual_space, num_comp));
    113:     PetscCheck(num_dual_basis_vectors % *num_comp == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for number of dual basis vectors %" PetscInt_FMT " not divisible by %" PetscInt_FMT " components", num_dual_basis_vectors, *num_comp);
    114:     *cell_size = num_dual_basis_vectors / *num_comp;
    115:   }
    116:   PetscInt restr_size = (*num_cells) * (*cell_size);
    117:   PetscCall(PetscMalloc1(restr_size, &restr_indices));
    118:   PetscInt cell_offset = 0;

    120:   PetscInt P = dim - height ? (PetscInt)PetscPowReal(*cell_size, 1.0 / (dim - height)) : 0;
    121:   for (PetscInt p = 0; p < *num_cells; p++) {
    122:     PetscBool flip = PETSC_FALSE;
    123:     PetscInt  c    = iter_indices[p];
    124:     PetscInt  num_indices, *indices;
    125:     PetscInt  field_offsets[17]; // max number of fields plus 1
    126:     PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
    127:     if (height > 0) {
    128:       PetscInt        num_cells_support, num_faces, start = -1;
    129:       const PetscInt *orients, *faces, *cells;
    130:       PetscCall(DMPlexGetSupport(dm, c, &cells));
    131:       PetscCall(DMPlexGetSupportSize(dm, c, &num_cells_support));
    132:       PetscCheck(num_cells_support == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Expected one cell in support of exterior face, but got %" PetscInt_FMT " cells", num_cells_support);
    133:       PetscCall(DMPlexGetCone(dm, cells[0], &faces));
    134:       PetscCall(DMPlexGetConeSize(dm, cells[0], &num_faces));
    135:       for (PetscInt i = 0; i < num_faces; i++) {
    136:         if (faces[i] == c) start = i;
    137:       }
    138:       PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Could not find face %" PetscInt_FMT " in cone of its support", c);
    139:       PetscCall(DMPlexGetConeOrientation(dm, cells[0], &orients));
    140:       if (orients[start] < 0) flip = PETSC_TRUE;
    141:     }

    143:     for (PetscInt i = 0; i < *cell_size; i++) {
    144:       PetscInt ii = i;
    145:       if (flip) {
    146:         if (*cell_size == P) ii = *cell_size - 1 - i;
    147:         else if (*cell_size == P * P) {
    148:           PetscInt row = i / P, col = i % P;
    149:           ii = row + col * P;
    150:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for flipping point with cell size %" PetscInt_FMT " != P (%" PetscInt_FMT ") or P^2", *cell_size, P);
    151:       }
    152:       // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
    153:       PetscInt loc                 = indices[field_offsets[dm_field] + ii * (*num_comp)];
    154:       loc                          = loc < 0 ? -(loc + 1) : loc;
    155:       restr_indices[cell_offset++] = loc;
    156:       PetscCheck(loc >= 0 && loc < *l_size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Location %" PetscInt_FMT " not in [0, %" PetscInt_FMT ") local vector", loc, *l_size);
    157:     }
    158:     PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
    159:   }
    160:   PetscCheck(cell_offset == restr_size, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, offsets array of shape (%" PetscInt_FMT ", %" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", *num_cells, *cell_size, cell_offset);
    161:   if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
    162:   PetscCall(ISDestroy(&iter_is));

    164:   *offsets = restr_indices;
    165:   PetscCall(PetscLogEventEnd(DMPLEX_GetLocalOffsets, dm, 0, 0, 0));
    166:   PetscFunctionReturn(PETSC_SUCCESS);
    167: }

    169: /*@C
    170:   DMPlexGetLocalOffsetsSupport - Allocate and populate arrays of local offsets for each face support.

    172:   Not collective

    174:   Input Parameters:
    175: + dm           - The `DMPLEX` object
    176: . domain_label - label for `DMPLEX` domain, or NULL for whole domain
    177: - label_value  - Stratum value

    179:   Output Parameters:
    180: + num_faces  - Number of local, non-boundary faces
    181: . num_comp   - Number of components per dof
    182: . l_size     - Size of local vector
    183: . offsetsNeg - Allocated offsets array for cells on the inward normal side of each face
    184: - offsetsPos - Allocated offsets array for cells on the outward normal side of each face

    186:   Level: developer

    188:   Notes:
    189:   Allocate and populate array of shape [num_cells, num_comp] defining offsets for each cell for local vector of the `DMPLEX` field. All offsets are in the range [0, l_size - 1].

    191:   Caller is responsible for freeing the offsets array using `PetscFree()`.

    193: .seealso: [](ch_unstructured), `DMPlexGetLocalOffsets()`, `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetClosureIndices()`, `DMPlexSetClosurePermutationTensor()`, `DMPlexGetCeedRestriction()`
    194: @*/
    195: PetscErrorCode DMPlexGetLocalOffsetsSupport(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt *num_faces, PetscInt *num_comp, PetscInt *l_size, PetscInt **offsetsNeg, PetscInt **offsetsPos)
    196: {
    197:   PetscDS         ds = NULL;
    198:   PetscFV         fv;
    199:   PetscSection    section;
    200:   PetscInt        dim, height = 1, dm_field = 0, ds_field = 0, Nf, NfInt = 0, Nc;
    201:   PetscInt       *restr_indices_neg, *restr_indices_pos;
    202:   const PetscInt *iter_indices;
    203:   IS              iter_is;

    205:   PetscFunctionBeginUser;
    207:   PetscCall(DMGetDimension(dm, &dim));
    208:   PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL));
    209:   PetscCall(DMGetLocalSection(dm, &section));
    210:   PetscCall(PetscSectionGetStorageSize(section, l_size));

    212:   PetscCall(DMGetPoints_Internal(dm, domain_label, label_value, height, &iter_is));
    213:   if (iter_is) {
    214:     PetscCall(ISGetIndices(iter_is, &iter_indices));
    215:     PetscCall(ISGetLocalSize(iter_is, &Nf));
    216:     for (PetscInt p = 0, Ns; p < Nf; ++p) {
    217:       PetscCall(DMPlexGetSupportSize(dm, iter_indices[p], &Ns));
    218:       if (Ns == 2) ++NfInt;
    219:     }
    220:     *num_faces = NfInt;
    221:   } else {
    222:     *num_faces   = 0;
    223:     iter_indices = NULL;
    224:   }

    226:   PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fv));
    227:   PetscCall(PetscFVGetNumComponents(fv, &Nc));
    228:   PetscCall(PetscMalloc1(NfInt, &restr_indices_neg));
    229:   PetscCall(PetscMalloc1(NfInt, &restr_indices_pos));
    230:   PetscInt face_offset_neg = 0, face_offset_pos = 0;

    232:   for (PetscInt p = 0; p < Nf; ++p) {
    233:     const PetscInt  face = iter_indices[p];
    234:     PetscInt        num_indices, *indices;
    235:     PetscInt        field_offsets[17]; // max number of fields plus 1
    236:     const PetscInt *supp;
    237:     PetscInt        Ns, loc;

    239:     PetscCall(DMPlexGetSupport(dm, face, &supp));
    240:     PetscCall(DMPlexGetSupportSize(dm, face, &Ns));
    241:     // Ignore boundary faces
    242:     //   TODO check for face on parallel boundary
    243:     if (Ns == 2) {
    244:       // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode.
    245:       PetscCall(DMPlexGetClosureIndices(dm, section, section, supp[0], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
    246:       PetscCheck(num_indices == Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of closure indices %" PetscInt_FMT " != %" PetscInt_FMT " number of FV components", num_indices, Nc);
    247:       loc                                  = indices[field_offsets[dm_field]];
    248:       loc                                  = loc < 0 ? -(loc + 1) : loc;
    249:       restr_indices_neg[face_offset_neg++] = loc;
    250:       PetscCheck(loc >= 0 && loc + Nc <= *l_size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Location %" PetscInt_FMT " + Nc not in [0, %" PetscInt_FMT ") local vector", loc, *l_size);
    251:       PetscCall(DMPlexRestoreClosureIndices(dm, section, section, supp[0], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
    252:       PetscCall(DMPlexGetClosureIndices(dm, section, section, supp[1], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
    253:       PetscCheck(num_indices == Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of closure indices %" PetscInt_FMT " != %" PetscInt_FMT " number of FV components", num_indices, Nc);
    254:       loc                                  = indices[field_offsets[dm_field]];
    255:       loc                                  = loc < 0 ? -(loc + 1) : loc;
    256:       restr_indices_pos[face_offset_pos++] = loc;
    257:       PetscCheck(loc >= 0 && loc + Nc <= *l_size, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Location %" PetscInt_FMT " + Nc not in [0, %" PetscInt_FMT ") local vector", loc, *l_size);
    258:       PetscCall(DMPlexRestoreClosureIndices(dm, section, section, supp[1], PETSC_TRUE, &num_indices, &indices, field_offsets, NULL));
    259:     }
    260:   }
    261:   PetscCheck(face_offset_neg == NfInt, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, neg offsets array of shape (%" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", NfInt, face_offset_neg);
    262:   PetscCheck(face_offset_pos == NfInt, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, pos offsets array of shape (%" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", NfInt, face_offset_pos);
    263:   if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
    264:   PetscCall(ISDestroy(&iter_is));

    266:   *num_comp   = Nc;
    267:   *offsetsNeg = restr_indices_neg;
    268:   *offsetsPos = restr_indices_pos;
    269:   PetscFunctionReturn(PETSC_SUCCESS);
    270: }

    272: #if defined(PETSC_HAVE_LIBCEED)
    273: #include <petscdmplexceed.h>

    275: // Consumes the input petsc_indices and provides the output ceed_indices; no-copy when the sizes match.
    276: static PetscErrorCode PetscIntArrayIntoCeedInt_Private(PetscInt length, PetscInt max_bound, const char *max_bound_name, PetscInt **petsc_indices, CeedInt **ceed_indices)
    277: {
    278:   PetscFunctionBegin;
    279:   if (length) PetscAssertPointer(petsc_indices, 3);
    280:   PetscAssertPointer(ceed_indices, 4);
    281:   #if defined(PETSC_USE_64BIT_INDICES)
    282:   PetscCheck(max_bound <= PETSC_INT32_MAX, PETSC_COMM_SELF, PETSC_ERR_SUP, "%s %" PetscInt_FMT " does not fit in int32_t", max_bound_name, max_bound);
    283:   {
    284:     CeedInt *ceed_ind;
    285:     PetscCall(PetscMalloc1(length, &ceed_ind));
    286:     for (PetscInt i = 0; i < length; i++) ceed_ind[i] = (*petsc_indices)[i];
    287:     *ceed_indices = ceed_ind;
    288:     PetscCall(PetscFree(*petsc_indices));
    289:   }
    290:   #else
    291:   *ceed_indices  = *petsc_indices;
    292:   *petsc_indices = NULL;
    293:   #endif
    294:   PetscFunctionReturn(PETSC_SUCCESS);
    295: }

    297: /*@C
    298:   DMPlexGetCeedRestriction - Define the libCEED map from the local vector (Lvector) to the cells (Evector)

    300:   Input Parameters:
    301: +  dm - The `DMPLEX` object
    302: .  domain_label - label for `DMPLEX` domain, or NULL for the whole domain
    303: .  label_value - Stratum value
    304: .  height - Height of target cells in `DMPLEX` topology
    305: -  dm_field - Index of `DMPLEX` field

    307:   Output Parameter:
    308: .  ERestrict - libCEED restriction from local vector to the cells

    310:   Level: developer

    312: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexGetLocalOffsets()`
    313: @*/
    314: PetscErrorCode DMPlexGetCeedRestriction(DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, CeedElemRestriction *ERestrict)
    315: {
    316:   PetscFunctionBeginUser;
    318:   PetscAssertPointer(ERestrict, 6);
    319:   if (!dm->ceedERestrict) {
    320:     PetscInt            num_cells, cell_size, num_comp, lvec_size, *restr_indices;
    321:     CeedElemRestriction elem_restr;
    322:     Ceed                ceed;

    324:     PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, dm_field, &num_cells, &cell_size, &num_comp, &lvec_size, &restr_indices));
    325:     PetscCall(DMGetCeed(dm, &ceed));
    326:     {
    327:       CeedInt *ceed_indices;
    328:       PetscCall(PetscIntArrayIntoCeedInt_Private(num_cells * cell_size, lvec_size, "lvec_size", &restr_indices, &ceed_indices));
    329:       PetscCallCEED(CeedElemRestrictionCreate(ceed, num_cells, cell_size, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, ceed_indices, &elem_restr));
    330:       PetscCall(PetscFree(ceed_indices));
    331:     }
    332:     dm->ceedERestrict = elem_restr;
    333:   }
    334:   *ERestrict = dm->ceedERestrict;
    335:   PetscFunctionReturn(PETSC_SUCCESS);
    336: }

    338: PetscErrorCode DMPlexCreateCeedRestrictionFVM(DM dm, CeedElemRestriction *erL, CeedElemRestriction *erR)
    339: {
    340:   Ceed      ceed;
    341:   PetscInt *offL, *offR;
    342:   PetscInt  num_faces, num_comp, lvec_size;

    344:   PetscFunctionBeginUser;
    346:   PetscAssertPointer(erL, 2);
    347:   PetscAssertPointer(erR, 3);
    348:   PetscCall(DMGetCeed(dm, &ceed));
    349:   PetscCall(DMPlexGetLocalOffsetsSupport(dm, NULL, 0, &num_faces, &num_comp, &lvec_size, &offL, &offR));
    350:   {
    351:     CeedInt *ceed_off;
    352:     PetscCall(PetscIntArrayIntoCeedInt_Private(num_faces * 1, lvec_size, "lvec_size", &offL, &ceed_off));
    353:     PetscCallCEED(CeedElemRestrictionCreate(ceed, num_faces, 1, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, ceed_off, erL));
    354:     PetscCall(PetscFree(ceed_off));
    355:     PetscCall(PetscIntArrayIntoCeedInt_Private(num_faces * 1, lvec_size, "lvec_size", &offR, &ceed_off));
    356:     PetscCallCEED(CeedElemRestrictionCreate(ceed, num_faces, 1, num_comp, 1, lvec_size, CEED_MEM_HOST, CEED_COPY_VALUES, ceed_off, erR));
    357:     PetscCall(PetscFree(ceed_off));
    358:   }
    359:   PetscFunctionReturn(PETSC_SUCCESS);
    360: }

    362: // TODO DMPlexComputeGeometryFVM() also computes centroids and minimum radius
    363: // TODO DMPlexComputeGeometryFVM() flips normal to match support orientation
    364: // This function computes area-weights normals
    365: PetscErrorCode DMPlexCeedComputeGeometryFVM(DM dm, CeedVector qd)
    366: {
    367:   DMLabel         domain_label = NULL;
    368:   PetscInt        label_value = 0, height = 1, Nf, cdim;
    369:   const PetscInt *iter_indices;
    370:   IS              iter_is;
    371:   CeedScalar     *qdata;

    373:   PetscFunctionBegin;
    374:   PetscCall(DMGetCoordinateDim(dm, &cdim));
    375:   PetscCall(DMGetPoints_Internal(dm, domain_label, label_value, height, &iter_is));
    376:   if (iter_is) {
    377:     PetscCall(ISGetIndices(iter_is, &iter_indices));
    378:     PetscCall(ISGetLocalSize(iter_is, &Nf));
    379:     for (PetscInt p = 0, Ns; p < Nf; ++p) PetscCall(DMPlexGetSupportSize(dm, iter_indices[p], &Ns));
    380:   } else {
    381:     iter_indices = NULL;
    382:   }

    384:   PetscCallCEED(CeedVectorSetValue(qd, 0.));
    385:   PetscCallCEED(CeedVectorGetArray(qd, CEED_MEM_HOST, &qdata));
    386:   for (PetscInt p = 0, off = 0; p < Nf; ++p) {
    387:     const PetscInt  face = iter_indices[p];
    388:     const PetscInt *supp;
    389:     PetscInt        suppSize;

    391:     PetscCall(DMPlexGetSupport(dm, face, &supp));
    392:     PetscCall(DMPlexGetSupportSize(dm, face, &suppSize));
    393:     // Ignore boundary faces
    394:     //   TODO check for face on parallel boundary
    395:     if (suppSize == 2) {
    396:       DMPolytopeType ct;
    397:       PetscReal      area, fcentroid[3], centroids[2][3];

    399:       PetscCall(DMPlexComputeCellGeometryFVM(dm, face, &area, fcentroid, &qdata[off]));
    400:       for (PetscInt d = 0; d < cdim; ++d) qdata[off + d] *= area;
    401:       off += cdim;
    402:       for (PetscInt s = 0; s < suppSize; ++s) {
    403:         PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
    404:         if (ct == DM_POLYTOPE_FV_GHOST) continue;
    405:         PetscCall(DMPlexComputeCellGeometryFVM(dm, supp[s], &qdata[off + s], centroids[s], NULL));
    406:       }
    407:       // Give FV ghosts the same volume as the opposite cell
    408:       for (PetscInt s = 0; s < suppSize; ++s) {
    409:         PetscCall(DMPlexGetCellType(dm, supp[s], &ct));
    410:         if (ct != DM_POLYTOPE_FV_GHOST) continue;
    411:         qdata[off + s] = qdata[off + (1 - s)];
    412:         for (PetscInt d = 0; d < cdim; ++d) centroids[s][d] = fcentroid[d];
    413:       }
    414:       // Flip normal orientation if necessary to match ordering in support
    415:       {
    416:         CeedScalar *normal = &qdata[off - cdim];
    417:         PetscReal   l[3], r[3], v[3];

    419:         PetscCall(DMLocalizeCoordinateReal_Internal(dm, cdim, fcentroid, centroids[0], l));
    420:         PetscCall(DMLocalizeCoordinateReal_Internal(dm, cdim, fcentroid, centroids[1], r));
    421:         DMPlex_WaxpyD_Internal(cdim, -1, l, r, v);
    422:         if (DMPlex_DotRealD_Internal(cdim, normal, v) < 0) {
    423:           for (PetscInt d = 0; d < cdim; ++d) normal[d] = -normal[d];
    424:         }
    425:         if (DMPlex_DotRealD_Internal(cdim, normal, v) <= 0) {
    426:           PetscCheck(cdim != 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed, normal (%g,%g) v (%g,%g)", face, (double)normal[0], (double)normal[1], (double)v[0], (double)v[1]);
    427:           PetscCheck(cdim != 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", face, (double)normal[0], (double)normal[1], (double)normal[2], (double)v[0], (double)v[1], (double)v[2]);
    428:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", face);
    429:         }
    430:       }
    431:       off += suppSize;
    432:     }
    433:   }
    434:   PetscCallCEED(CeedVectorRestoreArray(qd, &qdata));
    435:   if (iter_is) PetscCall(ISRestoreIndices(iter_is, &iter_indices));
    436:   PetscCall(ISDestroy(&iter_is));
    437:   PetscFunctionReturn(PETSC_SUCCESS);
    438: }

    440: #endif