Actual source code: plextree.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/isimpl.h>
      3: #include <petsc/private/petscfeimpl.h>
      4: #include <petscsf.h>
      5: #include <petscds.h>

      7: /* hierarchy routines */

      9: /*@
     10:   DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes.

     12:   Not Collective

     14:   Input Parameters:
     15: + dm  - The `DMPLEX` object
     16: - ref - The reference tree `DMPLEX` object

     18:   Level: intermediate

     20: .seealso: [](ch_unstructured), `DM`, `DMPLEX`,`DMPlexGetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`
     21: @*/
     22: PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref)
     23: {
     24:   DM_Plex *mesh = (DM_Plex *)dm->data;

     26:   PetscFunctionBegin;
     29:   PetscCall(PetscObjectReference((PetscObject)ref));
     30:   PetscCall(DMDestroy(&mesh->referenceTree));
     31:   mesh->referenceTree = ref;
     32:   PetscFunctionReturn(PETSC_SUCCESS);
     33: }

     35: /*@
     36:   DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes.

     38:   Not Collective

     40:   Input Parameter:
     41: . dm - The `DMPLEX` object

     43:   Output Parameter:
     44: . ref - The reference tree `DMPLEX` object

     46:   Level: intermediate

     48:   Developer Notes:
     49:   The reference tree is shallow copied during `DMClone()`, thus it is may be shared by different `DM`s.
     50:   It is not a topological-only object, since some parts of the library use its local section to compute
     51:   interpolation and injection matrices. This may lead to unexpected failures during those calls.

     53: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`
     54: @*/
     55: PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref)
     56: {
     57:   DM_Plex *mesh = (DM_Plex *)dm->data;

     59:   PetscFunctionBegin;
     61:   PetscAssertPointer(ref, 2);
     62:   *ref = mesh->referenceTree;
     63:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
     67: {
     68:   PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert;

     70:   PetscFunctionBegin;
     71:   if (parentOrientA == parentOrientB) {
     72:     if (childOrientB) *childOrientB = childOrientA;
     73:     if (childB) *childB = childA;
     74:     PetscFunctionReturn(PETSC_SUCCESS);
     75:   }
     76:   for (dim = 0; dim < 3; dim++) {
     77:     PetscCall(DMPlexGetDepthStratum(dm, dim, &dStart, &dEnd));
     78:     if (parent >= dStart && parent <= dEnd) break;
     79:   }
     80:   PetscCheck(dim <= 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform child symmetry for %" PetscInt_FMT "-cells", dim);
     81:   PetscCheck(dim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A vertex has no children");
     82:   if (childA < dStart || childA >= dEnd) {
     83:     /* this is a lower-dimensional child: bootstrap */
     84:     PetscInt        size, i, sA = -1, sB, sOrientB, sConeSize;
     85:     const PetscInt *supp, *coneA, *coneB, *oA, *oB;

     87:     PetscCall(DMPlexGetSupportSize(dm, childA, &size));
     88:     PetscCall(DMPlexGetSupport(dm, childA, &supp));

     90:     /* find a point sA in supp(childA) that has the same parent */
     91:     for (i = 0; i < size; i++) {
     92:       PetscInt sParent;

     94:       sA = supp[i];
     95:       if (sA == parent) continue;
     96:       PetscCall(DMPlexGetTreeParent(dm, sA, &sParent, NULL));
     97:       if (sParent == parent) break;
     98:     }
     99:     PetscCheck(i != size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "could not find support in children");
    100:     /* find out which point sB is in an equivalent position to sA under
    101:      * parentOrientB */
    102:     PetscCall(DMPlexReferenceTreeGetChildSymmetry_Default(dm, parent, parentOrientA, 0, sA, parentOrientB, &sOrientB, &sB));
    103:     PetscCall(DMPlexGetConeSize(dm, sA, &sConeSize));
    104:     PetscCall(DMPlexGetCone(dm, sA, &coneA));
    105:     PetscCall(DMPlexGetCone(dm, sB, &coneB));
    106:     PetscCall(DMPlexGetConeOrientation(dm, sA, &oA));
    107:     PetscCall(DMPlexGetConeOrientation(dm, sB, &oB));
    108:     /* step through the cone of sA in natural order */
    109:     for (i = 0; i < sConeSize; i++) {
    110:       if (coneA[i] == childA) {
    111:         /* if childA is at position i in coneA,
    112:          * then we want the point that is at sOrientB*i in coneB */
    113:         PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize - (sOrientB + 1) - i) % sConeSize);
    114:         if (childB) *childB = coneB[j];
    115:         if (childOrientB) {
    116:           DMPolytopeType ct;
    117:           PetscInt       oBtrue;

    119:           PetscCall(DMPlexGetConeSize(dm, childA, &coneSize));
    120:           /* compose sOrientB and oB[j] */
    121:           PetscCheck(coneSize == 0 || coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected a vertex or an edge");
    122:           ct = coneSize ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;
    123:           /* we may have to flip an edge */
    124:           oBtrue        = (sOrientB >= 0) ? oB[j] : DMPolytopeTypeComposeOrientation(ct, -1, oB[j]);
    125:           oBtrue        = DMPolytopeConvertNewOrientation_Internal(ct, oBtrue);
    126:           ABswap        = DihedralSwap(coneSize, DMPolytopeConvertNewOrientation_Internal(ct, oA[i]), oBtrue);
    127:           *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
    128:         }
    129:         break;
    130:       }
    131:     }
    132:     PetscCheck(i != sConeSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "support cone mismatch");
    133:     PetscFunctionReturn(PETSC_SUCCESS);
    134:   }
    135:   /* get the cone size and symmetry swap */
    136:   PetscCall(DMPlexGetConeSize(dm, parent, &coneSize));
    137:   ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB);
    138:   if (dim == 2) {
    139:     /* orientations refer to cones: we want them to refer to vertices:
    140:      * if it's a rotation, they are the same, but if the order is reversed, a
    141:      * permutation that puts side i first does *not* put vertex i first */
    142:     oAvert     = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1);
    143:     oBvert     = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1);
    144:     ABswapVert = DihedralSwap(coneSize, oAvert, oBvert);
    145:   } else {
    146:     ABswapVert = ABswap;
    147:   }
    148:   if (childB) {
    149:     /* assume that each child corresponds to a vertex, in the same order */
    150:     PetscInt        p, posA = -1, numChildren, i;
    151:     const PetscInt *children;

    153:     /* count which position the child is in */
    154:     PetscCall(DMPlexGetTreeChildren(dm, parent, &numChildren, &children));
    155:     for (i = 0; i < numChildren; i++) {
    156:       p = children[i];
    157:       if (p == childA) {
    158:         posA = i;
    159:         break;
    160:       }
    161:     }
    162:     if (posA >= coneSize) {
    163:       /* this is the triangle in the middle of a uniformly refined triangle: it
    164:        * is invariant */
    165:       PetscCheck(dim == 2 && posA == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected a middle triangle, got something else");
    166:       *childB = childA;
    167:     } else {
    168:       /* figure out position B by applying ABswapVert */
    169:       PetscInt posB;

    171:       posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize - (ABswapVert + 1) - posA) % coneSize);
    172:       if (childB) *childB = children[posB];
    173:     }
    174:   }
    175:   if (childOrientB) *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: /*@
    180:   DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another

    182:   Input Parameters:
    183: + dm            - the reference tree `DMPLEX` object
    184: . parent        - the parent point
    185: . parentOrientA - the reference orientation for describing the parent
    186: . childOrientA  - the reference orientation for describing the child
    187: . childA        - the reference childID for describing the child
    188: - parentOrientB - the new orientation for describing the parent

    190:   Output Parameters:
    191: + childOrientB - if not `NULL`, set to the new orientation for describing the child
    192: - childB       - if not `NULL`, the new childID for describing the child

    194:   Level: developer

    196: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetReferenceTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetTree()`
    197: @*/
    198: PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
    199: {
    200:   DM_Plex *mesh = (DM_Plex *)dm->data;

    202:   PetscFunctionBegin;
    204:   PetscCheck(mesh->getchildsymmetry, PETSC_COMM_SELF, PETSC_ERR_SUP, "DMPlexReferenceTreeGetChildSymmetry not implemented");
    205:   PetscCall(mesh->getchildsymmetry(dm, parent, parentOrientA, childOrientA, childA, parentOrientB, childOrientB, childB));
    206:   PetscFunctionReturn(PETSC_SUCCESS);
    207: }

    209: static PetscErrorCode DMPlexSetTree_Internal(DM, PetscSection, PetscInt *, PetscInt *, PetscBool, PetscBool);

    211: PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
    212: {
    213:   PetscFunctionBegin;
    214:   PetscCall(DMPlexSetTree_Internal(dm, parentSection, parents, childIDs, PETSC_TRUE, PETSC_FALSE));
    215:   PetscFunctionReturn(PETSC_SUCCESS);
    216: }

    218: PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref)
    219: {
    220:   MPI_Comm     comm;
    221:   PetscInt     dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs;
    222:   PetscInt    *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts;
    223:   DMLabel      identity, identityRef;
    224:   PetscSection unionSection, unionConeSection, parentSection;
    225:   PetscScalar *unionCoords;
    226:   IS           perm;

    228:   PetscFunctionBegin;
    229:   comm = PetscObjectComm((PetscObject)K);
    230:   PetscCall(DMGetDimension(K, &dim));
    231:   PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
    232:   PetscCall(DMGetLabel(K, labelName, &identity));
    233:   PetscCall(DMGetLabel(Kref, labelName, &identityRef));
    234:   PetscCall(DMPlexGetChart(Kref, &pRefStart, &pRefEnd));
    235:   PetscCall(PetscSectionCreate(comm, &unionSection));
    236:   PetscCall(PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart)));
    237:   /* count points that will go in the union */
    238:   for (p = pStart; p < pEnd; p++) PetscCall(PetscSectionSetDof(unionSection, p - pStart, 1));
    239:   for (p = pRefStart; p < pRefEnd; p++) {
    240:     PetscInt q, qSize;
    241:     PetscCall(DMLabelGetValue(identityRef, p, &q));
    242:     PetscCall(DMLabelGetStratumSize(identityRef, q, &qSize));
    243:     if (qSize > 1) PetscCall(PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1));
    244:   }
    245:   PetscCall(PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart, &permvals));
    246:   offset = 0;
    247:   /* stratify points in the union by topological dimension */
    248:   for (d = 0; d <= dim; d++) {
    249:     PetscInt cStart, cEnd, c;

    251:     PetscCall(DMPlexGetHeightStratum(K, d, &cStart, &cEnd));
    252:     for (c = cStart; c < cEnd; c++) permvals[offset++] = c;

    254:     PetscCall(DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd));
    255:     for (c = cStart; c < cEnd; c++) permvals[offset++] = c + (pEnd - pStart);
    256:   }
    257:   PetscCall(ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm));
    258:   PetscCall(PetscSectionSetPermutation(unionSection, perm));
    259:   PetscCall(PetscSectionSetUp(unionSection));
    260:   PetscCall(PetscSectionGetStorageSize(unionSection, &numUnionPoints));
    261:   PetscCall(PetscMalloc2(numUnionPoints, &coneSizes, dim + 1, &numDimPoints));
    262:   /* count dimension points */
    263:   for (d = 0; d <= dim; d++) {
    264:     PetscInt cStart, cOff, cOff2;
    265:     PetscCall(DMPlexGetHeightStratum(K, d, &cStart, NULL));
    266:     PetscCall(PetscSectionGetOffset(unionSection, cStart - pStart, &cOff));
    267:     if (d < dim) {
    268:       PetscCall(DMPlexGetHeightStratum(K, d + 1, &cStart, NULL));
    269:       PetscCall(PetscSectionGetOffset(unionSection, cStart - pStart, &cOff2));
    270:     } else {
    271:       cOff2 = numUnionPoints;
    272:     }
    273:     numDimPoints[dim - d] = cOff2 - cOff;
    274:   }
    275:   PetscCall(PetscSectionCreate(comm, &unionConeSection));
    276:   PetscCall(PetscSectionSetChart(unionConeSection, 0, numUnionPoints));
    277:   /* count the cones in the union */
    278:   for (p = pStart; p < pEnd; p++) {
    279:     PetscInt dof, uOff;

    281:     PetscCall(DMPlexGetConeSize(K, p, &dof));
    282:     PetscCall(PetscSectionGetOffset(unionSection, p - pStart, &uOff));
    283:     PetscCall(PetscSectionSetDof(unionConeSection, uOff, dof));
    284:     coneSizes[uOff] = dof;
    285:   }
    286:   for (p = pRefStart; p < pRefEnd; p++) {
    287:     PetscInt dof, uDof, uOff;

    289:     PetscCall(DMPlexGetConeSize(Kref, p, &dof));
    290:     PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
    291:     PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
    292:     if (uDof) {
    293:       PetscCall(PetscSectionSetDof(unionConeSection, uOff, dof));
    294:       coneSizes[uOff] = dof;
    295:     }
    296:   }
    297:   PetscCall(PetscSectionSetUp(unionConeSection));
    298:   PetscCall(PetscSectionGetStorageSize(unionConeSection, &numCones));
    299:   PetscCall(PetscMalloc2(numCones, &unionCones, numCones, &unionOrientations));
    300:   /* write the cones in the union */
    301:   for (p = pStart; p < pEnd; p++) {
    302:     PetscInt        dof, uOff, c, cOff;
    303:     const PetscInt *cone, *orientation;

    305:     PetscCall(DMPlexGetConeSize(K, p, &dof));
    306:     PetscCall(DMPlexGetCone(K, p, &cone));
    307:     PetscCall(DMPlexGetConeOrientation(K, p, &orientation));
    308:     PetscCall(PetscSectionGetOffset(unionSection, p - pStart, &uOff));
    309:     PetscCall(PetscSectionGetOffset(unionConeSection, uOff, &cOff));
    310:     for (c = 0; c < dof; c++) {
    311:       PetscInt e, eOff;
    312:       e = cone[c];
    313:       PetscCall(PetscSectionGetOffset(unionSection, e - pStart, &eOff));
    314:       unionCones[cOff + c]        = eOff;
    315:       unionOrientations[cOff + c] = orientation[c];
    316:     }
    317:   }
    318:   for (p = pRefStart; p < pRefEnd; p++) {
    319:     PetscInt        dof, uDof, uOff, c, cOff;
    320:     const PetscInt *cone, *orientation;

    322:     PetscCall(DMPlexGetConeSize(Kref, p, &dof));
    323:     PetscCall(DMPlexGetCone(Kref, p, &cone));
    324:     PetscCall(DMPlexGetConeOrientation(Kref, p, &orientation));
    325:     PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
    326:     PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
    327:     if (uDof) {
    328:       PetscCall(PetscSectionGetOffset(unionConeSection, uOff, &cOff));
    329:       for (c = 0; c < dof; c++) {
    330:         PetscInt e, eOff, eDof;

    332:         e = cone[c];
    333:         PetscCall(PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart), &eDof));
    334:         if (eDof) {
    335:           PetscCall(PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff));
    336:         } else {
    337:           PetscCall(DMLabelGetValue(identityRef, e, &e));
    338:           PetscCall(PetscSectionGetOffset(unionSection, e - pStart, &eOff));
    339:         }
    340:         unionCones[cOff + c]        = eOff;
    341:         unionOrientations[cOff + c] = orientation[c];
    342:       }
    343:     }
    344:   }
    345:   /* get the coordinates */
    346:   {
    347:     PetscInt     vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff;
    348:     PetscSection KcoordsSec, KrefCoordsSec;
    349:     Vec          KcoordsVec, KrefCoordsVec;
    350:     PetscScalar *Kcoords;

    352:     PetscCall(DMGetCoordinateSection(K, &KcoordsSec));
    353:     PetscCall(DMGetCoordinatesLocal(K, &KcoordsVec));
    354:     PetscCall(DMGetCoordinateSection(Kref, &KrefCoordsSec));
    355:     PetscCall(DMGetCoordinatesLocal(Kref, &KrefCoordsVec));

    357:     numVerts = numDimPoints[0];
    358:     PetscCall(PetscMalloc1(numVerts * dim, &unionCoords));
    359:     PetscCall(DMPlexGetDepthStratum(K, 0, &vStart, &vEnd));

    361:     offset = 0;
    362:     for (v = vStart; v < vEnd; v++) {
    363:       PetscCall(PetscSectionGetOffset(unionSection, v - pStart, &vOff));
    364:       PetscCall(VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords));
    365:       for (d = 0; d < dim; d++) unionCoords[offset * dim + d] = Kcoords[d];
    366:       offset++;
    367:     }
    368:     PetscCall(DMPlexGetDepthStratum(Kref, 0, &vRefStart, &vRefEnd));
    369:     for (v = vRefStart; v < vRefEnd; v++) {
    370:       PetscCall(PetscSectionGetDof(unionSection, v - pRefStart + (pEnd - pStart), &vDof));
    371:       PetscCall(PetscSectionGetOffset(unionSection, v - pRefStart + (pEnd - pStart), &vOff));
    372:       PetscCall(VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords));
    373:       if (vDof) {
    374:         for (d = 0; d < dim; d++) unionCoords[offset * dim + d] = Kcoords[d];
    375:         offset++;
    376:       }
    377:     }
    378:   }
    379:   PetscCall(DMCreate(comm, ref));
    380:   PetscCall(DMSetType(*ref, DMPLEX));
    381:   PetscCall(DMSetDimension(*ref, dim));
    382:   PetscCall(DMPlexCreateFromDAG(*ref, dim, numDimPoints, coneSizes, unionCones, unionOrientations, unionCoords));
    383:   /* set the tree */
    384:   PetscCall(PetscSectionCreate(comm, &parentSection));
    385:   PetscCall(PetscSectionSetChart(parentSection, 0, numUnionPoints));
    386:   for (p = pRefStart; p < pRefEnd; p++) {
    387:     PetscInt uDof, uOff;

    389:     PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
    390:     PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
    391:     if (uDof) PetscCall(PetscSectionSetDof(parentSection, uOff, 1));
    392:   }
    393:   PetscCall(PetscSectionSetUp(parentSection));
    394:   PetscCall(PetscSectionGetStorageSize(parentSection, &parentSize));
    395:   PetscCall(PetscMalloc2(parentSize, &parents, parentSize, &childIDs));
    396:   for (p = pRefStart; p < pRefEnd; p++) {
    397:     PetscInt uDof, uOff;

    399:     PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
    400:     PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
    401:     if (uDof) {
    402:       PetscInt pOff, parent, parentU;
    403:       PetscCall(PetscSectionGetOffset(parentSection, uOff, &pOff));
    404:       PetscCall(DMLabelGetValue(identityRef, p, &parent));
    405:       PetscCall(PetscSectionGetOffset(unionSection, parent - pStart, &parentU));
    406:       parents[pOff]  = parentU;
    407:       childIDs[pOff] = uOff;
    408:     }
    409:   }
    410:   PetscCall(DMPlexCreateReferenceTree_SetTree(*ref, parentSection, parents, childIDs));
    411:   PetscCall(PetscSectionDestroy(&parentSection));
    412:   PetscCall(PetscFree2(parents, childIDs));

    414:   /* clean up */
    415:   PetscCall(PetscSectionDestroy(&unionSection));
    416:   PetscCall(PetscSectionDestroy(&unionConeSection));
    417:   PetscCall(ISDestroy(&perm));
    418:   PetscCall(PetscFree(unionCoords));
    419:   PetscCall(PetscFree2(unionCones, unionOrientations));
    420:   PetscCall(PetscFree2(coneSizes, numDimPoints));
    421:   PetscFunctionReturn(PETSC_SUCCESS);
    422: }

    424: /*@
    425:   DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement.

    427:   Collective

    429:   Input Parameters:
    430: + comm    - the MPI communicator
    431: . dim     - the spatial dimension
    432: - simplex - Flag for simplex, otherwise use a tensor-product cell

    434:   Output Parameter:
    435: . ref - the reference tree `DMPLEX` object

    437:   Level: intermediate

    439: .seealso: `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`
    440: @*/
    441: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
    442: {
    443:   DM_Plex *mesh;
    444:   DM       K, Kref;
    445:   PetscInt p, pStart, pEnd;
    446:   DMLabel  identity;

    448:   PetscFunctionBegin;
    449: #if 1
    450:   comm = PETSC_COMM_SELF;
    451: #endif
    452:   /* create a reference element */
    453:   PetscCall(DMPlexCreateReferenceCell(comm, DMPolytopeTypeSimpleShape(dim, simplex), &K));
    454:   PetscCall(DMCreateLabel(K, "identity"));
    455:   PetscCall(DMGetLabel(K, "identity", &identity));
    456:   PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
    457:   for (p = pStart; p < pEnd; p++) PetscCall(DMLabelSetValue(identity, p, p));
    458:   /* refine it */
    459:   PetscCall(DMRefine(K, comm, &Kref));

    461:   /* the reference tree is the union of these two, without duplicating
    462:    * points that appear in both */
    463:   PetscCall(DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref));
    464:   mesh                   = (DM_Plex *)(*ref)->data;
    465:   mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
    466:   PetscCall(DMDestroy(&K));
    467:   PetscCall(DMDestroy(&Kref));
    468:   PetscFunctionReturn(PETSC_SUCCESS);
    469: }

    471: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
    472: {
    473:   DM_Plex     *mesh = (DM_Plex *)dm->data;
    474:   PetscSection childSec, pSec;
    475:   PetscInt     p, pSize, cSize, parMax = PETSC_INT_MIN, parMin = PETSC_INT_MAX;
    476:   PetscInt    *offsets, *children, pStart, pEnd;

    478:   PetscFunctionBegin;
    480:   PetscCall(PetscSectionDestroy(&mesh->childSection));
    481:   PetscCall(PetscFree(mesh->children));
    482:   pSec = mesh->parentSection;
    483:   if (!pSec) PetscFunctionReturn(PETSC_SUCCESS);
    484:   PetscCall(PetscSectionGetStorageSize(pSec, &pSize));
    485:   for (p = 0; p < pSize; p++) {
    486:     PetscInt par = mesh->parents[p];

    488:     parMax = PetscMax(parMax, par + 1);
    489:     parMin = PetscMin(parMin, par);
    490:   }
    491:   if (parMin > parMax) {
    492:     parMin = -1;
    493:     parMax = -1;
    494:   }
    495:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)pSec), &childSec));
    496:   PetscCall(PetscSectionSetChart(childSec, parMin, parMax));
    497:   for (p = 0; p < pSize; p++) {
    498:     PetscInt par = mesh->parents[p];

    500:     PetscCall(PetscSectionAddDof(childSec, par, 1));
    501:   }
    502:   PetscCall(PetscSectionSetUp(childSec));
    503:   PetscCall(PetscSectionGetStorageSize(childSec, &cSize));
    504:   PetscCall(PetscMalloc1(cSize, &children));
    505:   PetscCall(PetscCalloc1(parMax - parMin, &offsets));
    506:   PetscCall(PetscSectionGetChart(pSec, &pStart, &pEnd));
    507:   for (p = pStart; p < pEnd; p++) {
    508:     PetscInt dof, off, i;

    510:     PetscCall(PetscSectionGetDof(pSec, p, &dof));
    511:     PetscCall(PetscSectionGetOffset(pSec, p, &off));
    512:     for (i = 0; i < dof; i++) {
    513:       PetscInt par = mesh->parents[off + i], cOff;

    515:       PetscCall(PetscSectionGetOffset(childSec, par, &cOff));
    516:       children[cOff + offsets[par - parMin]++] = p;
    517:     }
    518:   }
    519:   mesh->childSection = childSec;
    520:   mesh->children     = children;
    521:   PetscCall(PetscFree(offsets));
    522:   PetscFunctionReturn(PETSC_SUCCESS);
    523: }

    525: static PetscErrorCode AnchorsFlatten(PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
    526: {
    527:   PetscInt        pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
    528:   const PetscInt *vals;
    529:   PetscSection    secNew;
    530:   PetscBool       anyNew, globalAnyNew;
    531:   PetscBool       compress;

    533:   PetscFunctionBegin;
    534:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
    535:   PetscCall(ISGetLocalSize(is, &size));
    536:   PetscCall(ISGetIndices(is, &vals));
    537:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &secNew));
    538:   PetscCall(PetscSectionSetChart(secNew, pStart, pEnd));
    539:   for (i = 0; i < size; i++) {
    540:     PetscInt dof;

    542:     p = vals[i];
    543:     if (p < pStart || p >= pEnd) continue;
    544:     PetscCall(PetscSectionGetDof(section, p, &dof));
    545:     if (dof) break;
    546:   }
    547:   if (i == size) {
    548:     PetscCall(PetscSectionSetUp(secNew));
    549:     anyNew   = PETSC_FALSE;
    550:     compress = PETSC_FALSE;
    551:     sizeNew  = 0;
    552:   } else {
    553:     anyNew = PETSC_TRUE;
    554:     for (p = pStart; p < pEnd; p++) {
    555:       PetscInt dof, off;

    557:       PetscCall(PetscSectionGetDof(section, p, &dof));
    558:       PetscCall(PetscSectionGetOffset(section, p, &off));
    559:       for (i = 0; i < dof; i++) {
    560:         PetscInt q = vals[off + i], qDof = 0;

    562:         if (q >= pStart && q < pEnd) PetscCall(PetscSectionGetDof(section, q, &qDof));
    563:         if (qDof) PetscCall(PetscSectionAddDof(secNew, p, qDof));
    564:         else PetscCall(PetscSectionAddDof(secNew, p, 1));
    565:       }
    566:     }
    567:     PetscCall(PetscSectionSetUp(secNew));
    568:     PetscCall(PetscSectionGetStorageSize(secNew, &sizeNew));
    569:     PetscCall(PetscMalloc1(sizeNew, &valsNew));
    570:     compress = PETSC_FALSE;
    571:     for (p = pStart; p < pEnd; p++) {
    572:       PetscInt dof, off, count, offNew, dofNew;

    574:       PetscCall(PetscSectionGetDof(section, p, &dof));
    575:       PetscCall(PetscSectionGetOffset(section, p, &off));
    576:       PetscCall(PetscSectionGetDof(secNew, p, &dofNew));
    577:       PetscCall(PetscSectionGetOffset(secNew, p, &offNew));
    578:       count = 0;
    579:       for (i = 0; i < dof; i++) {
    580:         PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;

    582:         if (q >= pStart && q < pEnd) {
    583:           PetscCall(PetscSectionGetDof(section, q, &qDof));
    584:           PetscCall(PetscSectionGetOffset(section, q, &qOff));
    585:         }
    586:         if (qDof) {
    587:           PetscInt oldCount = count;

    589:           for (j = 0; j < qDof; j++) {
    590:             PetscInt k, r = vals[qOff + j];

    592:             for (k = 0; k < oldCount; k++) {
    593:               if (valsNew[offNew + k] == r) break;
    594:             }
    595:             if (k == oldCount) valsNew[offNew + count++] = r;
    596:           }
    597:         } else {
    598:           PetscInt k, oldCount = count;

    600:           for (k = 0; k < oldCount; k++) {
    601:             if (valsNew[offNew + k] == q) break;
    602:           }
    603:           if (k == oldCount) valsNew[offNew + count++] = q;
    604:         }
    605:       }
    606:       if (count < dofNew) {
    607:         PetscCall(PetscSectionSetDof(secNew, p, count));
    608:         compress = PETSC_TRUE;
    609:       }
    610:     }
    611:   }
    612:   PetscCall(ISRestoreIndices(is, &vals));
    613:   PetscCallMPI(MPIU_Allreduce(&anyNew, &globalAnyNew, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)secNew)));
    614:   if (!globalAnyNew) {
    615:     PetscCall(PetscSectionDestroy(&secNew));
    616:     *sectionNew = NULL;
    617:     *isNew      = NULL;
    618:   } else {
    619:     if (compress) {
    620:       PetscSection secComp;
    621:       PetscInt    *valsComp = NULL;

    623:       PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &secComp));
    624:       PetscCall(PetscSectionSetChart(secComp, pStart, pEnd));
    625:       for (p = pStart; p < pEnd; p++) {
    626:         PetscInt dof;

    628:         PetscCall(PetscSectionGetDof(secNew, p, &dof));
    629:         PetscCall(PetscSectionSetDof(secComp, p, dof));
    630:       }
    631:       PetscCall(PetscSectionSetUp(secComp));
    632:       PetscCall(PetscSectionGetStorageSize(secComp, &sizeNew));
    633:       PetscCall(PetscMalloc1(sizeNew, &valsComp));
    634:       for (p = pStart; p < pEnd; p++) {
    635:         PetscInt dof, off, offNew, j;

    637:         PetscCall(PetscSectionGetDof(secNew, p, &dof));
    638:         PetscCall(PetscSectionGetOffset(secNew, p, &off));
    639:         PetscCall(PetscSectionGetOffset(secComp, p, &offNew));
    640:         for (j = 0; j < dof; j++) valsComp[offNew + j] = valsNew[off + j];
    641:       }
    642:       PetscCall(PetscSectionDestroy(&secNew));
    643:       secNew = secComp;
    644:       PetscCall(PetscFree(valsNew));
    645:       valsNew = valsComp;
    646:     }
    647:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), sizeNew, valsNew, PETSC_OWN_POINTER, isNew));
    648:   }
    649:   PetscFunctionReturn(PETSC_SUCCESS);
    650: }

    652: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
    653: {
    654:   PetscInt     p, pStart, pEnd, *anchors, size;
    655:   PetscInt     aMin = PETSC_INT_MAX, aMax = PETSC_INT_MIN;
    656:   PetscSection aSec;
    657:   DMLabel      canonLabel;
    658:   IS           aIS;

    660:   PetscFunctionBegin;
    662:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    663:   PetscCall(DMGetLabel(dm, "canonical", &canonLabel));
    664:   for (p = pStart; p < pEnd; p++) {
    665:     PetscInt parent;

    667:     if (canonLabel) {
    668:       PetscInt canon;

    670:       PetscCall(DMLabelGetValue(canonLabel, p, &canon));
    671:       if (p != canon) continue;
    672:     }
    673:     PetscCall(DMPlexGetTreeParent(dm, p, &parent, NULL));
    674:     if (parent != p) {
    675:       aMin = PetscMin(aMin, p);
    676:       aMax = PetscMax(aMax, p + 1);
    677:     }
    678:   }
    679:   if (aMin > aMax) {
    680:     aMin = -1;
    681:     aMax = -1;
    682:   }
    683:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &aSec));
    684:   PetscCall(PetscSectionSetChart(aSec, aMin, aMax));
    685:   for (p = aMin; p < aMax; p++) {
    686:     PetscInt parent, ancestor = p;

    688:     if (canonLabel) {
    689:       PetscInt canon;

    691:       PetscCall(DMLabelGetValue(canonLabel, p, &canon));
    692:       if (p != canon) continue;
    693:     }
    694:     PetscCall(DMPlexGetTreeParent(dm, p, &parent, NULL));
    695:     while (parent != ancestor) {
    696:       ancestor = parent;
    697:       PetscCall(DMPlexGetTreeParent(dm, ancestor, &parent, NULL));
    698:     }
    699:     if (ancestor != p) {
    700:       PetscInt closureSize, *closure = NULL;

    702:       PetscCall(DMPlexGetTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
    703:       PetscCall(PetscSectionSetDof(aSec, p, closureSize));
    704:       PetscCall(DMPlexRestoreTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
    705:     }
    706:   }
    707:   PetscCall(PetscSectionSetUp(aSec));
    708:   PetscCall(PetscSectionGetStorageSize(aSec, &size));
    709:   PetscCall(PetscMalloc1(size, &anchors));
    710:   for (p = aMin; p < aMax; p++) {
    711:     PetscInt parent, ancestor = p;

    713:     if (canonLabel) {
    714:       PetscInt canon;

    716:       PetscCall(DMLabelGetValue(canonLabel, p, &canon));
    717:       if (p != canon) continue;
    718:     }
    719:     PetscCall(DMPlexGetTreeParent(dm, p, &parent, NULL));
    720:     while (parent != ancestor) {
    721:       ancestor = parent;
    722:       PetscCall(DMPlexGetTreeParent(dm, ancestor, &parent, NULL));
    723:     }
    724:     if (ancestor != p) {
    725:       PetscInt j, closureSize, *closure = NULL, aOff;

    727:       PetscCall(PetscSectionGetOffset(aSec, p, &aOff));

    729:       PetscCall(DMPlexGetTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
    730:       for (j = 0; j < closureSize; j++) anchors[aOff + j] = closure[2 * j];
    731:       PetscCall(DMPlexRestoreTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
    732:     }
    733:   }
    734:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, anchors, PETSC_OWN_POINTER, &aIS));
    735:   {
    736:     PetscSection aSecNew = aSec;
    737:     IS           aISNew  = aIS;

    739:     PetscCall(PetscObjectReference((PetscObject)aSec));
    740:     PetscCall(PetscObjectReference((PetscObject)aIS));
    741:     while (aSecNew) {
    742:       PetscCall(PetscSectionDestroy(&aSec));
    743:       PetscCall(ISDestroy(&aIS));
    744:       aSec    = aSecNew;
    745:       aIS     = aISNew;
    746:       aSecNew = NULL;
    747:       aISNew  = NULL;
    748:       PetscCall(AnchorsFlatten(aSec, aIS, &aSecNew, &aISNew));
    749:     }
    750:   }
    751:   PetscCall(DMPlexSetAnchors(dm, aSec, aIS));
    752:   PetscCall(PetscSectionDestroy(&aSec));
    753:   PetscCall(ISDestroy(&aIS));
    754:   PetscFunctionReturn(PETSC_SUCCESS);
    755: }

    757: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm, PetscInt p, PetscInt *dof, PetscInt *numTrueSupp)
    758: {
    759:   PetscFunctionBegin;
    760:   if (numTrueSupp[p] == -1) {
    761:     PetscInt        i, alldof;
    762:     const PetscInt *supp;
    763:     PetscInt        count = 0;

    765:     PetscCall(DMPlexGetSupportSize(dm, p, &alldof));
    766:     PetscCall(DMPlexGetSupport(dm, p, &supp));
    767:     for (i = 0; i < alldof; i++) {
    768:       PetscInt        q = supp[i], numCones, j;
    769:       const PetscInt *cone;

    771:       PetscCall(DMPlexGetConeSize(dm, q, &numCones));
    772:       PetscCall(DMPlexGetCone(dm, q, &cone));
    773:       for (j = 0; j < numCones; j++) {
    774:         if (cone[j] == p) break;
    775:       }
    776:       if (j < numCones) count++;
    777:     }
    778:     numTrueSupp[p] = count;
    779:   }
    780:   *dof = numTrueSupp[p];
    781:   PetscFunctionReturn(PETSC_SUCCESS);
    782: }

    784: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
    785: {
    786:   DM_Plex     *mesh = (DM_Plex *)dm->data;
    787:   PetscSection newSupportSection;
    788:   PetscInt     newSize, *newSupports, pStart, pEnd, p, d, depth;
    789:   PetscInt    *numTrueSupp;
    790:   PetscInt    *offsets;

    792:   PetscFunctionBegin;
    794:   /* symmetrize the hierarchy */
    795:   PetscCall(DMPlexGetDepth(dm, &depth));
    796:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)mesh->supportSection), &newSupportSection));
    797:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    798:   PetscCall(PetscSectionSetChart(newSupportSection, pStart, pEnd));
    799:   PetscCall(PetscCalloc1(pEnd, &offsets));
    800:   PetscCall(PetscMalloc1(pEnd, &numTrueSupp));
    801:   for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
    802:   /* if a point is in the (true) support of q, it should be in the support of
    803:    * parent(q) */
    804:   for (d = 0; d <= depth; d++) {
    805:     PetscCall(DMPlexGetHeightStratum(dm, d, &pStart, &pEnd));
    806:     for (p = pStart; p < pEnd; ++p) {
    807:       PetscInt dof, q, qdof, parent;

    809:       PetscCall(DMPlexGetTrueSupportSize(dm, p, &dof, numTrueSupp));
    810:       PetscCall(PetscSectionAddDof(newSupportSection, p, dof));
    811:       q = p;
    812:       PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
    813:       while (parent != q && parent >= pStart && parent < pEnd) {
    814:         q = parent;

    816:         PetscCall(DMPlexGetTrueSupportSize(dm, q, &qdof, numTrueSupp));
    817:         PetscCall(PetscSectionAddDof(newSupportSection, p, qdof));
    818:         PetscCall(PetscSectionAddDof(newSupportSection, q, dof));
    819:         PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
    820:       }
    821:     }
    822:   }
    823:   PetscCall(PetscSectionSetUp(newSupportSection));
    824:   PetscCall(PetscSectionGetStorageSize(newSupportSection, &newSize));
    825:   PetscCall(PetscMalloc1(newSize, &newSupports));
    826:   for (d = 0; d <= depth; d++) {
    827:     PetscCall(DMPlexGetHeightStratum(dm, d, &pStart, &pEnd));
    828:     for (p = pStart; p < pEnd; p++) {
    829:       PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;

    831:       PetscCall(PetscSectionGetDof(mesh->supportSection, p, &dof));
    832:       PetscCall(PetscSectionGetOffset(mesh->supportSection, p, &off));
    833:       PetscCall(PetscSectionGetDof(newSupportSection, p, &newDof));
    834:       PetscCall(PetscSectionGetOffset(newSupportSection, p, &newOff));
    835:       for (i = 0; i < dof; i++) {
    836:         PetscInt        numCones, j;
    837:         const PetscInt *cone;
    838:         PetscInt        q = mesh->supports[off + i];

    840:         PetscCall(DMPlexGetConeSize(dm, q, &numCones));
    841:         PetscCall(DMPlexGetCone(dm, q, &cone));
    842:         for (j = 0; j < numCones; j++) {
    843:           if (cone[j] == p) break;
    844:         }
    845:         if (j < numCones) newSupports[newOff + offsets[p]++] = q;
    846:       }

    848:       q = p;
    849:       PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
    850:       while (parent != q && parent >= pStart && parent < pEnd) {
    851:         q = parent;
    852:         PetscCall(PetscSectionGetDof(mesh->supportSection, q, &qdof));
    853:         PetscCall(PetscSectionGetOffset(mesh->supportSection, q, &qoff));
    854:         PetscCall(PetscSectionGetOffset(newSupportSection, q, &newqOff));
    855:         for (i = 0; i < qdof; i++) {
    856:           PetscInt        numCones, j;
    857:           const PetscInt *cone;
    858:           PetscInt        r = mesh->supports[qoff + i];

    860:           PetscCall(DMPlexGetConeSize(dm, r, &numCones));
    861:           PetscCall(DMPlexGetCone(dm, r, &cone));
    862:           for (j = 0; j < numCones; j++) {
    863:             if (cone[j] == q) break;
    864:           }
    865:           if (j < numCones) newSupports[newOff + offsets[p]++] = r;
    866:         }
    867:         for (i = 0; i < dof; i++) {
    868:           PetscInt        numCones, j;
    869:           const PetscInt *cone;
    870:           PetscInt        r = mesh->supports[off + i];

    872:           PetscCall(DMPlexGetConeSize(dm, r, &numCones));
    873:           PetscCall(DMPlexGetCone(dm, r, &cone));
    874:           for (j = 0; j < numCones; j++) {
    875:             if (cone[j] == p) break;
    876:           }
    877:           if (j < numCones) newSupports[newqOff + offsets[q]++] = r;
    878:         }
    879:         PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
    880:       }
    881:     }
    882:   }
    883:   PetscCall(PetscSectionDestroy(&mesh->supportSection));
    884:   mesh->supportSection = newSupportSection;
    885:   PetscCall(PetscFree(mesh->supports));
    886:   mesh->supports = newSupports;
    887:   PetscCall(PetscFree(offsets));
    888:   PetscCall(PetscFree(numTrueSupp));
    889:   PetscFunctionReturn(PETSC_SUCCESS);
    890: }

    892: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM, PetscSection, PetscSection, Mat);
    893: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM, PetscSection, PetscSection, Mat);

    895: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
    896: {
    897:   DM_Plex *mesh = (DM_Plex *)dm->data;
    898:   DM       refTree;
    899:   PetscInt size;

    901:   PetscFunctionBegin;
    904:   PetscCall(PetscObjectReference((PetscObject)parentSection));
    905:   PetscCall(PetscSectionDestroy(&mesh->parentSection));
    906:   mesh->parentSection = parentSection;
    907:   PetscCall(PetscSectionGetStorageSize(parentSection, &size));
    908:   if (parents != mesh->parents) {
    909:     PetscCall(PetscFree(mesh->parents));
    910:     PetscCall(PetscMalloc1(size, &mesh->parents));
    911:     PetscCall(PetscArraycpy(mesh->parents, parents, size));
    912:   }
    913:   if (childIDs != mesh->childIDs) {
    914:     PetscCall(PetscFree(mesh->childIDs));
    915:     PetscCall(PetscMalloc1(size, &mesh->childIDs));
    916:     PetscCall(PetscArraycpy(mesh->childIDs, childIDs, size));
    917:   }
    918:   PetscCall(DMPlexGetReferenceTree(dm, &refTree));
    919:   if (refTree) {
    920:     DMLabel canonLabel;

    922:     PetscCall(DMGetLabel(refTree, "canonical", &canonLabel));
    923:     if (canonLabel) {
    924:       PetscInt i;

    926:       for (i = 0; i < size; i++) {
    927:         PetscInt canon;
    928:         PetscCall(DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon));
    929:         if (canon >= 0) mesh->childIDs[i] = canon;
    930:       }
    931:     }
    932:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
    933:   } else {
    934:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
    935:   }
    936:   PetscCall(DMPlexTreeSymmetrize(dm));
    937:   if (computeCanonical) {
    938:     PetscInt d, dim;

    940:     /* add the canonical label */
    941:     PetscCall(DMGetDimension(dm, &dim));
    942:     PetscCall(DMCreateLabel(dm, "canonical"));
    943:     for (d = 0; d <= dim; d++) {
    944:       PetscInt        p, dStart, dEnd, canon = -1, cNumChildren;
    945:       const PetscInt *cChildren;

    947:       PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
    948:       for (p = dStart; p < dEnd; p++) {
    949:         PetscCall(DMPlexGetTreeChildren(dm, p, &cNumChildren, &cChildren));
    950:         if (cNumChildren) {
    951:           canon = p;
    952:           break;
    953:         }
    954:       }
    955:       if (canon == -1) continue;
    956:       for (p = dStart; p < dEnd; p++) {
    957:         PetscInt        numChildren, i;
    958:         const PetscInt *children;

    960:         PetscCall(DMPlexGetTreeChildren(dm, p, &numChildren, &children));
    961:         if (numChildren) {
    962:           PetscCheck(numChildren == cNumChildren, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "All parent points in a stratum should have the same number of children: %" PetscInt_FMT " != %" PetscInt_FMT, numChildren, cNumChildren);
    963:           PetscCall(DMSetLabelValue(dm, "canonical", p, canon));
    964:           for (i = 0; i < numChildren; i++) PetscCall(DMSetLabelValue(dm, "canonical", children[i], cChildren[i]));
    965:         }
    966:       }
    967:     }
    968:   }
    969:   if (exchangeSupports) PetscCall(DMPlexTreeExchangeSupports(dm));
    970:   mesh->createanchors = DMPlexCreateAnchors_Tree;
    971:   /* reset anchors */
    972:   PetscCall(DMPlexSetAnchors(dm, NULL, NULL));
    973:   PetscFunctionReturn(PETSC_SUCCESS);
    974: }

    976: /*@
    977:   DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points.  This routine also creates
    978:   the point-to-point constraints determined by the tree: a point is constrained to the points in the closure of its
    979:   tree root.

    981:   Collective

    983:   Input Parameters:
    984: + dm            - the `DMPLEX` object
    985: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
    986:                   offset indexes the parent and childID list; the reference count of parentSection is incremented
    987: . parents       - a list of the point parents; copied, can be destroyed
    988: - childIDs      - identifies the relationship of the child point to the parent point; if there is a reference tree, then
    989:              the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed

    991:   Level: intermediate

    993: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()`
    994: @*/
    995: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
    996: {
    997:   PetscFunctionBegin;
    998:   PetscCall(DMPlexSetTree_Internal(dm, parentSection, parents, childIDs, PETSC_FALSE, PETSC_TRUE));
    999:   PetscFunctionReturn(PETSC_SUCCESS);
    1000: }

    1002: /*@
    1003:   DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
    1004:   Collective

    1006:   Input Parameter:
    1007: . dm - the `DMPLEX` object

    1009:   Output Parameters:
    1010: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
    1011:                   offset indexes the parent and childID list
    1012: . parents       - a list of the point parents
    1013: . childIDs      - identifies the relationship of the child point to the parent point; if there is a reference tree, then
    1014:              the child corresponds to the point in the reference tree with index childID
    1015: . childSection  - the inverse of the parent section
    1016: - children      - a list of the point children

    1018:   Level: intermediate

    1020: .seealso: [](ch_unstructured), `DM`, `DMPLEX`,`DMPlexSetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()`
    1021: @*/
    1022: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
    1023: {
    1024:   DM_Plex *mesh = (DM_Plex *)dm->data;

    1026:   PetscFunctionBegin;
    1028:   if (parentSection) *parentSection = mesh->parentSection;
    1029:   if (parents) *parents = mesh->parents;
    1030:   if (childIDs) *childIDs = mesh->childIDs;
    1031:   if (childSection) *childSection = mesh->childSection;
    1032:   if (children) *children = mesh->children;
    1033:   PetscFunctionReturn(PETSC_SUCCESS);
    1034: }

    1036: /*@
    1037:   DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the DAG)

    1039:   Input Parameters:
    1040: + dm    - the `DMPLEX` object
    1041: - point - the query point

    1043:   Output Parameters:
    1044: + parent  - if not `NULL`, set to the parent of the point, or the point itself if the point does not have a parent
    1045: - childID - if not `NULL`, set to the child ID of the point with respect to its parent, or 0 if the point
    1046:             does not have a parent

    1048:   Level: intermediate

    1050: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeChildren()`
    1051: @*/
    1052: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
    1053: {
    1054:   DM_Plex     *mesh = (DM_Plex *)dm->data;
    1055:   PetscSection pSec;

    1057:   PetscFunctionBegin;
    1059:   pSec = mesh->parentSection;
    1060:   if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
    1061:     PetscInt dof;

    1063:     PetscCall(PetscSectionGetDof(pSec, point, &dof));
    1064:     if (dof) {
    1065:       PetscInt off;

    1067:       PetscCall(PetscSectionGetOffset(pSec, point, &off));
    1068:       if (parent) *parent = mesh->parents[off];
    1069:       if (childID) *childID = mesh->childIDs[off];
    1070:       PetscFunctionReturn(PETSC_SUCCESS);
    1071:     }
    1072:   }
    1073:   if (parent) *parent = point;
    1074:   if (childID) *childID = 0;
    1075:   PetscFunctionReturn(PETSC_SUCCESS);
    1076: }

    1078: /*@C
    1079:   DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the DAG)

    1081:   Input Parameters:
    1082: + dm    - the `DMPLEX` object
    1083: - point - the query point

    1085:   Output Parameters:
    1086: + numChildren - if not `NULL`, set to the number of children
    1087: - children    - if not `NULL`, set to a list children, or set to `NULL` if the point has no children

    1089:   Level: intermediate

    1091: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeParent()`
    1092: @*/
    1093: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
    1094: {
    1095:   DM_Plex     *mesh = (DM_Plex *)dm->data;
    1096:   PetscSection childSec;
    1097:   PetscInt     dof = 0;

    1099:   PetscFunctionBegin;
    1101:   childSec = mesh->childSection;
    1102:   if (childSec && point >= childSec->pStart && point < childSec->pEnd) PetscCall(PetscSectionGetDof(childSec, point, &dof));
    1103:   if (numChildren) *numChildren = dof;
    1104:   if (children) {
    1105:     if (dof) {
    1106:       PetscInt off;

    1108:       PetscCall(PetscSectionGetOffset(childSec, point, &off));
    1109:       *children = &mesh->children[off];
    1110:     } else {
    1111:       *children = NULL;
    1112:     }
    1113:   }
    1114:   PetscFunctionReturn(PETSC_SUCCESS);
    1115: }

    1117: static PetscErrorCode EvaluateBasis(PetscSpace space, PetscInt nBasis, PetscInt nFunctionals, PetscInt nComps, PetscInt nPoints, const PetscInt *pointsPerFn, const PetscReal *points, const PetscReal *weights, PetscReal *work, Mat basisAtPoints)
    1118: {
    1119:   PetscInt f, b, p, c, offset, qPoints;

    1121:   PetscFunctionBegin;
    1122:   PetscCall(PetscSpaceEvaluate(space, nPoints, points, work, NULL, NULL));
    1123:   for (f = 0, offset = 0; f < nFunctionals; f++) {
    1124:     qPoints = pointsPerFn[f];
    1125:     for (b = 0; b < nBasis; b++) {
    1126:       PetscScalar val = 0.;

    1128:       for (p = 0; p < qPoints; p++) {
    1129:         for (c = 0; c < nComps; c++) val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c];
    1130:       }
    1131:       PetscCall(MatSetValue(basisAtPoints, b, f, val, INSERT_VALUES));
    1132:     }
    1133:     offset += qPoints;
    1134:   }
    1135:   PetscCall(MatAssemblyBegin(basisAtPoints, MAT_FINAL_ASSEMBLY));
    1136:   PetscCall(MatAssemblyEnd(basisAtPoints, MAT_FINAL_ASSEMBLY));
    1137:   PetscFunctionReturn(PETSC_SUCCESS);
    1138: }

    1140: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
    1141: {
    1142:   PetscDS         ds;
    1143:   PetscInt        spdim;
    1144:   PetscInt        numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
    1145:   const PetscInt *anchors;
    1146:   PetscSection    aSec;
    1147:   PetscReal      *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
    1148:   IS              aIS;

    1150:   PetscFunctionBegin;
    1151:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    1152:   PetscCall(DMGetDS(dm, &ds));
    1153:   PetscCall(PetscDSGetNumFields(ds, &numFields));
    1154:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    1155:   PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS));
    1156:   PetscCall(ISGetIndices(aIS, &anchors));
    1157:   PetscCall(PetscSectionGetChart(cSec, &conStart, &conEnd));
    1158:   PetscCall(DMGetDimension(dm, &spdim));
    1159:   PetscCall(PetscMalloc6(spdim, &v0, spdim, &v0parent, spdim, &vtmp, spdim * spdim, &J, spdim * spdim, &Jparent, spdim * spdim, &invJparent));

    1161:   for (f = 0; f < numFields; f++) {
    1162:     PetscObject          disc;
    1163:     PetscClassId         id;
    1164:     PetscSpace           bspace;
    1165:     PetscDualSpace       dspace;
    1166:     PetscInt             i, j, k, nPoints, Nc, offset;
    1167:     PetscInt             fSize, maxDof;
    1168:     PetscReal           *weights, *pointsRef, *pointsReal, *work;
    1169:     PetscScalar         *scwork;
    1170:     const PetscScalar   *X;
    1171:     PetscInt            *sizes, *workIndRow, *workIndCol;
    1172:     Mat                  Amat, Bmat, Xmat;
    1173:     const PetscInt      *numDof = NULL;
    1174:     const PetscInt    ***perms  = NULL;
    1175:     const PetscScalar ***flips  = NULL;

    1177:     PetscCall(PetscDSGetDiscretization(ds, f, &disc));
    1178:     PetscCall(PetscObjectGetClassId(disc, &id));
    1179:     if (id == PETSCFE_CLASSID) {
    1180:       PetscFE fe = (PetscFE)disc;

    1182:       PetscCall(PetscFEGetBasisSpace(fe, &bspace));
    1183:       PetscCall(PetscFEGetDualSpace(fe, &dspace));
    1184:       PetscCall(PetscDualSpaceGetDimension(dspace, &fSize));
    1185:       PetscCall(PetscFEGetNumComponents(fe, &Nc));
    1186:     } else if (id == PETSCFV_CLASSID) {
    1187:       PetscFV fv = (PetscFV)disc;

    1189:       PetscCall(PetscFVGetNumComponents(fv, &Nc));
    1190:       PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)fv), &bspace));
    1191:       PetscCall(PetscSpaceSetType(bspace, PETSCSPACEPOLYNOMIAL));
    1192:       PetscCall(PetscSpaceSetDegree(bspace, 0, PETSC_DETERMINE));
    1193:       PetscCall(PetscSpaceSetNumComponents(bspace, Nc));
    1194:       PetscCall(PetscSpaceSetNumVariables(bspace, spdim));
    1195:       PetscCall(PetscSpaceSetUp(bspace));
    1196:       PetscCall(PetscFVGetDualSpace(fv, &dspace));
    1197:       PetscCall(PetscDualSpaceGetDimension(dspace, &fSize));
    1198:     } else SETERRQ(PetscObjectComm(disc), PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
    1199:     PetscCall(PetscDualSpaceGetNumDof(dspace, &numDof));
    1200:     for (i = 0, maxDof = 0; i <= spdim; i++) maxDof = PetscMax(maxDof, numDof[i]);
    1201:     PetscCall(PetscDualSpaceGetSymmetries(dspace, &perms, &flips));

    1203:     PetscCall(MatCreate(PETSC_COMM_SELF, &Amat));
    1204:     PetscCall(MatSetSizes(Amat, fSize, fSize, fSize, fSize));
    1205:     PetscCall(MatSetType(Amat, MATSEQDENSE));
    1206:     PetscCall(MatSetUp(Amat));
    1207:     PetscCall(MatDuplicate(Amat, MAT_DO_NOT_COPY_VALUES, &Bmat));
    1208:     PetscCall(MatDuplicate(Amat, MAT_DO_NOT_COPY_VALUES, &Xmat));
    1209:     nPoints = 0;
    1210:     for (i = 0; i < fSize; i++) {
    1211:       PetscInt        qPoints, thisNc;
    1212:       PetscQuadrature quad;

    1214:       PetscCall(PetscDualSpaceGetFunctional(dspace, i, &quad));
    1215:       PetscCall(PetscQuadratureGetData(quad, NULL, &thisNc, &qPoints, NULL, NULL));
    1216:       PetscCheck(thisNc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT, thisNc, Nc);
    1217:       nPoints += qPoints;
    1218:     }
    1219:     PetscCall(PetscMalloc7(fSize, &sizes, nPoints * Nc, &weights, spdim * nPoints, &pointsRef, spdim * nPoints, &pointsReal, nPoints * fSize * Nc, &work, maxDof, &workIndRow, maxDof, &workIndCol));
    1220:     PetscCall(PetscMalloc1(maxDof * maxDof, &scwork));
    1221:     offset = 0;
    1222:     for (i = 0; i < fSize; i++) {
    1223:       PetscInt         qPoints;
    1224:       const PetscReal *p, *w;
    1225:       PetscQuadrature  quad;

    1227:       PetscCall(PetscDualSpaceGetFunctional(dspace, i, &quad));
    1228:       PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &qPoints, &p, &w));
    1229:       PetscCall(PetscArraycpy(weights + Nc * offset, w, Nc * qPoints));
    1230:       PetscCall(PetscArraycpy(pointsRef + spdim * offset, p, spdim * qPoints));
    1231:       sizes[i] = qPoints;
    1232:       offset += qPoints;
    1233:     }
    1234:     PetscCall(EvaluateBasis(bspace, fSize, fSize, Nc, nPoints, sizes, pointsRef, weights, work, Amat));
    1235:     PetscCall(MatLUFactor(Amat, NULL, NULL, NULL));
    1236:     for (c = cStart; c < cEnd; c++) {
    1237:       PetscInt  parent;
    1238:       PetscInt  closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
    1239:       PetscInt *childOffsets, *parentOffsets;

    1241:       PetscCall(DMPlexGetTreeParent(dm, c, &parent, NULL));
    1242:       if (parent == c) continue;
    1243:       PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    1244:       for (i = 0; i < closureSize; i++) {
    1245:         PetscInt p = closure[2 * i];
    1246:         PetscInt conDof;

    1248:         if (p < conStart || p >= conEnd) continue;
    1249:         if (numFields) {
    1250:           PetscCall(PetscSectionGetFieldDof(cSec, p, f, &conDof));
    1251:         } else {
    1252:           PetscCall(PetscSectionGetDof(cSec, p, &conDof));
    1253:         }
    1254:         if (conDof) break;
    1255:       }
    1256:       if (i == closureSize) {
    1257:         PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    1258:         continue;
    1259:       }

    1261:       PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ));
    1262:       PetscCall(DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent));
    1263:       for (i = 0; i < nPoints; i++) {
    1264:         const PetscReal xi0[3] = {-1., -1., -1.};

    1266:         CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i * spdim], vtmp);
    1267:         CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i * spdim]);
    1268:       }
    1269:       PetscCall(EvaluateBasis(bspace, fSize, fSize, Nc, nPoints, sizes, pointsReal, weights, work, Bmat));
    1270:       PetscCall(MatMatSolve(Amat, Bmat, Xmat));
    1271:       PetscCall(MatDenseGetArrayRead(Xmat, &X));
    1272:       PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSizeP, &closureP));
    1273:       PetscCall(PetscMalloc2(closureSize + 1, &childOffsets, closureSizeP + 1, &parentOffsets));
    1274:       childOffsets[0] = 0;
    1275:       for (i = 0; i < closureSize; i++) {
    1276:         PetscInt p = closure[2 * i];
    1277:         PetscInt dof;

    1279:         if (numFields) {
    1280:           PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
    1281:         } else {
    1282:           PetscCall(PetscSectionGetDof(section, p, &dof));
    1283:         }
    1284:         childOffsets[i + 1] = childOffsets[i] + dof;
    1285:       }
    1286:       parentOffsets[0] = 0;
    1287:       for (i = 0; i < closureSizeP; i++) {
    1288:         PetscInt p = closureP[2 * i];
    1289:         PetscInt dof;

    1291:         if (numFields) {
    1292:           PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
    1293:         } else {
    1294:           PetscCall(PetscSectionGetDof(section, p, &dof));
    1295:         }
    1296:         parentOffsets[i + 1] = parentOffsets[i] + dof;
    1297:       }
    1298:       for (i = 0; i < closureSize; i++) {
    1299:         PetscInt           conDof, conOff, aDof, aOff, nWork;
    1300:         PetscInt           p = closure[2 * i];
    1301:         PetscInt           o = closure[2 * i + 1];
    1302:         const PetscInt    *perm;
    1303:         const PetscScalar *flip;

    1305:         if (p < conStart || p >= conEnd) continue;
    1306:         if (numFields) {
    1307:           PetscCall(PetscSectionGetFieldDof(cSec, p, f, &conDof));
    1308:           PetscCall(PetscSectionGetFieldOffset(cSec, p, f, &conOff));
    1309:         } else {
    1310:           PetscCall(PetscSectionGetDof(cSec, p, &conDof));
    1311:           PetscCall(PetscSectionGetOffset(cSec, p, &conOff));
    1312:         }
    1313:         if (!conDof) continue;
    1314:         perm = (perms && perms[i]) ? perms[i][o] : NULL;
    1315:         flip = (flips && flips[i]) ? flips[i][o] : NULL;
    1316:         PetscCall(PetscSectionGetDof(aSec, p, &aDof));
    1317:         PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
    1318:         nWork = childOffsets[i + 1] - childOffsets[i];
    1319:         for (k = 0; k < aDof; k++) {
    1320:           PetscInt a = anchors[aOff + k];
    1321:           PetscInt aSecDof, aSecOff;

    1323:           if (numFields) {
    1324:             PetscCall(PetscSectionGetFieldDof(section, a, f, &aSecDof));
    1325:             PetscCall(PetscSectionGetFieldOffset(section, a, f, &aSecOff));
    1326:           } else {
    1327:             PetscCall(PetscSectionGetDof(section, a, &aSecDof));
    1328:             PetscCall(PetscSectionGetOffset(section, a, &aSecOff));
    1329:           }
    1330:           if (!aSecDof) continue;

    1332:           for (j = 0; j < closureSizeP; j++) {
    1333:             PetscInt q  = closureP[2 * j];
    1334:             PetscInt oq = closureP[2 * j + 1];

    1336:             if (q == a) {
    1337:               PetscInt           r, s, nWorkP;
    1338:               const PetscInt    *permP;
    1339:               const PetscScalar *flipP;

    1341:               permP  = (perms && perms[j]) ? perms[j][oq] : NULL;
    1342:               flipP  = (flips && flips[j]) ? flips[j][oq] : NULL;
    1343:               nWorkP = parentOffsets[j + 1] - parentOffsets[j];
    1344:               /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the
    1345:                * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArrayRead() is
    1346:                * column-major, so transpose-transpose = do nothing */
    1347:               for (r = 0; r < nWork; r++) {
    1348:                 for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])];
    1349:               }
    1350:               for (r = 0; r < nWork; r++) workIndRow[perm ? perm[r] : r] = conOff + r;
    1351:               for (s = 0; s < nWorkP; s++) workIndCol[permP ? permP[s] : s] = aSecOff + s;
    1352:               if (flip) {
    1353:                 for (r = 0; r < nWork; r++) {
    1354:                   for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] *= flip[r];
    1355:                 }
    1356:               }
    1357:               if (flipP) {
    1358:                 for (r = 0; r < nWork; r++) {
    1359:                   for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] *= flipP[s];
    1360:                 }
    1361:               }
    1362:               PetscCall(MatSetValues(cMat, nWork, workIndRow, nWorkP, workIndCol, scwork, INSERT_VALUES));
    1363:               break;
    1364:             }
    1365:           }
    1366:         }
    1367:       }
    1368:       PetscCall(MatDenseRestoreArrayRead(Xmat, &X));
    1369:       PetscCall(PetscFree2(childOffsets, parentOffsets));
    1370:       PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    1371:       PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSizeP, &closureP));
    1372:     }
    1373:     PetscCall(MatDestroy(&Amat));
    1374:     PetscCall(MatDestroy(&Bmat));
    1375:     PetscCall(MatDestroy(&Xmat));
    1376:     PetscCall(PetscFree(scwork));
    1377:     PetscCall(PetscFree7(sizes, weights, pointsRef, pointsReal, work, workIndRow, workIndCol));
    1378:     if (id == PETSCFV_CLASSID) PetscCall(PetscSpaceDestroy(&bspace));
    1379:   }
    1380:   PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
    1381:   PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
    1382:   PetscCall(PetscFree6(v0, v0parent, vtmp, J, Jparent, invJparent));
    1383:   PetscCall(ISRestoreIndices(aIS, &anchors));
    1384:   PetscFunctionReturn(PETSC_SUCCESS);
    1385: }

    1387: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
    1388: {
    1389:   Mat                 refCmat;
    1390:   PetscDS             ds;
    1391:   PetscInt            numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
    1392:   PetscScalar      ***refPointFieldMats;
    1393:   PetscSection        refConSec, refAnSec, refSection;
    1394:   IS                  refAnIS;
    1395:   const PetscInt     *refAnchors;
    1396:   const PetscInt    **perms;
    1397:   const PetscScalar **flips;

    1399:   PetscFunctionBegin;
    1400:   PetscCall(DMGetDS(refTree, &ds));
    1401:   PetscCall(PetscDSGetNumFields(ds, &numFields));
    1402:   maxFields = PetscMax(1, numFields);
    1403:   PetscCall(DMGetDefaultConstraints(refTree, &refConSec, &refCmat, NULL));
    1404:   PetscCall(DMPlexGetAnchors(refTree, &refAnSec, &refAnIS));
    1405:   PetscCall(ISGetIndices(refAnIS, &refAnchors));
    1406:   PetscCall(DMGetLocalSection(refTree, &refSection));
    1407:   PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
    1408:   PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats));
    1409:   PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldN));
    1410:   PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
    1411:   PetscCall(PetscSectionGetMaxDof(refAnSec, &maxAnDof));
    1412:   PetscCall(PetscMalloc1(maxDof, &rows));
    1413:   PetscCall(PetscMalloc1(maxDof * maxAnDof, &cols));
    1414:   for (p = pRefStart; p < pRefEnd; p++) {
    1415:     PetscInt parent, closureSize, *closure = NULL, pDof;

    1417:     PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
    1418:     PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
    1419:     if (!pDof || parent == p) continue;

    1421:     PetscCall(PetscMalloc1(maxFields, &refPointFieldMats[p - pRefStart]));
    1422:     PetscCall(PetscCalloc1(maxFields, &refPointFieldN[p - pRefStart]));
    1423:     PetscCall(DMPlexGetTransitiveClosure(refTree, parent, PETSC_TRUE, &closureSize, &closure));
    1424:     for (f = 0; f < maxFields; f++) {
    1425:       PetscInt cDof, cOff, numCols, r, i;

    1427:       if (f < numFields) {
    1428:         PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
    1429:         PetscCall(PetscSectionGetFieldOffset(refConSec, p, f, &cOff));
    1430:         PetscCall(PetscSectionGetFieldPointSyms(refSection, f, closureSize, closure, &perms, &flips));
    1431:       } else {
    1432:         PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
    1433:         PetscCall(PetscSectionGetOffset(refConSec, p, &cOff));
    1434:         PetscCall(PetscSectionGetPointSyms(refSection, closureSize, closure, &perms, &flips));
    1435:       }

    1437:       for (r = 0; r < cDof; r++) rows[r] = cOff + r;
    1438:       numCols = 0;
    1439:       for (i = 0; i < closureSize; i++) {
    1440:         PetscInt        q = closure[2 * i];
    1441:         PetscInt        aDof, aOff, j;
    1442:         const PetscInt *perm = perms ? perms[i] : NULL;

    1444:         if (numFields) {
    1445:           PetscCall(PetscSectionGetFieldDof(refSection, q, f, &aDof));
    1446:           PetscCall(PetscSectionGetFieldOffset(refSection, q, f, &aOff));
    1447:         } else {
    1448:           PetscCall(PetscSectionGetDof(refSection, q, &aDof));
    1449:           PetscCall(PetscSectionGetOffset(refSection, q, &aOff));
    1450:         }

    1452:         for (j = 0; j < aDof; j++) cols[numCols++] = aOff + (perm ? perm[j] : j);
    1453:       }
    1454:       refPointFieldN[p - pRefStart][f] = numCols;
    1455:       PetscCall(PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]));
    1456:       PetscCall(MatGetValues(refCmat, cDof, rows, numCols, cols, refPointFieldMats[p - pRefStart][f]));
    1457:       if (flips) {
    1458:         PetscInt colOff = 0;

    1460:         for (i = 0; i < closureSize; i++) {
    1461:           PetscInt           q = closure[2 * i];
    1462:           PetscInt           aDof, aOff, j;
    1463:           const PetscScalar *flip = flips ? flips[i] : NULL;

    1465:           if (numFields) {
    1466:             PetscCall(PetscSectionGetFieldDof(refSection, q, f, &aDof));
    1467:             PetscCall(PetscSectionGetFieldOffset(refSection, q, f, &aOff));
    1468:           } else {
    1469:             PetscCall(PetscSectionGetDof(refSection, q, &aDof));
    1470:             PetscCall(PetscSectionGetOffset(refSection, q, &aOff));
    1471:           }
    1472:           if (flip) {
    1473:             PetscInt k;
    1474:             for (k = 0; k < cDof; k++) {
    1475:               for (j = 0; j < aDof; j++) refPointFieldMats[p - pRefStart][f][k * numCols + colOff + j] *= flip[j];
    1476:             }
    1477:           }
    1478:           colOff += aDof;
    1479:         }
    1480:       }
    1481:       if (numFields) {
    1482:         PetscCall(PetscSectionRestoreFieldPointSyms(refSection, f, closureSize, closure, &perms, &flips));
    1483:       } else {
    1484:         PetscCall(PetscSectionRestorePointSyms(refSection, closureSize, closure, &perms, &flips));
    1485:       }
    1486:     }
    1487:     PetscCall(DMPlexRestoreTransitiveClosure(refTree, parent, PETSC_TRUE, &closureSize, &closure));
    1488:   }
    1489:   *childrenMats = refPointFieldMats;
    1490:   *childrenN    = refPointFieldN;
    1491:   PetscCall(ISRestoreIndices(refAnIS, &refAnchors));
    1492:   PetscCall(PetscFree(rows));
    1493:   PetscCall(PetscFree(cols));
    1494:   PetscFunctionReturn(PETSC_SUCCESS);
    1495: }

    1497: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
    1498: {
    1499:   PetscDS        ds;
    1500:   PetscInt     **refPointFieldN;
    1501:   PetscScalar ***refPointFieldMats;
    1502:   PetscInt       numFields, maxFields, pRefStart, pRefEnd, p, f;
    1503:   PetscSection   refConSec;

    1505:   PetscFunctionBegin;
    1506:   refPointFieldN    = *childrenN;
    1507:   *childrenN        = NULL;
    1508:   refPointFieldMats = *childrenMats;
    1509:   *childrenMats     = NULL;
    1510:   PetscCall(DMGetDS(refTree, &ds));
    1511:   PetscCall(PetscDSGetNumFields(ds, &numFields));
    1512:   maxFields = PetscMax(1, numFields);
    1513:   PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
    1514:   PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
    1515:   for (p = pRefStart; p < pRefEnd; p++) {
    1516:     PetscInt parent, pDof;

    1518:     PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
    1519:     PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
    1520:     if (!pDof || parent == p) continue;

    1522:     for (f = 0; f < maxFields; f++) {
    1523:       PetscInt cDof;

    1525:       if (numFields) {
    1526:         PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
    1527:       } else {
    1528:         PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
    1529:       }

    1531:       PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f]));
    1532:     }
    1533:     PetscCall(PetscFree(refPointFieldMats[p - pRefStart]));
    1534:     PetscCall(PetscFree(refPointFieldN[p - pRefStart]));
    1535:   }
    1536:   PetscCall(PetscFree(refPointFieldMats));
    1537:   PetscCall(PetscFree(refPointFieldN));
    1538:   PetscFunctionReturn(PETSC_SUCCESS);
    1539: }

    1541: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
    1542: {
    1543:   DM              refTree;
    1544:   PetscDS         ds;
    1545:   Mat             refCmat;
    1546:   PetscInt        numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
    1547:   PetscScalar  ***refPointFieldMats, *pointWork;
    1548:   PetscSection    refConSec, refAnSec, anSec;
    1549:   IS              refAnIS, anIS;
    1550:   const PetscInt *anchors;

    1552:   PetscFunctionBegin;
    1554:   PetscCall(DMGetDS(dm, &ds));
    1555:   PetscCall(PetscDSGetNumFields(ds, &numFields));
    1556:   maxFields = PetscMax(1, numFields);
    1557:   PetscCall(DMPlexGetReferenceTree(dm, &refTree));
    1558:   PetscCall(DMCopyDisc(dm, refTree));
    1559:   PetscCall(DMSetLocalSection(refTree, NULL));
    1560:   PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
    1561:   PetscCall(DMGetDefaultConstraints(refTree, &refConSec, &refCmat, NULL));
    1562:   PetscCall(DMPlexGetAnchors(refTree, &refAnSec, &refAnIS));
    1563:   PetscCall(DMPlexGetAnchors(dm, &anSec, &anIS));
    1564:   PetscCall(ISGetIndices(anIS, &anchors));
    1565:   PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
    1566:   PetscCall(PetscSectionGetChart(conSec, &conStart, &conEnd));
    1567:   PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
    1568:   PetscCall(PetscSectionGetMaxDof(refAnSec, &maxAnDof));
    1569:   PetscCall(PetscMalloc1(maxDof * maxDof * maxAnDof, &pointWork));

    1571:   /* step 1: get submats for every constrained point in the reference tree */
    1572:   PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));

    1574:   /* step 2: compute the preorder */
    1575:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    1576:   PetscCall(PetscMalloc2(pEnd - pStart, &perm, pEnd - pStart, &iperm));
    1577:   for (p = pStart; p < pEnd; p++) {
    1578:     perm[p - pStart]  = p;
    1579:     iperm[p - pStart] = p - pStart;
    1580:   }
    1581:   for (p = 0; p < pEnd - pStart;) {
    1582:     PetscInt point = perm[p];
    1583:     PetscInt parent;

    1585:     PetscCall(DMPlexGetTreeParent(dm, point, &parent, NULL));
    1586:     if (parent == point) {
    1587:       p++;
    1588:     } else {
    1589:       PetscInt size, closureSize, *closure = NULL, i;

    1591:       PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
    1592:       for (i = 0; i < closureSize; i++) {
    1593:         PetscInt q = closure[2 * i];
    1594:         if (iperm[q - pStart] > iperm[point - pStart]) {
    1595:           /* swap */
    1596:           perm[p]                 = q;
    1597:           perm[iperm[q - pStart]] = point;
    1598:           iperm[point - pStart]   = iperm[q - pStart];
    1599:           iperm[q - pStart]       = p;
    1600:           break;
    1601:         }
    1602:       }
    1603:       size = closureSize;
    1604:       PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
    1605:       if (i == size) p++;
    1606:     }
    1607:   }

    1609:   /* step 3: fill the constraint matrix */
    1610:   /* we are going to use a preorder progressive fill strategy.  Mat doesn't
    1611:    * allow progressive fill without assembly, so we are going to set up the
    1612:    * values outside of the Mat first.
    1613:    */
    1614:   {
    1615:     PetscInt        nRows, row, nnz;
    1616:     PetscBool       done;
    1617:     PetscInt        secStart, secEnd;
    1618:     const PetscInt *ia, *ja;
    1619:     PetscScalar    *vals;

    1621:     PetscCall(PetscSectionGetChart(section, &secStart, &secEnd));
    1622:     PetscCall(MatGetRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done));
    1623:     PetscCheck(done, PetscObjectComm((PetscObject)cMat), PETSC_ERR_PLIB, "Could not get RowIJ of constraint matrix");
    1624:     nnz = ia[nRows];
    1625:     /* malloc and then zero rows right before we fill them: this way valgrind
    1626:      * can tell if we are doing progressive fill in the wrong order */
    1627:     PetscCall(PetscMalloc1(nnz, &vals));
    1628:     for (p = 0; p < pEnd - pStart; p++) {
    1629:       PetscInt parent, childid, closureSize, *closure = NULL;
    1630:       PetscInt point = perm[p], pointDof;

    1632:       PetscCall(DMPlexGetTreeParent(dm, point, &parent, &childid));
    1633:       if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
    1634:       PetscCall(PetscSectionGetDof(conSec, point, &pointDof));
    1635:       if (!pointDof) continue;
    1636:       PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
    1637:       for (f = 0; f < maxFields; f++) {
    1638:         PetscInt            cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
    1639:         PetscScalar        *pointMat;
    1640:         const PetscInt    **perms;
    1641:         const PetscScalar **flips;

    1643:         if (numFields) {
    1644:           PetscCall(PetscSectionGetFieldDof(conSec, point, f, &cDof));
    1645:           PetscCall(PetscSectionGetFieldOffset(conSec, point, f, &cOff));
    1646:         } else {
    1647:           PetscCall(PetscSectionGetDof(conSec, point, &cDof));
    1648:           PetscCall(PetscSectionGetOffset(conSec, point, &cOff));
    1649:         }
    1650:         if (!cDof) continue;
    1651:         if (numFields) PetscCall(PetscSectionGetFieldPointSyms(section, f, closureSize, closure, &perms, &flips));
    1652:         else PetscCall(PetscSectionGetPointSyms(section, closureSize, closure, &perms, &flips));

    1654:         /* make sure that every row for this point is the same size */
    1655:         if (PetscDefined(USE_DEBUG)) {
    1656:           for (r = 0; r < cDof; r++) {
    1657:             if (cDof > 1 && r) {
    1658:               PetscCheck((ia[cOff + r + 1] - ia[cOff + r]) == (ia[cOff + r] - ia[cOff + r - 1]), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two point rows have different nnz: %" PetscInt_FMT " vs. %" PetscInt_FMT, ia[cOff + r + 1] - ia[cOff + r], ia[cOff + r] - ia[cOff + r - 1]);
    1659:             }
    1660:           }
    1661:         }
    1662:         /* zero rows */
    1663:         for (i = ia[cOff]; i < ia[cOff + cDof]; i++) vals[i] = 0.;
    1664:         matOffset   = ia[cOff];
    1665:         numFillCols = ia[cOff + 1] - matOffset;
    1666:         pointMat    = refPointFieldMats[childid - pRefStart][f];
    1667:         numCols     = refPointFieldN[childid - pRefStart][f];
    1668:         offset      = 0;
    1669:         for (i = 0; i < closureSize; i++) {
    1670:           PetscInt           q = closure[2 * i];
    1671:           PetscInt           aDof, aOff, j, k, qConDof, qConOff;
    1672:           const PetscInt    *perm = perms ? perms[i] : NULL;
    1673:           const PetscScalar *flip = flips ? flips[i] : NULL;

    1675:           qConDof = qConOff = 0;
    1676:           if (q < secStart || q >= secEnd) continue;
    1677:           if (numFields) {
    1678:             PetscCall(PetscSectionGetFieldDof(section, q, f, &aDof));
    1679:             PetscCall(PetscSectionGetFieldOffset(section, q, f, &aOff));
    1680:             if (q >= conStart && q < conEnd) {
    1681:               PetscCall(PetscSectionGetFieldDof(conSec, q, f, &qConDof));
    1682:               PetscCall(PetscSectionGetFieldOffset(conSec, q, f, &qConOff));
    1683:             }
    1684:           } else {
    1685:             PetscCall(PetscSectionGetDof(section, q, &aDof));
    1686:             PetscCall(PetscSectionGetOffset(section, q, &aOff));
    1687:             if (q >= conStart && q < conEnd) {
    1688:               PetscCall(PetscSectionGetDof(conSec, q, &qConDof));
    1689:               PetscCall(PetscSectionGetOffset(conSec, q, &qConOff));
    1690:             }
    1691:           }
    1692:           if (!aDof) continue;
    1693:           if (qConDof) {
    1694:             /* this point has anchors: its rows of the matrix should already
    1695:              * be filled, thanks to preordering */
    1696:             /* first multiply into pointWork, then set in matrix */
    1697:             PetscInt aMatOffset   = ia[qConOff];
    1698:             PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
    1699:             for (r = 0; r < cDof; r++) {
    1700:               for (j = 0; j < aNumFillCols; j++) {
    1701:                 PetscScalar inVal = 0;
    1702:                 for (k = 0; k < aDof; k++) {
    1703:                   PetscInt col = perm ? perm[k] : k;

    1705:                   inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.);
    1706:                 }
    1707:                 pointWork[r * aNumFillCols + j] = inVal;
    1708:               }
    1709:             }
    1710:             /* assume that the columns are sorted, spend less time searching */
    1711:             for (j = 0, k = 0; j < aNumFillCols; j++) {
    1712:               PetscInt col = ja[aMatOffset + j];
    1713:               for (; k < numFillCols; k++) {
    1714:                 if (ja[matOffset + k] == col) break;
    1715:               }
    1716:               PetscCheck(k != numFillCols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, col);
    1717:               for (r = 0; r < cDof; r++) vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
    1718:             }
    1719:           } else {
    1720:             /* find where to put this portion of pointMat into the matrix */
    1721:             for (k = 0; k < numFillCols; k++) {
    1722:               if (ja[matOffset + k] == aOff) break;
    1723:             }
    1724:             PetscCheck(k != numFillCols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, aOff);
    1725:             for (r = 0; r < cDof; r++) {
    1726:               for (j = 0; j < aDof; j++) {
    1727:                 PetscInt col = perm ? perm[j] : j;

    1729:                 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
    1730:               }
    1731:             }
    1732:           }
    1733:           offset += aDof;
    1734:         }
    1735:         if (numFields) {
    1736:           PetscCall(PetscSectionRestoreFieldPointSyms(section, f, closureSize, closure, &perms, &flips));
    1737:         } else {
    1738:           PetscCall(PetscSectionRestorePointSyms(section, closureSize, closure, &perms, &flips));
    1739:         }
    1740:       }
    1741:       PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
    1742:     }
    1743:     for (row = 0; row < nRows; row++) PetscCall(MatSetValues(cMat, 1, &row, ia[row + 1] - ia[row], &ja[ia[row]], &vals[ia[row]], INSERT_VALUES));
    1744:     PetscCall(MatRestoreRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done));
    1745:     PetscCheck(done, PetscObjectComm((PetscObject)cMat), PETSC_ERR_PLIB, "Could not restore RowIJ of constraint matrix");
    1746:     PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
    1747:     PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
    1748:     PetscCall(PetscFree(vals));
    1749:   }

    1751:   /* clean up */
    1752:   PetscCall(ISRestoreIndices(anIS, &anchors));
    1753:   PetscCall(PetscFree2(perm, iperm));
    1754:   PetscCall(PetscFree(pointWork));
    1755:   PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
    1756:   PetscFunctionReturn(PETSC_SUCCESS);
    1757: }

    1759: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
    1760:  * a non-conforming mesh.  Local refinement comes later */
    1761: PetscErrorCode DMPlexTreeRefineCell(DM dm, PetscInt cell, DM *ncdm)
    1762: {
    1763:   DM           K;
    1764:   PetscMPIInt  rank;
    1765:   PetscInt     dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
    1766:   PetscInt     numNewCones, *newConeSizes, *newCones, *newOrientations;
    1767:   PetscInt    *Kembedding;
    1768:   PetscInt    *cellClosure = NULL, nc;
    1769:   PetscScalar *newVertexCoords;
    1770:   PetscInt     numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
    1771:   PetscSection parentSection;

    1773:   PetscFunctionBegin;
    1774:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    1775:   PetscCall(DMGetDimension(dm, &dim));
    1776:   PetscCall(DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm));
    1777:   PetscCall(DMSetDimension(*ncdm, dim));

    1779:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    1780:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &parentSection));
    1781:   PetscCall(DMPlexGetReferenceTree(dm, &K));
    1782:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
    1783:   if (rank == 0) {
    1784:     /* compute the new charts */
    1785:     PetscCall(PetscMalloc5(dim + 1, &pNewCount, dim + 1, &pNewStart, dim + 1, &pNewEnd, dim + 1, &pOldStart, dim + 1, &pOldEnd));
    1786:     offset = 0;
    1787:     for (d = 0; d <= dim; d++) {
    1788:       PetscInt pOldCount, kStart, kEnd, k;

    1790:       pNewStart[d] = offset;
    1791:       PetscCall(DMPlexGetHeightStratum(dm, d, &pOldStart[d], &pOldEnd[d]));
    1792:       PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
    1793:       pOldCount = pOldEnd[d] - pOldStart[d];
    1794:       /* adding the new points */
    1795:       pNewCount[d] = pOldCount + kEnd - kStart;
    1796:       if (!d) {
    1797:         /* removing the cell */
    1798:         pNewCount[d]--;
    1799:       }
    1800:       for (k = kStart; k < kEnd; k++) {
    1801:         PetscInt parent;
    1802:         PetscCall(DMPlexGetTreeParent(K, k, &parent, NULL));
    1803:         if (parent == k) {
    1804:           /* avoid double counting points that won't actually be new */
    1805:           pNewCount[d]--;
    1806:         }
    1807:       }
    1808:       pNewEnd[d] = pNewStart[d] + pNewCount[d];
    1809:       offset     = pNewEnd[d];
    1810:     }
    1811:     PetscCheck(cell >= pOldStart[0] && cell < pOldEnd[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " not in cell range [%" PetscInt_FMT ", %" PetscInt_FMT ")", cell, pOldStart[0], pOldEnd[0]);
    1812:     /* get the current closure of the cell that we are removing */
    1813:     PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));

    1815:     PetscCall(PetscMalloc1(pNewEnd[dim], &newConeSizes));
    1816:     {
    1817:       DMPolytopeType pct, qct;
    1818:       PetscInt       kStart, kEnd, k, closureSizeK, *closureK = NULL, j;

    1820:       PetscCall(DMPlexGetChart(K, &kStart, &kEnd));
    1821:       PetscCall(PetscMalloc4(kEnd - kStart, &Kembedding, kEnd - kStart, &perm, kEnd - kStart, &iperm, kEnd - kStart, &preOrient));

    1823:       for (k = kStart; k < kEnd; k++) {
    1824:         perm[k - kStart]      = k;
    1825:         iperm[k - kStart]     = k - kStart;
    1826:         preOrient[k - kStart] = 0;
    1827:       }

    1829:       PetscCall(DMPlexGetTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
    1830:       for (j = 1; j < closureSizeK; j++) {
    1831:         PetscInt parentOrientA = closureK[2 * j + 1];
    1832:         PetscInt parentOrientB = cellClosure[2 * j + 1];
    1833:         PetscInt p, q;

    1835:         p = closureK[2 * j];
    1836:         q = cellClosure[2 * j];
    1837:         PetscCall(DMPlexGetCellType(K, p, &pct));
    1838:         PetscCall(DMPlexGetCellType(dm, q, &qct));
    1839:         for (d = 0; d <= dim; d++) {
    1840:           if (q >= pOldStart[d] && q < pOldEnd[d]) Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
    1841:         }
    1842:         parentOrientA = DMPolytopeConvertNewOrientation_Internal(pct, parentOrientA);
    1843:         parentOrientB = DMPolytopeConvertNewOrientation_Internal(qct, parentOrientB);
    1844:         if (parentOrientA != parentOrientB) {
    1845:           PetscInt        numChildren, i;
    1846:           const PetscInt *children;

    1848:           PetscCall(DMPlexGetTreeChildren(K, p, &numChildren, &children));
    1849:           for (i = 0; i < numChildren; i++) {
    1850:             PetscInt kPerm, oPerm;

    1852:             k = children[i];
    1853:             PetscCall(DMPlexReferenceTreeGetChildSymmetry(K, p, parentOrientA, 0, k, parentOrientB, &oPerm, &kPerm));
    1854:             /* perm = what refTree position I'm in */
    1855:             perm[kPerm - kStart] = k;
    1856:             /* iperm = who is at this position */
    1857:             iperm[k - kStart]         = kPerm - kStart;
    1858:             preOrient[kPerm - kStart] = oPerm;
    1859:           }
    1860:         }
    1861:       }
    1862:       PetscCall(DMPlexRestoreTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
    1863:     }
    1864:     PetscCall(PetscSectionSetChart(parentSection, 0, pNewEnd[dim]));
    1865:     offset      = 0;
    1866:     numNewCones = 0;
    1867:     for (d = 0; d <= dim; d++) {
    1868:       PetscInt kStart, kEnd, k;
    1869:       PetscInt p;
    1870:       PetscInt size;

    1872:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
    1873:         /* skip cell 0 */
    1874:         if (p == cell) continue;
    1875:         /* old cones to new cones */
    1876:         PetscCall(DMPlexGetConeSize(dm, p, &size));
    1877:         newConeSizes[offset++] = size;
    1878:         numNewCones += size;
    1879:       }

    1881:       PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
    1882:       for (k = kStart; k < kEnd; k++) {
    1883:         PetscInt kParent;

    1885:         PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
    1886:         if (kParent != k) {
    1887:           Kembedding[k] = offset;
    1888:           PetscCall(DMPlexGetConeSize(K, k, &size));
    1889:           newConeSizes[offset++] = size;
    1890:           numNewCones += size;
    1891:           if (kParent != 0) PetscCall(PetscSectionSetDof(parentSection, Kembedding[k], 1));
    1892:         }
    1893:       }
    1894:     }

    1896:     PetscCall(PetscSectionSetUp(parentSection));
    1897:     PetscCall(PetscSectionGetStorageSize(parentSection, &numPointsWithParents));
    1898:     PetscCall(PetscMalloc2(numNewCones, &newCones, numNewCones, &newOrientations));
    1899:     PetscCall(PetscMalloc2(numPointsWithParents, &parents, numPointsWithParents, &childIDs));

    1901:     /* fill new cones */
    1902:     offset = 0;
    1903:     for (d = 0; d <= dim; d++) {
    1904:       PetscInt        kStart, kEnd, k, l;
    1905:       PetscInt        p;
    1906:       PetscInt        size;
    1907:       const PetscInt *cone, *orientation;

    1909:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
    1910:         /* skip cell 0 */
    1911:         if (p == cell) continue;
    1912:         /* old cones to new cones */
    1913:         PetscCall(DMPlexGetConeSize(dm, p, &size));
    1914:         PetscCall(DMPlexGetCone(dm, p, &cone));
    1915:         PetscCall(DMPlexGetConeOrientation(dm, p, &orientation));
    1916:         for (l = 0; l < size; l++) {
    1917:           newCones[offset]          = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
    1918:           newOrientations[offset++] = orientation[l];
    1919:         }
    1920:       }

    1922:       PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
    1923:       for (k = kStart; k < kEnd; k++) {
    1924:         PetscInt kPerm = perm[k], kParent;
    1925:         PetscInt preO  = preOrient[k];

    1927:         PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
    1928:         if (kParent != k) {
    1929:           /* embed new cones */
    1930:           PetscCall(DMPlexGetConeSize(K, k, &size));
    1931:           PetscCall(DMPlexGetCone(K, kPerm, &cone));
    1932:           PetscCall(DMPlexGetConeOrientation(K, kPerm, &orientation));
    1933:           for (l = 0; l < size; l++) {
    1934:             PetscInt       q, m = (preO >= 0) ? ((preO + l) % size) : ((size - (preO + 1) - l) % size);
    1935:             PetscInt       newO, lSize, oTrue;
    1936:             DMPolytopeType ct = DM_NUM_POLYTOPES;

    1938:             q                = iperm[cone[m]];
    1939:             newCones[offset] = Kembedding[q];
    1940:             PetscCall(DMPlexGetConeSize(K, q, &lSize));
    1941:             if (lSize == 2) ct = DM_POLYTOPE_SEGMENT;
    1942:             else if (lSize == 4) ct = DM_POLYTOPE_QUADRILATERAL;
    1943:             oTrue                     = DMPolytopeConvertNewOrientation_Internal(ct, orientation[m]);
    1944:             oTrue                     = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
    1945:             newO                      = DihedralCompose(lSize, oTrue, preOrient[q]);
    1946:             newOrientations[offset++] = DMPolytopeConvertOldOrientation_Internal(ct, newO);
    1947:           }
    1948:           if (kParent != 0) {
    1949:             PetscInt newPoint = Kembedding[kParent];
    1950:             PetscCall(PetscSectionGetOffset(parentSection, Kembedding[k], &pOffset));
    1951:             parents[pOffset]  = newPoint;
    1952:             childIDs[pOffset] = k;
    1953:           }
    1954:         }
    1955:       }
    1956:     }

    1958:     PetscCall(PetscMalloc1(dim * (pNewEnd[dim] - pNewStart[dim]), &newVertexCoords));

    1960:     /* fill coordinates */
    1961:     offset = 0;
    1962:     {
    1963:       PetscInt     kStart, kEnd, l;
    1964:       PetscSection vSection;
    1965:       PetscInt     v;
    1966:       Vec          coords;
    1967:       PetscScalar *coordvals;
    1968:       PetscInt     dof, off;
    1969:       PetscReal    v0[3], J[9], detJ;

    1971:       if (PetscDefined(USE_DEBUG)) {
    1972:         PetscInt k;
    1973:         PetscCall(DMPlexGetHeightStratum(K, 0, &kStart, &kEnd));
    1974:         for (k = kStart; k < kEnd; k++) {
    1975:           PetscCall(DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ));
    1976:           PetscCheck(detJ > 0., PETSC_COMM_SELF, PETSC_ERR_PLIB, "reference tree cell %" PetscInt_FMT " has bad determinant", k);
    1977:         }
    1978:       }
    1979:       PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ));
    1980:       PetscCall(DMGetCoordinateSection(dm, &vSection));
    1981:       PetscCall(DMGetCoordinatesLocal(dm, &coords));
    1982:       PetscCall(VecGetArray(coords, &coordvals));
    1983:       for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
    1984:         PetscCall(PetscSectionGetDof(vSection, v, &dof));
    1985:         PetscCall(PetscSectionGetOffset(vSection, v, &off));
    1986:         for (l = 0; l < dof; l++) newVertexCoords[offset++] = coordvals[off + l];
    1987:       }
    1988:       PetscCall(VecRestoreArray(coords, &coordvals));

    1990:       PetscCall(DMGetCoordinateSection(K, &vSection));
    1991:       PetscCall(DMGetCoordinatesLocal(K, &coords));
    1992:       PetscCall(VecGetArray(coords, &coordvals));
    1993:       PetscCall(DMPlexGetDepthStratum(K, 0, &kStart, &kEnd));
    1994:       for (v = kStart; v < kEnd; v++) {
    1995:         PetscReal       coord[3], newCoord[3];
    1996:         PetscInt        vPerm = perm[v];
    1997:         PetscInt        kParent;
    1998:         const PetscReal xi0[3] = {-1., -1., -1.};

    2000:         PetscCall(DMPlexGetTreeParent(K, v, &kParent, NULL));
    2001:         if (kParent != v) {
    2002:           /* this is a new vertex */
    2003:           PetscCall(PetscSectionGetOffset(vSection, vPerm, &off));
    2004:           for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off + l]);
    2005:           CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
    2006:           for (l = 0; l < dim; ++l) newVertexCoords[offset + l] = newCoord[l];
    2007:           offset += dim;
    2008:         }
    2009:       }
    2010:       PetscCall(VecRestoreArray(coords, &coordvals));
    2011:     }

    2013:     /* need to reverse the order of pNewCount: vertices first, cells last */
    2014:     for (d = 0; d < (dim + 1) / 2; d++) {
    2015:       PetscInt tmp;

    2017:       tmp                = pNewCount[d];
    2018:       pNewCount[d]       = pNewCount[dim - d];
    2019:       pNewCount[dim - d] = tmp;
    2020:     }

    2022:     PetscCall(DMPlexCreateFromDAG(*ncdm, dim, pNewCount, newConeSizes, newCones, newOrientations, newVertexCoords));
    2023:     PetscCall(DMPlexSetReferenceTree(*ncdm, K));
    2024:     PetscCall(DMPlexSetTree(*ncdm, parentSection, parents, childIDs));

    2026:     /* clean up */
    2027:     PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));
    2028:     PetscCall(PetscFree5(pNewCount, pNewStart, pNewEnd, pOldStart, pOldEnd));
    2029:     PetscCall(PetscFree(newConeSizes));
    2030:     PetscCall(PetscFree2(newCones, newOrientations));
    2031:     PetscCall(PetscFree(newVertexCoords));
    2032:     PetscCall(PetscFree2(parents, childIDs));
    2033:     PetscCall(PetscFree4(Kembedding, perm, iperm, preOrient));
    2034:   } else {
    2035:     PetscInt     p, counts[4];
    2036:     PetscInt    *coneSizes, *cones, *orientations;
    2037:     Vec          coordVec;
    2038:     PetscScalar *coords;

    2040:     for (d = 0; d <= dim; d++) {
    2041:       PetscInt dStart, dEnd;

    2043:       PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
    2044:       counts[d] = dEnd - dStart;
    2045:     }
    2046:     PetscCall(PetscMalloc1(pEnd - pStart, &coneSizes));
    2047:     for (p = pStart; p < pEnd; p++) PetscCall(DMPlexGetConeSize(dm, p, &coneSizes[p - pStart]));
    2048:     PetscCall(DMPlexGetCones(dm, &cones));
    2049:     PetscCall(DMPlexGetConeOrientations(dm, &orientations));
    2050:     PetscCall(DMGetCoordinatesLocal(dm, &coordVec));
    2051:     PetscCall(VecGetArray(coordVec, &coords));

    2053:     PetscCall(PetscSectionSetChart(parentSection, pStart, pEnd));
    2054:     PetscCall(PetscSectionSetUp(parentSection));
    2055:     PetscCall(DMPlexCreateFromDAG(*ncdm, dim, counts, coneSizes, cones, orientations, NULL));
    2056:     PetscCall(DMPlexSetReferenceTree(*ncdm, K));
    2057:     PetscCall(DMPlexSetTree(*ncdm, parentSection, NULL, NULL));
    2058:     PetscCall(VecRestoreArray(coordVec, &coords));
    2059:   }
    2060:   PetscCall(PetscSectionDestroy(&parentSection));
    2061:   PetscFunctionReturn(PETSC_SUCCESS);
    2062: }

    2064: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
    2065: {
    2066:   PetscSF              coarseToFineEmbedded;
    2067:   PetscSection         globalCoarse, globalFine;
    2068:   PetscSection         localCoarse, localFine;
    2069:   PetscSection         aSec, cSec;
    2070:   PetscSection         rootIndicesSec, rootMatricesSec;
    2071:   PetscSection         leafIndicesSec, leafMatricesSec;
    2072:   PetscInt            *rootIndices, *leafIndices;
    2073:   PetscScalar         *rootMatrices, *leafMatrices;
    2074:   IS                   aIS;
    2075:   const PetscInt      *anchors;
    2076:   Mat                  cMat;
    2077:   PetscInt             numFields, maxFields;
    2078:   PetscInt             pStartC, pEndC, pStartF, pEndF, p;
    2079:   PetscInt             aStart, aEnd, cStart, cEnd;
    2080:   PetscInt            *maxChildIds;
    2081:   PetscInt            *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
    2082:   const PetscInt    ***perms;
    2083:   const PetscScalar ***flips;

    2085:   PetscFunctionBegin;
    2086:   PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
    2087:   PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
    2088:   PetscCall(DMGetGlobalSection(fine, &globalFine));
    2089:   { /* winnow fine points that don't have global dofs out of the sf */
    2090:     PetscInt        dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
    2091:     const PetscInt *leaves;

    2093:     PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
    2094:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
    2095:       p = leaves ? leaves[l] : l;
    2096:       PetscCall(PetscSectionGetDof(globalFine, p, &dof));
    2097:       PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
    2098:       if ((dof - cdof) > 0) numPointsWithDofs++;
    2099:     }
    2100:     PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
    2101:     for (l = 0, offset = 0; l < nleaves; l++) {
    2102:       p = leaves ? leaves[l] : l;
    2103:       PetscCall(PetscSectionGetDof(globalFine, p, &dof));
    2104:       PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
    2105:       if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
    2106:     }
    2107:     PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
    2108:     PetscCall(PetscFree(pointsWithDofs));
    2109:   }
    2110:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
    2111:   PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
    2112:   for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
    2113:   PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));
    2114:   PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));

    2116:   PetscCall(DMGetLocalSection(coarse, &localCoarse));
    2117:   PetscCall(DMGetGlobalSection(coarse, &globalCoarse));

    2119:   PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
    2120:   PetscCall(ISGetIndices(aIS, &anchors));
    2121:   PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));

    2123:   PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
    2124:   PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));

    2126:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
    2127:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
    2128:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootMatricesSec));
    2129:   PetscCall(PetscSectionSetChart(rootIndicesSec, pStartC, pEndC));
    2130:   PetscCall(PetscSectionSetChart(rootMatricesSec, pStartC, pEndC));
    2131:   PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
    2132:   maxFields = PetscMax(1, numFields);
    2133:   PetscCall(PetscMalloc7(maxFields + 1, &offsets, maxFields + 1, &offsetsCopy, maxFields + 1, &newOffsets, maxFields + 1, &newOffsetsCopy, maxFields + 1, &rowOffsets, maxFields + 1, &numD, maxFields + 1, &numO));
    2134:   PetscCall(PetscMalloc2(maxFields + 1, (PetscInt ****)&perms, maxFields + 1, (PetscScalar ****)&flips));
    2135:   PetscCall(PetscMemzero((void *)perms, (maxFields + 1) * sizeof(const PetscInt **)));
    2136:   PetscCall(PetscMemzero((void *)flips, (maxFields + 1) * sizeof(const PetscScalar **)));

    2138:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
    2139:     PetscInt dof, matSize = 0;
    2140:     PetscInt aDof          = 0;
    2141:     PetscInt cDof          = 0;
    2142:     PetscInt maxChildId    = maxChildIds[p - pStartC];
    2143:     PetscInt numRowIndices = 0;
    2144:     PetscInt numColIndices = 0;
    2145:     PetscInt f;

    2147:     PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
    2148:     if (dof < 0) dof = -(dof + 1);
    2149:     if (p >= aStart && p < aEnd) PetscCall(PetscSectionGetDof(aSec, p, &aDof));
    2150:     if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(cSec, p, &cDof));
    2151:     for (f = 0; f <= numFields; f++) offsets[f] = 0;
    2152:     for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
    2153:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
    2154:       PetscInt *closure = NULL, closureSize, cl;

    2156:       PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
    2157:       for (cl = 0; cl < closureSize; cl++) { /* get the closure */
    2158:         PetscInt c = closure[2 * cl], clDof;

    2160:         PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
    2161:         numRowIndices += clDof;
    2162:         for (f = 0; f < numFields; f++) {
    2163:           PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &clDof));
    2164:           offsets[f + 1] += clDof;
    2165:         }
    2166:       }
    2167:       for (f = 0; f < numFields; f++) {
    2168:         offsets[f + 1] += offsets[f];
    2169:         newOffsets[f + 1] = offsets[f + 1];
    2170:       }
    2171:       /* get the number of indices needed and their field offsets */
    2172:       PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, NULL, NULL, NULL, &numColIndices, NULL, NULL, newOffsets, PETSC_FALSE));
    2173:       PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
    2174:       if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
    2175:         numColIndices = numRowIndices;
    2176:         matSize       = 0;
    2177:       } else if (numFields) { /* we send one submat for each field: sum their sizes */
    2178:         matSize = 0;
    2179:         for (f = 0; f < numFields; f++) {
    2180:           PetscInt numRow, numCol;

    2182:           numRow = offsets[f + 1] - offsets[f];
    2183:           numCol = newOffsets[f + 1] - newOffsets[f];
    2184:           matSize += numRow * numCol;
    2185:         }
    2186:       } else {
    2187:         matSize = numRowIndices * numColIndices;
    2188:       }
    2189:     } else if (maxChildId == -1) {
    2190:       if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
    2191:         PetscInt aOff, a;

    2193:         PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
    2194:         for (f = 0; f < numFields; f++) {
    2195:           PetscInt fDof;

    2197:           PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
    2198:           offsets[f + 1] = fDof;
    2199:         }
    2200:         for (a = 0; a < aDof; a++) {
    2201:           PetscInt anchor = anchors[a + aOff], aLocalDof;

    2203:           PetscCall(PetscSectionGetDof(localCoarse, anchor, &aLocalDof));
    2204:           numColIndices += aLocalDof;
    2205:           for (f = 0; f < numFields; f++) {
    2206:             PetscInt fDof;

    2208:             PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
    2209:             newOffsets[f + 1] += fDof;
    2210:           }
    2211:         }
    2212:         if (numFields) {
    2213:           matSize = 0;
    2214:           for (f = 0; f < numFields; f++) matSize += offsets[f + 1] * newOffsets[f + 1];
    2215:         } else {
    2216:           matSize = numColIndices * dof;
    2217:         }
    2218:       } else { /* no children, and no constraints on dofs: just get the global indices */
    2219:         numColIndices = dof;
    2220:         matSize       = 0;
    2221:       }
    2222:     }
    2223:     /* we will pack the column indices with the field offsets */
    2224:     PetscCall(PetscSectionSetDof(rootIndicesSec, p, numColIndices ? numColIndices + 2 * numFields : 0));
    2225:     PetscCall(PetscSectionSetDof(rootMatricesSec, p, matSize));
    2226:   }
    2227:   PetscCall(PetscSectionSetUp(rootIndicesSec));
    2228:   PetscCall(PetscSectionSetUp(rootMatricesSec));
    2229:   {
    2230:     PetscInt numRootIndices, numRootMatrices;

    2232:     PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
    2233:     PetscCall(PetscSectionGetStorageSize(rootMatricesSec, &numRootMatrices));
    2234:     PetscCall(PetscMalloc2(numRootIndices, &rootIndices, numRootMatrices, &rootMatrices));
    2235:     for (p = pStartC; p < pEndC; p++) {
    2236:       PetscInt     numRowIndices = 0, numColIndices, matSize, dof;
    2237:       PetscInt     pIndOff, pMatOff, f;
    2238:       PetscInt    *pInd;
    2239:       PetscInt     maxChildId = maxChildIds[p - pStartC];
    2240:       PetscScalar *pMat       = NULL;

    2242:       PetscCall(PetscSectionGetDof(rootIndicesSec, p, &numColIndices));
    2243:       if (!numColIndices) continue;
    2244:       for (f = 0; f <= numFields; f++) {
    2245:         offsets[f]        = 0;
    2246:         newOffsets[f]     = 0;
    2247:         offsetsCopy[f]    = 0;
    2248:         newOffsetsCopy[f] = 0;
    2249:       }
    2250:       numColIndices -= 2 * numFields;
    2251:       PetscCall(PetscSectionGetOffset(rootIndicesSec, p, &pIndOff));
    2252:       pInd = &rootIndices[pIndOff];
    2253:       PetscCall(PetscSectionGetDof(rootMatricesSec, p, &matSize));
    2254:       if (matSize) {
    2255:         PetscCall(PetscSectionGetOffset(rootMatricesSec, p, &pMatOff));
    2256:         pMat = &rootMatrices[pMatOff];
    2257:       }
    2258:       PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
    2259:       if (dof < 0) dof = -(dof + 1);
    2260:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
    2261:         PetscInt i, j;

    2263:         if (matSize == 0) { /* don't need to calculate the mat, just the indices */
    2264:           PetscInt numIndices, *indices;
    2265:           PetscCall(DMPlexGetClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
    2266:           PetscCheck(numIndices == numColIndices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "mismatching constraint indices calculations");
    2267:           for (i = 0; i < numColIndices; i++) pInd[i] = indices[i];
    2268:           for (i = 0; i < numFields; i++) {
    2269:             pInd[numColIndices + i]             = offsets[i + 1];
    2270:             pInd[numColIndices + numFields + i] = offsets[i + 1];
    2271:           }
    2272:           PetscCall(DMPlexRestoreClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
    2273:         } else {
    2274:           PetscInt     closureSize, *closure = NULL, cl;
    2275:           PetscScalar *pMatIn, *pMatModified;
    2276:           PetscInt     numPoints, *points;

    2278:           {
    2279:             PetscInt *closure = NULL, closureSize, cl;

    2281:             PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
    2282:             for (cl = 0; cl < closureSize; cl++) { /* get the closure */
    2283:               PetscInt c = closure[2 * cl], clDof;

    2285:               PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
    2286:               numRowIndices += clDof;
    2287:             }
    2288:             PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
    2289:           }

    2291:           PetscCall(DMGetWorkArray(coarse, numRowIndices * numRowIndices, MPIU_SCALAR, &pMatIn));
    2292:           for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
    2293:             for (j = 0; j < numRowIndices; j++) pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
    2294:           }
    2295:           PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
    2296:           for (f = 0; f < maxFields; f++) {
    2297:             if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
    2298:             else PetscCall(PetscSectionGetPointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
    2299:           }
    2300:           if (numFields) {
    2301:             for (cl = 0; cl < closureSize; cl++) {
    2302:               PetscInt c = closure[2 * cl];

    2304:               for (f = 0; f < numFields; f++) {
    2305:                 PetscInt fDof;

    2307:                 PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &fDof));
    2308:                 offsets[f + 1] += fDof;
    2309:               }
    2310:             }
    2311:             for (f = 0; f < numFields; f++) {
    2312:               offsets[f + 1] += offsets[f];
    2313:               newOffsets[f + 1] = offsets[f + 1];
    2314:             }
    2315:           }
    2316:           /* TODO : flips here ? */
    2317:           /* apply hanging node constraints on the right, get the new points and the new offsets */
    2318:           PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, perms, pMatIn, &numPoints, NULL, &points, &pMatModified, newOffsets, PETSC_FALSE));
    2319:           for (f = 0; f < maxFields; f++) {
    2320:             if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
    2321:             else PetscCall(PetscSectionRestorePointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
    2322:           }
    2323:           for (f = 0; f < maxFields; f++) {
    2324:             if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
    2325:             else PetscCall(PetscSectionGetPointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
    2326:           }
    2327:           if (!numFields) {
    2328:             for (i = 0; i < numRowIndices * numColIndices; i++) pMat[i] = pMatModified[i];
    2329:           } else {
    2330:             PetscInt i, j, count;
    2331:             for (f = 0, count = 0; f < numFields; f++) {
    2332:               for (i = offsets[f]; i < offsets[f + 1]; i++) {
    2333:                 for (j = newOffsets[f]; j < newOffsets[f + 1]; j++, count++) pMat[count] = pMatModified[i * numColIndices + j];
    2334:               }
    2335:             }
    2336:           }
    2337:           PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatModified));
    2338:           PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
    2339:           PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatIn));
    2340:           if (numFields) {
    2341:             for (f = 0; f < numFields; f++) {
    2342:               pInd[numColIndices + f]             = offsets[f + 1];
    2343:               pInd[numColIndices + numFields + f] = newOffsets[f + 1];
    2344:             }
    2345:             for (cl = 0; cl < numPoints; cl++) {
    2346:               PetscInt globalOff, c = points[2 * cl];
    2347:               PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
    2348:               PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd));
    2349:             }
    2350:           } else {
    2351:             for (cl = 0; cl < numPoints; cl++) {
    2352:               PetscInt        c    = points[2 * cl], globalOff;
    2353:               const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;

    2355:               PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
    2356:               PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd));
    2357:             }
    2358:           }
    2359:           for (f = 0; f < maxFields; f++) {
    2360:             if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
    2361:             else PetscCall(PetscSectionRestorePointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
    2362:           }
    2363:           PetscCall(DMRestoreWorkArray(coarse, numPoints, MPIU_SCALAR, &points));
    2364:         }
    2365:       } else if (matSize) {
    2366:         PetscInt  cOff;
    2367:         PetscInt *rowIndices, *colIndices, a, aDof = 0, aOff;

    2369:         numRowIndices = dof;
    2370:         PetscCall(DMGetWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
    2371:         PetscCall(DMGetWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
    2372:         PetscCall(PetscSectionGetOffset(cSec, p, &cOff));
    2373:         PetscCall(PetscSectionGetDof(aSec, p, &aDof));
    2374:         PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
    2375:         if (numFields) {
    2376:           for (f = 0; f < numFields; f++) {
    2377:             PetscInt fDof;

    2379:             PetscCall(PetscSectionGetFieldDof(cSec, p, f, &fDof));
    2380:             offsets[f + 1] = fDof;
    2381:             for (a = 0; a < aDof; a++) {
    2382:               PetscInt anchor = anchors[a + aOff];
    2383:               PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
    2384:               newOffsets[f + 1] += fDof;
    2385:             }
    2386:           }
    2387:           for (f = 0; f < numFields; f++) {
    2388:             offsets[f + 1] += offsets[f];
    2389:             offsetsCopy[f + 1] = offsets[f + 1];
    2390:             newOffsets[f + 1] += newOffsets[f];
    2391:             newOffsetsCopy[f + 1] = newOffsets[f + 1];
    2392:           }
    2393:           PetscCall(DMPlexGetIndicesPointFields_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, -1, NULL, rowIndices));
    2394:           for (a = 0; a < aDof; a++) {
    2395:             PetscInt anchor = anchors[a + aOff], lOff;
    2396:             PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
    2397:             PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, -1, NULL, colIndices));
    2398:           }
    2399:         } else {
    2400:           PetscCall(DMPlexGetIndicesPoint_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, NULL, rowIndices));
    2401:           for (a = 0; a < aDof; a++) {
    2402:             PetscInt anchor = anchors[a + aOff], lOff;
    2403:             PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
    2404:             PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, NULL, colIndices));
    2405:           }
    2406:         }
    2407:         if (numFields) {
    2408:           PetscInt count, a;

    2410:           for (f = 0, count = 0; f < numFields; f++) {
    2411:             PetscInt iSize = offsets[f + 1] - offsets[f];
    2412:             PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
    2413:             PetscCall(MatGetValues(cMat, iSize, &rowIndices[offsets[f]], jSize, &colIndices[newOffsets[f]], &pMat[count]));
    2414:             count += iSize * jSize;
    2415:             pInd[numColIndices + f]             = offsets[f + 1];
    2416:             pInd[numColIndices + numFields + f] = newOffsets[f + 1];
    2417:           }
    2418:           for (a = 0; a < aDof; a++) {
    2419:             PetscInt anchor = anchors[a + aOff];
    2420:             PetscInt gOff;
    2421:             PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
    2422:             PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, -1, NULL, pInd));
    2423:           }
    2424:         } else {
    2425:           PetscInt a;
    2426:           PetscCall(MatGetValues(cMat, numRowIndices, rowIndices, numColIndices, colIndices, pMat));
    2427:           for (a = 0; a < aDof; a++) {
    2428:             PetscInt anchor = anchors[a + aOff];
    2429:             PetscInt gOff;
    2430:             PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
    2431:             PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, NULL, pInd));
    2432:           }
    2433:         }
    2434:         PetscCall(DMRestoreWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
    2435:         PetscCall(DMRestoreWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
    2436:       } else {
    2437:         PetscInt gOff;

    2439:         PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
    2440:         if (numFields) {
    2441:           for (f = 0; f < numFields; f++) {
    2442:             PetscInt fDof;
    2443:             PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
    2444:             offsets[f + 1] = fDof + offsets[f];
    2445:           }
    2446:           for (f = 0; f < numFields; f++) {
    2447:             pInd[numColIndices + f]             = offsets[f + 1];
    2448:             pInd[numColIndices + numFields + f] = offsets[f + 1];
    2449:           }
    2450:           PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
    2451:         } else {
    2452:           PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
    2453:         }
    2454:       }
    2455:     }
    2456:     PetscCall(PetscFree(maxChildIds));
    2457:   }
    2458:   {
    2459:     PetscSF   indicesSF, matricesSF;
    2460:     PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;

    2462:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
    2463:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafMatricesSec));
    2464:     PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootIndicesSec, &remoteOffsetsIndices, leafIndicesSec));
    2465:     PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootMatricesSec, &remoteOffsetsMatrices, leafMatricesSec));
    2466:     PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootIndicesSec, remoteOffsetsIndices, leafIndicesSec, &indicesSF));
    2467:     PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootMatricesSec, remoteOffsetsMatrices, leafMatricesSec, &matricesSF));
    2468:     PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
    2469:     PetscCall(PetscFree(remoteOffsetsIndices));
    2470:     PetscCall(PetscFree(remoteOffsetsMatrices));
    2471:     PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numLeafIndices));
    2472:     PetscCall(PetscSectionGetStorageSize(leafMatricesSec, &numLeafMatrices));
    2473:     PetscCall(PetscMalloc2(numLeafIndices, &leafIndices, numLeafMatrices, &leafMatrices));
    2474:     PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
    2475:     PetscCall(PetscSFBcastBegin(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
    2476:     PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
    2477:     PetscCall(PetscSFBcastEnd(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
    2478:     PetscCall(PetscSFDestroy(&matricesSF));
    2479:     PetscCall(PetscSFDestroy(&indicesSF));
    2480:     PetscCall(PetscFree2(rootIndices, rootMatrices));
    2481:     PetscCall(PetscSectionDestroy(&rootIndicesSec));
    2482:     PetscCall(PetscSectionDestroy(&rootMatricesSec));
    2483:   }
    2484:   /* count to preallocate */
    2485:   PetscCall(DMGetLocalSection(fine, &localFine));
    2486:   {
    2487:     PetscInt       nGlobal;
    2488:     PetscInt      *dnnz, *onnz;
    2489:     PetscLayout    rowMap, colMap;
    2490:     PetscInt       rowStart, rowEnd, colStart, colEnd;
    2491:     PetscInt       maxDof;
    2492:     PetscInt      *rowIndices;
    2493:     DM             refTree;
    2494:     PetscInt     **refPointFieldN;
    2495:     PetscScalar ***refPointFieldMats;
    2496:     PetscSection   refConSec, refAnSec;
    2497:     PetscInt       pRefStart, pRefEnd, maxConDof, maxColumns, leafStart, leafEnd;
    2498:     PetscScalar   *pointWork;

    2500:     PetscCall(PetscSectionGetConstrainedStorageSize(globalFine, &nGlobal));
    2501:     PetscCall(PetscCalloc2(nGlobal, &dnnz, nGlobal, &onnz));
    2502:     PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
    2503:     PetscCall(PetscLayoutSetUp(rowMap));
    2504:     PetscCall(PetscLayoutSetUp(colMap));
    2505:     PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
    2506:     PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
    2507:     PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
    2508:     PetscCall(PetscSectionGetChart(leafIndicesSec, &leafStart, &leafEnd));
    2509:     PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
    2510:     for (p = leafStart; p < leafEnd; p++) {
    2511:       PetscInt gDof, gcDof, gOff;
    2512:       PetscInt numColIndices, pIndOff, *pInd;
    2513:       PetscInt matSize;
    2514:       PetscInt i;

    2516:       PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
    2517:       PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
    2518:       if ((gDof - gcDof) <= 0) continue;
    2519:       PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
    2520:       PetscCheck(gOff >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I though having global dofs meant a non-negative offset");
    2521:       PetscCheck(gOff >= rowStart && (gOff + gDof - gcDof) <= rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I thought the row map would constrain the global dofs");
    2522:       PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices));
    2523:       PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff));
    2524:       numColIndices -= 2 * numFields;
    2525:       PetscCheck(numColIndices > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "global fine dof with no dofs to interpolate from");
    2526:       pInd              = &leafIndices[pIndOff];
    2527:       offsets[0]        = 0;
    2528:       offsetsCopy[0]    = 0;
    2529:       newOffsets[0]     = 0;
    2530:       newOffsetsCopy[0] = 0;
    2531:       if (numFields) {
    2532:         PetscInt f;
    2533:         for (f = 0; f < numFields; f++) {
    2534:           PetscInt rowDof;

    2536:           PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
    2537:           offsets[f + 1]     = offsets[f] + rowDof;
    2538:           offsetsCopy[f + 1] = offsets[f + 1];
    2539:           newOffsets[f + 1]  = pInd[numColIndices + numFields + f];
    2540:           numD[f]            = 0;
    2541:           numO[f]            = 0;
    2542:         }
    2543:         PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
    2544:         for (f = 0; f < numFields; f++) {
    2545:           PetscInt colOffset    = newOffsets[f];
    2546:           PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];

    2548:           for (i = 0; i < numFieldCols; i++) {
    2549:             PetscInt gInd = pInd[i + colOffset];

    2551:             if (gInd >= colStart && gInd < colEnd) {
    2552:               numD[f]++;
    2553:             } else if (gInd >= 0) { /* negative means non-entry */
    2554:               numO[f]++;
    2555:             }
    2556:           }
    2557:         }
    2558:       } else {
    2559:         PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
    2560:         numD[0] = 0;
    2561:         numO[0] = 0;
    2562:         for (i = 0; i < numColIndices; i++) {
    2563:           PetscInt gInd = pInd[i];

    2565:           if (gInd >= colStart && gInd < colEnd) {
    2566:             numD[0]++;
    2567:           } else if (gInd >= 0) { /* negative means non-entry */
    2568:             numO[0]++;
    2569:           }
    2570:         }
    2571:       }
    2572:       PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize));
    2573:       if (!matSize) { /* incoming matrix is identity */
    2574:         PetscInt childId;

    2576:         childId = childIds[p - pStartF];
    2577:         if (childId < 0) { /* no child interpolation: one nnz per */
    2578:           if (numFields) {
    2579:             PetscInt f;
    2580:             for (f = 0; f < numFields; f++) {
    2581:               PetscInt numRows = offsets[f + 1] - offsets[f], row;
    2582:               for (row = 0; row < numRows; row++) {
    2583:                 PetscInt gIndCoarse = pInd[newOffsets[f] + row];
    2584:                 PetscInt gIndFine   = rowIndices[offsets[f] + row];
    2585:                 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
    2586:                   PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2587:                   dnnz[gIndFine - rowStart] = 1;
    2588:                 } else if (gIndCoarse >= 0) { /* remote */
    2589:                   PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2590:                   onnz[gIndFine - rowStart] = 1;
    2591:                 } else { /* constrained */
    2592:                   PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2593:                 }
    2594:               }
    2595:             }
    2596:           } else {
    2597:             PetscInt i;
    2598:             for (i = 0; i < gDof; i++) {
    2599:               PetscInt gIndCoarse = pInd[i];
    2600:               PetscInt gIndFine   = rowIndices[i];
    2601:               if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
    2602:                 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2603:                 dnnz[gIndFine - rowStart] = 1;
    2604:               } else if (gIndCoarse >= 0) { /* remote */
    2605:                 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2606:                 onnz[gIndFine - rowStart] = 1;
    2607:               } else { /* constrained */
    2608:                 PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2609:               }
    2610:             }
    2611:           }
    2612:         } else { /* interpolate from all */
    2613:           if (numFields) {
    2614:             PetscInt f;
    2615:             for (f = 0; f < numFields; f++) {
    2616:               PetscInt numRows = offsets[f + 1] - offsets[f], row;
    2617:               for (row = 0; row < numRows; row++) {
    2618:                 PetscInt gIndFine = rowIndices[offsets[f] + row];
    2619:                 if (gIndFine >= 0) {
    2620:                   PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2621:                   dnnz[gIndFine - rowStart] = numD[f];
    2622:                   onnz[gIndFine - rowStart] = numO[f];
    2623:                 }
    2624:               }
    2625:             }
    2626:           } else {
    2627:             PetscInt i;
    2628:             for (i = 0; i < gDof; i++) {
    2629:               PetscInt gIndFine = rowIndices[i];
    2630:               if (gIndFine >= 0) {
    2631:                 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2632:                 dnnz[gIndFine - rowStart] = numD[0];
    2633:                 onnz[gIndFine - rowStart] = numO[0];
    2634:               }
    2635:             }
    2636:           }
    2637:         }
    2638:       } else { /* interpolate from all */
    2639:         if (numFields) {
    2640:           PetscInt f;
    2641:           for (f = 0; f < numFields; f++) {
    2642:             PetscInt numRows = offsets[f + 1] - offsets[f], row;
    2643:             for (row = 0; row < numRows; row++) {
    2644:               PetscInt gIndFine = rowIndices[offsets[f] + row];
    2645:               if (gIndFine >= 0) {
    2646:                 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2647:                 dnnz[gIndFine - rowStart] = numD[f];
    2648:                 onnz[gIndFine - rowStart] = numO[f];
    2649:               }
    2650:             }
    2651:           }
    2652:         } else { /* every dof get a full row */
    2653:           PetscInt i;
    2654:           for (i = 0; i < gDof; i++) {
    2655:             PetscInt gIndFine = rowIndices[i];
    2656:             if (gIndFine >= 0) {
    2657:               PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
    2658:               dnnz[gIndFine - rowStart] = numD[0];
    2659:               onnz[gIndFine - rowStart] = numO[0];
    2660:             }
    2661:           }
    2662:         }
    2663:       }
    2664:     }
    2665:     PetscCall(MatXAIJSetPreallocation(mat, 1, dnnz, onnz, NULL, NULL));
    2666:     PetscCall(PetscFree2(dnnz, onnz));

    2668:     PetscCall(DMPlexGetReferenceTree(fine, &refTree));
    2669:     PetscCall(DMCopyDisc(fine, refTree));
    2670:     PetscCall(DMSetLocalSection(refTree, NULL));
    2671:     PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
    2672:     PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
    2673:     PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
    2674:     PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL));
    2675:     PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
    2676:     PetscCall(PetscSectionGetMaxDof(refConSec, &maxConDof));
    2677:     PetscCall(PetscSectionGetMaxDof(leafIndicesSec, &maxColumns));
    2678:     PetscCall(PetscMalloc1(maxConDof * maxColumns, &pointWork));
    2679:     for (p = leafStart; p < leafEnd; p++) {
    2680:       PetscInt gDof, gcDof, gOff;
    2681:       PetscInt numColIndices, pIndOff, *pInd;
    2682:       PetscInt matSize;
    2683:       PetscInt childId;

    2685:       PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
    2686:       PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
    2687:       if ((gDof - gcDof) <= 0) continue;
    2688:       childId = childIds[p - pStartF];
    2689:       PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
    2690:       PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices));
    2691:       PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff));
    2692:       numColIndices -= 2 * numFields;
    2693:       pInd              = &leafIndices[pIndOff];
    2694:       offsets[0]        = 0;
    2695:       offsetsCopy[0]    = 0;
    2696:       newOffsets[0]     = 0;
    2697:       newOffsetsCopy[0] = 0;
    2698:       rowOffsets[0]     = 0;
    2699:       if (numFields) {
    2700:         PetscInt f;
    2701:         for (f = 0; f < numFields; f++) {
    2702:           PetscInt rowDof;

    2704:           PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
    2705:           offsets[f + 1]     = offsets[f] + rowDof;
    2706:           offsetsCopy[f + 1] = offsets[f + 1];
    2707:           rowOffsets[f + 1]  = pInd[numColIndices + f];
    2708:           newOffsets[f + 1]  = pInd[numColIndices + numFields + f];
    2709:         }
    2710:         PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
    2711:       } else {
    2712:         PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
    2713:       }
    2714:       PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize));
    2715:       if (!matSize) {      /* incoming matrix is identity */
    2716:         if (childId < 0) { /* no child interpolation: scatter */
    2717:           if (numFields) {
    2718:             PetscInt f;
    2719:             for (f = 0; f < numFields; f++) {
    2720:               PetscInt numRows = offsets[f + 1] - offsets[f], row;
    2721:               for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[offsets[f] + row], pInd[newOffsets[f] + row], 1., INSERT_VALUES));
    2722:             }
    2723:           } else {
    2724:             PetscInt numRows = gDof, row;
    2725:             for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[row], pInd[row], 1., INSERT_VALUES));
    2726:           }
    2727:         } else { /* interpolate from all */
    2728:           if (numFields) {
    2729:             PetscInt f;
    2730:             for (f = 0; f < numFields; f++) {
    2731:               PetscInt numRows = offsets[f + 1] - offsets[f];
    2732:               PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
    2733:               PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], refPointFieldMats[childId - pRefStart][f], INSERT_VALUES));
    2734:             }
    2735:           } else {
    2736:             PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, refPointFieldMats[childId - pRefStart][0], INSERT_VALUES));
    2737:           }
    2738:         }
    2739:       } else { /* interpolate from all */
    2740:         PetscInt     pMatOff;
    2741:         PetscScalar *pMat;

    2743:         PetscCall(PetscSectionGetOffset(leafMatricesSec, p, &pMatOff));
    2744:         pMat = &leafMatrices[pMatOff];
    2745:         if (childId < 0) { /* copy the incoming matrix */
    2746:           if (numFields) {
    2747:             PetscInt f, count;
    2748:             for (f = 0, count = 0; f < numFields; f++) {
    2749:               PetscInt     numRows   = offsets[f + 1] - offsets[f];
    2750:               PetscInt     numCols   = newOffsets[f + 1] - newOffsets[f];
    2751:               PetscInt     numInRows = rowOffsets[f + 1] - rowOffsets[f];
    2752:               PetscScalar *inMat     = &pMat[count];

    2754:               PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], inMat, INSERT_VALUES));
    2755:               count += numCols * numInRows;
    2756:             }
    2757:           } else {
    2758:             PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, pMat, INSERT_VALUES));
    2759:           }
    2760:         } else { /* multiply the incoming matrix by the child interpolation */
    2761:           if (numFields) {
    2762:             PetscInt f, count;
    2763:             for (f = 0, count = 0; f < numFields; f++) {
    2764:               PetscInt     numRows   = offsets[f + 1] - offsets[f];
    2765:               PetscInt     numCols   = newOffsets[f + 1] - newOffsets[f];
    2766:               PetscInt     numInRows = rowOffsets[f + 1] - rowOffsets[f];
    2767:               PetscScalar *inMat     = &pMat[count];
    2768:               PetscInt     i, j, k;
    2769:               PetscCheck(refPointFieldN[childId - pRefStart][f] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
    2770:               for (i = 0; i < numRows; i++) {
    2771:                 for (j = 0; j < numCols; j++) {
    2772:                   PetscScalar val = 0.;
    2773:                   for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
    2774:                   pointWork[i * numCols + j] = val;
    2775:                 }
    2776:               }
    2777:               PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], pointWork, INSERT_VALUES));
    2778:               count += numCols * numInRows;
    2779:             }
    2780:           } else { /* every dof gets a full row */
    2781:             PetscInt numRows   = gDof;
    2782:             PetscInt numCols   = numColIndices;
    2783:             PetscInt numInRows = matSize / numColIndices;
    2784:             PetscInt i, j, k;
    2785:             PetscCheck(refPointFieldN[childId - pRefStart][0] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
    2786:             for (i = 0; i < numRows; i++) {
    2787:               for (j = 0; j < numCols; j++) {
    2788:                 PetscScalar val = 0.;
    2789:                 for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
    2790:                 pointWork[i * numCols + j] = val;
    2791:               }
    2792:             }
    2793:             PetscCall(MatSetValues(mat, numRows, rowIndices, numCols, pInd, pointWork, INSERT_VALUES));
    2794:           }
    2795:         }
    2796:       }
    2797:     }
    2798:     PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
    2799:     PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
    2800:     PetscCall(PetscFree(pointWork));
    2801:   }
    2802:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
    2803:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
    2804:   PetscCall(PetscSectionDestroy(&leafIndicesSec));
    2805:   PetscCall(PetscSectionDestroy(&leafMatricesSec));
    2806:   PetscCall(PetscFree2(leafIndices, leafMatrices));
    2807:   PetscCall(PetscFree2(*(PetscInt ****)&perms, *(PetscScalar ****)&flips));
    2808:   PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
    2809:   PetscCall(ISRestoreIndices(aIS, &anchors));
    2810:   PetscFunctionReturn(PETSC_SUCCESS);
    2811: }

    2813: /*
    2814:  * Assuming a nodal basis (w.r.t. the dual basis) basis:
    2815:  *
    2816:  * for each coarse dof \phi^c_i:
    2817:  *   for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
    2818:  *     for each fine dof \phi^f_j;
    2819:  *       a_{i,j} = 0;
    2820:  *       for each fine dof \phi^f_k:
    2821:  *         a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
    2822:  *                    [^^^ this is = \phi^c_i ^^^]
    2823:  */
    2824: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
    2825: {
    2826:   PetscDS      ds;
    2827:   PetscSection section, cSection;
    2828:   DMLabel      canonical, depth;
    2829:   Mat          cMat, mat;
    2830:   PetscInt    *nnz;
    2831:   PetscInt     f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
    2832:   PetscInt     m, n;
    2833:   PetscScalar *pointScalar;
    2834:   PetscReal   *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;

    2836:   PetscFunctionBegin;
    2837:   PetscCall(DMGetLocalSection(refTree, &section));
    2838:   PetscCall(DMGetDimension(refTree, &dim));
    2839:   PetscCall(PetscMalloc6(dim, &v0, dim, &v0parent, dim, &vtmp, dim * dim, &J, dim * dim, &Jparent, dim * dim, &invJ));
    2840:   PetscCall(PetscMalloc2(dim, &pointScalar, dim, &pointRef));
    2841:   PetscCall(DMGetDS(refTree, &ds));
    2842:   PetscCall(PetscDSGetNumFields(ds, &numFields));
    2843:   PetscCall(PetscSectionGetNumFields(section, &numSecFields));
    2844:   PetscCall(DMGetLabel(refTree, "canonical", &canonical));
    2845:   PetscCall(DMGetLabel(refTree, "depth", &depth));
    2846:   PetscCall(DMGetDefaultConstraints(refTree, &cSection, &cMat, NULL));
    2847:   PetscCall(DMPlexGetChart(refTree, &pStart, &pEnd));
    2848:   PetscCall(DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd));
    2849:   PetscCall(MatGetSize(cMat, &n, &m)); /* the injector has transpose sizes from the constraint matrix */
    2850:   /* Step 1: compute non-zero pattern.  A proper subset of constraint matrix non-zero */
    2851:   PetscCall(PetscCalloc1(m, &nnz));
    2852:   for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
    2853:     const PetscInt *children;
    2854:     PetscInt        numChildren;
    2855:     PetscInt        i, numChildDof, numSelfDof;

    2857:     if (canonical) {
    2858:       PetscInt pCanonical;
    2859:       PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
    2860:       if (p != pCanonical) continue;
    2861:     }
    2862:     PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
    2863:     if (!numChildren) continue;
    2864:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
    2865:       PetscInt child = children[i];
    2866:       PetscInt dof;

    2868:       PetscCall(PetscSectionGetDof(section, child, &dof));
    2869:       numChildDof += dof;
    2870:     }
    2871:     PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
    2872:     if (!numChildDof || !numSelfDof) continue;
    2873:     for (f = 0; f < numFields; f++) {
    2874:       PetscInt selfOff;

    2876:       if (numSecFields) { /* count the dofs for just this field */
    2877:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
    2878:           PetscInt child = children[i];
    2879:           PetscInt dof;

    2881:           PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
    2882:           numChildDof += dof;
    2883:         }
    2884:         PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
    2885:         PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
    2886:       } else {
    2887:         PetscCall(PetscSectionGetOffset(section, p, &selfOff));
    2888:       }
    2889:       for (i = 0; i < numSelfDof; i++) nnz[selfOff + i] = numChildDof;
    2890:     }
    2891:   }
    2892:   PetscCall(MatCreateAIJ(PETSC_COMM_SELF, m, n, m, n, -1, nnz, -1, NULL, &mat));
    2893:   PetscCall(PetscFree(nnz));
    2894:   /* Setp 2: compute entries */
    2895:   for (p = pStart; p < pEnd; p++) {
    2896:     const PetscInt *children;
    2897:     PetscInt        numChildren;
    2898:     PetscInt        i, numChildDof, numSelfDof;

    2900:     /* same conditions about when entries occur */
    2901:     if (canonical) {
    2902:       PetscInt pCanonical;
    2903:       PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
    2904:       if (p != pCanonical) continue;
    2905:     }
    2906:     PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
    2907:     if (!numChildren) continue;
    2908:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
    2909:       PetscInt child = children[i];
    2910:       PetscInt dof;

    2912:       PetscCall(PetscSectionGetDof(section, child, &dof));
    2913:       numChildDof += dof;
    2914:     }
    2915:     PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
    2916:     if (!numChildDof || !numSelfDof) continue;

    2918:     for (f = 0; f < numFields; f++) {
    2919:       PetscInt        pI = -1, cI = -1;
    2920:       PetscInt        selfOff, Nc, parentCell;
    2921:       PetscInt        cellShapeOff;
    2922:       PetscObject     disc;
    2923:       PetscDualSpace  dsp;
    2924:       PetscClassId    classId;
    2925:       PetscScalar    *pointMat;
    2926:       PetscInt       *matRows, *matCols;
    2927:       PetscInt        pO = PETSC_INT_MIN;
    2928:       const PetscInt *depthNumDof;

    2930:       if (numSecFields) {
    2931:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
    2932:           PetscInt child = children[i];
    2933:           PetscInt dof;

    2935:           PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
    2936:           numChildDof += dof;
    2937:         }
    2938:         PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
    2939:         PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
    2940:       } else {
    2941:         PetscCall(PetscSectionGetOffset(section, p, &selfOff));
    2942:       }

    2944:       /* find a cell whose closure contains p */
    2945:       if (p >= cStart && p < cEnd) {
    2946:         parentCell = p;
    2947:       } else {
    2948:         PetscInt *star = NULL;
    2949:         PetscInt  numStar;

    2951:         parentCell = -1;
    2952:         PetscCall(DMPlexGetTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
    2953:         for (i = numStar - 1; i >= 0; i--) {
    2954:           PetscInt c = star[2 * i];

    2956:           if (c >= cStart && c < cEnd) {
    2957:             parentCell = c;
    2958:             break;
    2959:           }
    2960:         }
    2961:         PetscCall(DMPlexRestoreTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
    2962:       }
    2963:       /* determine the offset of p's shape functions within parentCell's shape functions */
    2964:       PetscCall(PetscDSGetDiscretization(ds, f, &disc));
    2965:       PetscCall(PetscObjectGetClassId(disc, &classId));
    2966:       if (classId == PETSCFE_CLASSID) {
    2967:         PetscCall(PetscFEGetDualSpace((PetscFE)disc, &dsp));
    2968:       } else if (classId == PETSCFV_CLASSID) {
    2969:         PetscCall(PetscFVGetDualSpace((PetscFV)disc, &dsp));
    2970:       } else {
    2971:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported discretization object");
    2972:       }
    2973:       PetscCall(PetscDualSpaceGetNumDof(dsp, &depthNumDof));
    2974:       PetscCall(PetscDualSpaceGetNumComponents(dsp, &Nc));
    2975:       {
    2976:         PetscInt *closure = NULL;
    2977:         PetscInt  numClosure;

    2979:         PetscCall(DMPlexGetTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
    2980:         for (i = 0, pI = -1, cellShapeOff = 0; i < numClosure; i++) {
    2981:           PetscInt point = closure[2 * i], pointDepth;

    2983:           pO = closure[2 * i + 1];
    2984:           if (point == p) {
    2985:             pI = i;
    2986:             break;
    2987:           }
    2988:           PetscCall(DMLabelGetValue(depth, point, &pointDepth));
    2989:           cellShapeOff += depthNumDof[pointDepth];
    2990:         }
    2991:         PetscCall(DMPlexRestoreTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
    2992:       }

    2994:       PetscCall(DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
    2995:       PetscCall(DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
    2996:       matCols = matRows + numSelfDof;
    2997:       for (i = 0; i < numSelfDof; i++) matRows[i] = selfOff + i;
    2998:       for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
    2999:       {
    3000:         PetscInt colOff = 0;

    3002:         for (i = 0; i < numChildren; i++) {
    3003:           PetscInt child = children[i];
    3004:           PetscInt dof, off, j;

    3006:           if (numSecFields) {
    3007:             PetscCall(PetscSectionGetFieldDof(cSection, child, f, &dof));
    3008:             PetscCall(PetscSectionGetFieldOffset(cSection, child, f, &off));
    3009:           } else {
    3010:             PetscCall(PetscSectionGetDof(cSection, child, &dof));
    3011:             PetscCall(PetscSectionGetOffset(cSection, child, &off));
    3012:           }

    3014:           for (j = 0; j < dof; j++) matCols[colOff++] = off + j;
    3015:         }
    3016:       }
    3017:       if (classId == PETSCFE_CLASSID) {
    3018:         PetscFE              fe = (PetscFE)disc;
    3019:         PetscInt             fSize;
    3020:         const PetscInt    ***perms;
    3021:         const PetscScalar ***flips;
    3022:         const PetscInt      *pperms;

    3024:         PetscCall(PetscFEGetDualSpace(fe, &dsp));
    3025:         PetscCall(PetscDualSpaceGetDimension(dsp, &fSize));
    3026:         PetscCall(PetscDualSpaceGetSymmetries(dsp, &perms, &flips));
    3027:         pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL;
    3028:         for (i = 0; i < numSelfDof; i++) { /* for every shape function */
    3029:           PetscQuadrature  q;
    3030:           PetscInt         dim, thisNc, numPoints, j, k;
    3031:           const PetscReal *points;
    3032:           const PetscReal *weights;
    3033:           PetscInt        *closure = NULL;
    3034:           PetscInt         numClosure;
    3035:           PetscInt         iCell              = pperms ? pperms[i] : i;
    3036:           PetscInt         parentCellShapeDof = cellShapeOff + iCell;
    3037:           PetscTabulation  Tparent;

    3039:           PetscCall(PetscDualSpaceGetFunctional(dsp, parentCellShapeDof, &q));
    3040:           PetscCall(PetscQuadratureGetData(q, &dim, &thisNc, &numPoints, &points, &weights));
    3041:           PetscCheck(thisNc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT, thisNc, Nc);
    3042:           PetscCall(PetscFECreateTabulation(fe, 1, numPoints, points, 0, &Tparent)); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
    3043:           for (j = 0; j < numPoints; j++) {
    3044:             PetscInt           childCell = -1;
    3045:             PetscReal         *parentValAtPoint;
    3046:             const PetscReal    xi0[3]    = {-1., -1., -1.};
    3047:             const PetscReal   *pointReal = &points[dim * j];
    3048:             const PetscScalar *point;
    3049:             PetscTabulation    Tchild;
    3050:             PetscInt           childCellShapeOff, pointMatOff;
    3051: #if defined(PETSC_USE_COMPLEX)
    3052:             PetscInt d;

    3054:             for (d = 0; d < dim; d++) pointScalar[d] = points[dim * j + d];
    3055:             point = pointScalar;
    3056: #else
    3057:             point = pointReal;
    3058: #endif

    3060:             parentValAtPoint = &Tparent->T[0][(fSize * j + parentCellShapeDof) * Nc];

    3062:             for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
    3063:               PetscInt  child = children[k];
    3064:               PetscInt *star  = NULL;
    3065:               PetscInt  numStar, s;

    3067:               PetscCall(DMPlexGetTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
    3068:               for (s = numStar - 1; s >= 0; s--) {
    3069:                 PetscInt c = star[2 * s];

    3071:                 if (c < cStart || c >= cEnd) continue;
    3072:                 PetscCall(DMPlexLocatePoint_Internal(refTree, dim, point, c, &childCell));
    3073:                 if (childCell >= 0) break;
    3074:               }
    3075:               PetscCall(DMPlexRestoreTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
    3076:               if (childCell >= 0) break;
    3077:             }
    3078:             PetscCheck(childCell >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not locate quadrature point");
    3079:             PetscCall(DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ));
    3080:             PetscCall(DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent));
    3081:             CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
    3082:             CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);

    3084:             PetscCall(PetscFECreateTabulation(fe, 1, 1, pointRef, 0, &Tchild));
    3085:             PetscCall(DMPlexGetTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
    3086:             for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
    3087:               PetscInt        child = children[k], childDepth, childDof, childO = PETSC_INT_MIN;
    3088:               PetscInt        l;
    3089:               const PetscInt *cperms;

    3091:               PetscCall(DMLabelGetValue(depth, child, &childDepth));
    3092:               childDof = depthNumDof[childDepth];
    3093:               for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) {
    3094:                 PetscInt point = closure[2 * l];
    3095:                 PetscInt pointDepth;

    3097:                 childO = closure[2 * l + 1];
    3098:                 if (point == child) {
    3099:                   cI = l;
    3100:                   break;
    3101:                 }
    3102:                 PetscCall(DMLabelGetValue(depth, point, &pointDepth));
    3103:                 childCellShapeOff += depthNumDof[pointDepth];
    3104:               }
    3105:               if (l == numClosure) {
    3106:                 pointMatOff += childDof;
    3107:                 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
    3108:               }
    3109:               cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL;
    3110:               for (l = 0; l < childDof; l++) {
    3111:                 PetscInt   lCell        = cperms ? cperms[l] : l;
    3112:                 PetscInt   childCellDof = childCellShapeOff + lCell;
    3113:                 PetscReal *childValAtPoint;
    3114:                 PetscReal  val = 0.;

    3116:                 childValAtPoint = &Tchild->T[0][childCellDof * Nc];
    3117:                 for (m = 0; m < Nc; m++) val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];

    3119:                 pointMat[i * numChildDof + pointMatOff + l] += val;
    3120:               }
    3121:               pointMatOff += childDof;
    3122:             }
    3123:             PetscCall(DMPlexRestoreTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
    3124:             PetscCall(PetscTabulationDestroy(&Tchild));
    3125:           }
    3126:           PetscCall(PetscTabulationDestroy(&Tparent));
    3127:         }
    3128:       } else { /* just the volume-weighted averages of the children */
    3129:         PetscReal parentVol;
    3130:         PetscInt  childCell;

    3132:         PetscCall(DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL));
    3133:         for (i = 0, childCell = 0; i < numChildren; i++) {
    3134:           PetscInt  child = children[i], j;
    3135:           PetscReal childVol;

    3137:           if (child < cStart || child >= cEnd) continue;
    3138:           PetscCall(DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL));
    3139:           for (j = 0; j < Nc; j++) pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
    3140:           childCell++;
    3141:         }
    3142:       }
    3143:       /* Insert pointMat into mat */
    3144:       PetscCall(MatSetValues(mat, numSelfDof, matRows, numChildDof, matCols, pointMat, INSERT_VALUES));
    3145:       PetscCall(DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
    3146:       PetscCall(DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
    3147:     }
    3148:   }
    3149:   PetscCall(PetscFree6(v0, v0parent, vtmp, J, Jparent, invJ));
    3150:   PetscCall(PetscFree2(pointScalar, pointRef));
    3151:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
    3152:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
    3153:   *inj = mat;
    3154:   PetscFunctionReturn(PETSC_SUCCESS);
    3155: }

    3157: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
    3158: {
    3159:   PetscDS        ds;
    3160:   PetscInt       numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
    3161:   PetscScalar ***refPointFieldMats;
    3162:   PetscSection   refConSec, refSection;

    3164:   PetscFunctionBegin;
    3165:   PetscCall(DMGetDS(refTree, &ds));
    3166:   PetscCall(PetscDSGetNumFields(ds, &numFields));
    3167:   PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
    3168:   PetscCall(DMGetLocalSection(refTree, &refSection));
    3169:   PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
    3170:   PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats));
    3171:   PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
    3172:   PetscCall(PetscMalloc1(maxDof, &rows));
    3173:   PetscCall(PetscMalloc1(maxDof * maxDof, &cols));
    3174:   for (p = pRefStart; p < pRefEnd; p++) {
    3175:     PetscInt parent, pDof, parentDof;

    3177:     PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
    3178:     PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
    3179:     PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
    3180:     if (!pDof || !parentDof || parent == p) continue;

    3182:     PetscCall(PetscMalloc1(numFields, &refPointFieldMats[p - pRefStart]));
    3183:     for (f = 0; f < numFields; f++) {
    3184:       PetscInt cDof, cOff, numCols, r;

    3186:       if (numFields > 1) {
    3187:         PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
    3188:         PetscCall(PetscSectionGetFieldOffset(refConSec, p, f, &cOff));
    3189:       } else {
    3190:         PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
    3191:         PetscCall(PetscSectionGetOffset(refConSec, p, &cOff));
    3192:       }

    3194:       for (r = 0; r < cDof; r++) rows[r] = cOff + r;
    3195:       numCols = 0;
    3196:       {
    3197:         PetscInt aDof, aOff, j;

    3199:         if (numFields > 1) {
    3200:           PetscCall(PetscSectionGetFieldDof(refSection, parent, f, &aDof));
    3201:           PetscCall(PetscSectionGetFieldOffset(refSection, parent, f, &aOff));
    3202:         } else {
    3203:           PetscCall(PetscSectionGetDof(refSection, parent, &aDof));
    3204:           PetscCall(PetscSectionGetOffset(refSection, parent, &aOff));
    3205:         }

    3207:         for (j = 0; j < aDof; j++) cols[numCols++] = aOff + j;
    3208:       }
    3209:       PetscCall(PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]));
    3210:       /* transpose of constraint matrix */
    3211:       PetscCall(MatGetValues(inj, numCols, cols, cDof, rows, refPointFieldMats[p - pRefStart][f]));
    3212:     }
    3213:   }
    3214:   *childrenMats = refPointFieldMats;
    3215:   PetscCall(PetscFree(rows));
    3216:   PetscCall(PetscFree(cols));
    3217:   PetscFunctionReturn(PETSC_SUCCESS);
    3218: }

    3220: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
    3221: {
    3222:   PetscDS        ds;
    3223:   PetscScalar ***refPointFieldMats;
    3224:   PetscInt       numFields, pRefStart, pRefEnd, p, f;
    3225:   PetscSection   refConSec, refSection;

    3227:   PetscFunctionBegin;
    3228:   refPointFieldMats = *childrenMats;
    3229:   *childrenMats     = NULL;
    3230:   PetscCall(DMGetDS(refTree, &ds));
    3231:   PetscCall(DMGetLocalSection(refTree, &refSection));
    3232:   PetscCall(PetscDSGetNumFields(ds, &numFields));
    3233:   PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
    3234:   PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
    3235:   for (p = pRefStart; p < pRefEnd; p++) {
    3236:     PetscInt parent, pDof, parentDof;

    3238:     PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
    3239:     PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
    3240:     PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
    3241:     if (!pDof || !parentDof || parent == p) continue;

    3243:     for (f = 0; f < numFields; f++) {
    3244:       PetscInt cDof;

    3246:       if (numFields > 1) {
    3247:         PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
    3248:       } else {
    3249:         PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
    3250:       }

    3252:       PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f]));
    3253:     }
    3254:     PetscCall(PetscFree(refPointFieldMats[p - pRefStart]));
    3255:   }
    3256:   PetscCall(PetscFree(refPointFieldMats));
    3257:   PetscFunctionReturn(PETSC_SUCCESS);
    3258: }

    3260: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree, Mat *injRef)
    3261: {
    3262:   Mat         cMatRef;
    3263:   PetscObject injRefObj;

    3265:   PetscFunctionBegin;
    3266:   PetscCall(DMGetDefaultConstraints(refTree, NULL, &cMatRef, NULL));
    3267:   PetscCall(PetscObjectQuery((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", &injRefObj));
    3268:   *injRef = (Mat)injRefObj;
    3269:   if (!*injRef) {
    3270:     PetscCall(DMPlexComputeInjectorReferenceTree(refTree, injRef));
    3271:     PetscCall(PetscObjectCompose((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", (PetscObject)*injRef));
    3272:     /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
    3273:     PetscCall(PetscObjectDereference((PetscObject)*injRef));
    3274:   }
    3275:   PetscFunctionReturn(PETSC_SUCCESS);
    3276: }

    3278: static PetscErrorCode DMPlexTransferInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, const PetscInt *childIds, Vec fineVec, PetscInt numFields, PetscInt *offsets, PetscSection *rootMultiSec, PetscSection *multiLeafSec, PetscInt **gatheredIndices, PetscScalar **gatheredValues)
    3279: {
    3280:   PetscInt        pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
    3281:   PetscSection    globalCoarse, globalFine;
    3282:   PetscSection    localCoarse, localFine, leafIndicesSec;
    3283:   PetscSection    multiRootSec, rootIndicesSec;
    3284:   PetscInt       *leafInds, *rootInds = NULL;
    3285:   const PetscInt *rootDegrees;
    3286:   PetscScalar    *leafVals = NULL, *rootVals = NULL;
    3287:   PetscSF         coarseToFineEmbedded;

    3289:   PetscFunctionBegin;
    3290:   PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
    3291:   PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
    3292:   PetscCall(DMGetLocalSection(fine, &localFine));
    3293:   PetscCall(DMGetGlobalSection(fine, &globalFine));
    3294:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
    3295:   PetscCall(PetscSectionSetChart(leafIndicesSec, pStartF, pEndF));
    3296:   PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
    3297:   { /* winnow fine points that don't have global dofs out of the sf */
    3298:     PetscInt        l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
    3299:     const PetscInt *leaves;

    3301:     PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
    3302:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
    3303:       p = leaves ? leaves[l] : l;
    3304:       PetscCall(PetscSectionGetDof(globalFine, p, &dof));
    3305:       PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
    3306:       if ((dof - cdof) > 0) {
    3307:         numPointsWithDofs++;

    3309:         PetscCall(PetscSectionGetDof(localFine, p, &dof));
    3310:         PetscCall(PetscSectionSetDof(leafIndicesSec, p, dof + 1));
    3311:       }
    3312:     }
    3313:     PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
    3314:     PetscCall(PetscSectionSetUp(leafIndicesSec));
    3315:     PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numIndices));
    3316:     PetscCall(PetscMalloc1((gatheredIndices ? numIndices : (maxDof + 1)), &leafInds));
    3317:     if (gatheredValues) PetscCall(PetscMalloc1(numIndices, &leafVals));
    3318:     for (l = 0, offset = 0; l < nleaves; l++) {
    3319:       p = leaves ? leaves[l] : l;
    3320:       PetscCall(PetscSectionGetDof(globalFine, p, &dof));
    3321:       PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
    3322:       if ((dof - cdof) > 0) {
    3323:         PetscInt     off, gOff;
    3324:         PetscInt    *pInd;
    3325:         PetscScalar *pVal = NULL;

    3327:         pointsWithDofs[offset++] = l;

    3329:         PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));

    3331:         pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
    3332:         if (gatheredValues) {
    3333:           PetscInt i;

    3335:           pVal = &leafVals[off + 1];
    3336:           for (i = 0; i < dof; i++) pVal[i] = 0.;
    3337:         }
    3338:         PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));

    3340:         offsets[0] = 0;
    3341:         if (numFields) {
    3342:           PetscInt f;

    3344:           for (f = 0; f < numFields; f++) {
    3345:             PetscInt fDof;
    3346:             PetscCall(PetscSectionGetFieldDof(localFine, p, f, &fDof));
    3347:             offsets[f + 1] = fDof + offsets[f];
    3348:           }
    3349:           PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
    3350:         } else {
    3351:           PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
    3352:         }
    3353:         if (gatheredValues) PetscCall(VecGetValues(fineVec, dof, pInd, pVal));
    3354:       }
    3355:     }
    3356:     PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
    3357:     PetscCall(PetscFree(pointsWithDofs));
    3358:   }

    3360:   PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
    3361:   PetscCall(DMGetLocalSection(coarse, &localCoarse));
    3362:   PetscCall(DMGetGlobalSection(coarse, &globalCoarse));

    3364:   { /* there may be the case where an sf root has a parent: broadcast parents back to children */
    3365:     MPI_Datatype threeInt;
    3366:     PetscMPIInt  rank;
    3367:     PetscInt (*parentNodeAndIdCoarse)[3];
    3368:     PetscInt (*parentNodeAndIdFine)[3];
    3369:     PetscInt           p, nleaves, nleavesToParents;
    3370:     PetscSF            pointSF, sfToParents;
    3371:     const PetscInt    *ilocal;
    3372:     const PetscSFNode *iremote;
    3373:     PetscSFNode       *iremoteToParents;
    3374:     PetscInt          *ilocalToParents;

    3376:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
    3377:     PetscCallMPI(MPI_Type_contiguous(3, MPIU_INT, &threeInt));
    3378:     PetscCallMPI(MPI_Type_commit(&threeInt));
    3379:     PetscCall(PetscMalloc2(pEndC - pStartC, &parentNodeAndIdCoarse, pEndF - pStartF, &parentNodeAndIdFine));
    3380:     PetscCall(DMGetPointSF(coarse, &pointSF));
    3381:     PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, &iremote));
    3382:     for (p = pStartC; p < pEndC; p++) {
    3383:       PetscInt parent, childId;
    3384:       PetscCall(DMPlexGetTreeParent(coarse, p, &parent, &childId));
    3385:       parentNodeAndIdCoarse[p - pStartC][0] = rank;
    3386:       parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
    3387:       parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
    3388:       if (nleaves > 0) {
    3389:         PetscInt leaf = -1;

    3391:         if (ilocal) {
    3392:           PetscCall(PetscFindInt(parent, nleaves, ilocal, &leaf));
    3393:         } else {
    3394:           leaf = p - pStartC;
    3395:         }
    3396:         if (leaf >= 0) {
    3397:           parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
    3398:           parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
    3399:         }
    3400:       }
    3401:     }
    3402:     for (p = pStartF; p < pEndF; p++) {
    3403:       parentNodeAndIdFine[p - pStartF][0] = -1;
    3404:       parentNodeAndIdFine[p - pStartF][1] = -1;
    3405:       parentNodeAndIdFine[p - pStartF][2] = -1;
    3406:     }
    3407:     PetscCall(PetscSFBcastBegin(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
    3408:     PetscCall(PetscSFBcastEnd(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
    3409:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
    3410:       PetscInt dof;

    3412:       PetscCall(PetscSectionGetDof(leafIndicesSec, p, &dof));
    3413:       if (dof) {
    3414:         PetscInt off;

    3416:         PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));
    3417:         if (gatheredIndices) {
    3418:           leafInds[off] = PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
    3419:         } else if (gatheredValues) {
    3420:           leafVals[off] = (PetscScalar)PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
    3421:         }
    3422:       }
    3423:       if (parentNodeAndIdFine[p - pStartF][0] >= 0) nleavesToParents++;
    3424:     }
    3425:     PetscCall(PetscMalloc1(nleavesToParents, &ilocalToParents));
    3426:     PetscCall(PetscMalloc1(nleavesToParents, &iremoteToParents));
    3427:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
    3428:       if (parentNodeAndIdFine[p - pStartF][0] >= 0) {
    3429:         ilocalToParents[nleavesToParents] = p - pStartF;
    3430:         // FIXME PetscCall(PetscMPIIntCast(parentNodeAndIdFine[p - pStartF][0],&iremoteToParents[nleavesToParents].rank));
    3431:         iremoteToParents[nleavesToParents].rank  = parentNodeAndIdFine[p - pStartF][0];
    3432:         iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p - pStartF][1];
    3433:         nleavesToParents++;
    3434:       }
    3435:     }
    3436:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)coarse), &sfToParents));
    3437:     PetscCall(PetscSFSetGraph(sfToParents, pEndC - pStartC, nleavesToParents, ilocalToParents, PETSC_OWN_POINTER, iremoteToParents, PETSC_OWN_POINTER));
    3438:     PetscCall(PetscSFDestroy(&coarseToFineEmbedded));

    3440:     coarseToFineEmbedded = sfToParents;

    3442:     PetscCall(PetscFree2(parentNodeAndIdCoarse, parentNodeAndIdFine));
    3443:     PetscCallMPI(MPI_Type_free(&threeInt));
    3444:   }

    3446:   { /* winnow out coarse points that don't have dofs */
    3447:     PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
    3448:     PetscSF  sfDofsOnly;

    3450:     for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
    3451:       PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
    3452:       PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
    3453:       if ((dof - cdof) > 0) numPointsWithDofs++;
    3454:     }
    3455:     PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
    3456:     for (p = pStartC, offset = 0; p < pEndC; p++) {
    3457:       PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
    3458:       PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
    3459:       if ((dof - cdof) > 0) pointsWithDofs[offset++] = p - pStartC;
    3460:     }
    3461:     PetscCall(PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly));
    3462:     PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
    3463:     PetscCall(PetscFree(pointsWithDofs));
    3464:     coarseToFineEmbedded = sfDofsOnly;
    3465:   }

    3467:   /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
    3468:   PetscCall(PetscSFComputeDegreeBegin(coarseToFineEmbedded, &rootDegrees));
    3469:   PetscCall(PetscSFComputeDegreeEnd(coarseToFineEmbedded, &rootDegrees));
    3470:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &multiRootSec));
    3471:   PetscCall(PetscSectionSetChart(multiRootSec, pStartC, pEndC));
    3472:   for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionSetDof(multiRootSec, p, rootDegrees[p - pStartC]));
    3473:   PetscCall(PetscSectionSetUp(multiRootSec));
    3474:   PetscCall(PetscSectionGetStorageSize(multiRootSec, &numMulti));
    3475:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
    3476:   { /* distribute the leaf section */
    3477:     PetscSF   multi, multiInv, indicesSF;
    3478:     PetscInt *remoteOffsets, numRootIndices;

    3480:     PetscCall(PetscSFGetMultiSF(coarseToFineEmbedded, &multi));
    3481:     PetscCall(PetscSFCreateInverseSF(multi, &multiInv));
    3482:     PetscCall(PetscSFDistributeSection(multiInv, leafIndicesSec, &remoteOffsets, rootIndicesSec));
    3483:     PetscCall(PetscSFCreateSectionSF(multiInv, leafIndicesSec, remoteOffsets, rootIndicesSec, &indicesSF));
    3484:     PetscCall(PetscFree(remoteOffsets));
    3485:     PetscCall(PetscSFDestroy(&multiInv));
    3486:     PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
    3487:     if (gatheredIndices) {
    3488:       PetscCall(PetscMalloc1(numRootIndices, &rootInds));
    3489:       PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
    3490:       PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
    3491:     }
    3492:     if (gatheredValues) {
    3493:       PetscCall(PetscMalloc1(numRootIndices, &rootVals));
    3494:       PetscCall(PetscSFBcastBegin(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
    3495:       PetscCall(PetscSFBcastEnd(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
    3496:     }
    3497:     PetscCall(PetscSFDestroy(&indicesSF));
    3498:   }
    3499:   PetscCall(PetscSectionDestroy(&leafIndicesSec));
    3500:   PetscCall(PetscFree(leafInds));
    3501:   PetscCall(PetscFree(leafVals));
    3502:   PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
    3503:   *rootMultiSec = multiRootSec;
    3504:   *multiLeafSec = rootIndicesSec;
    3505:   if (gatheredIndices) *gatheredIndices = rootInds;
    3506:   if (gatheredValues) *gatheredValues = rootVals;
    3507:   PetscFunctionReturn(PETSC_SUCCESS);
    3508: }

    3510: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
    3511: {
    3512:   DM             refTree;
    3513:   PetscSection   multiRootSec, rootIndicesSec;
    3514:   PetscSection   globalCoarse, globalFine;
    3515:   PetscSection   localCoarse, localFine;
    3516:   PetscSection   cSecRef;
    3517:   PetscInt      *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd;
    3518:   Mat            injRef;
    3519:   PetscInt       numFields, maxDof;
    3520:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
    3521:   PetscInt      *offsets, *offsetsCopy, *rowOffsets;
    3522:   PetscLayout    rowMap, colMap;
    3523:   PetscInt       rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
    3524:   PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */

    3526:   PetscFunctionBegin;
    3527:   /* get the templates for the fine-to-coarse injection from the reference tree */
    3528:   PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
    3529:   PetscCall(DMCopyDisc(coarse, refTree));
    3530:   PetscCall(DMSetLocalSection(refTree, NULL));
    3531:   PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
    3532:   PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
    3533:   PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
    3534:   PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));

    3536:   PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
    3537:   PetscCall(DMGetLocalSection(fine, &localFine));
    3538:   PetscCall(DMGetGlobalSection(fine, &globalFine));
    3539:   PetscCall(PetscSectionGetNumFields(localFine, &numFields));
    3540:   PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
    3541:   PetscCall(DMGetLocalSection(coarse, &localCoarse));
    3542:   PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
    3543:   PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
    3544:   {
    3545:     PetscInt maxFields = PetscMax(1, numFields) + 1;
    3546:     PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
    3547:   }

    3549:   PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, childIds, NULL, numFields, offsets, &multiRootSec, &rootIndicesSec, &rootIndices, NULL));

    3551:   PetscCall(PetscMalloc1(maxDof, &parentIndices));

    3553:   /* count indices */
    3554:   PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
    3555:   PetscCall(PetscLayoutSetUp(rowMap));
    3556:   PetscCall(PetscLayoutSetUp(colMap));
    3557:   PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
    3558:   PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
    3559:   PetscCall(PetscCalloc2(rowEnd - rowStart, &nnzD, rowEnd - rowStart, &nnzO));
    3560:   for (p = pStartC; p < pEndC; p++) {
    3561:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

    3563:     PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
    3564:     PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
    3565:     if ((dof - cdof) <= 0) continue;
    3566:     PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));

    3568:     rowOffsets[0]  = 0;
    3569:     offsetsCopy[0] = 0;
    3570:     if (numFields) {
    3571:       PetscInt f;

    3573:       for (f = 0; f < numFields; f++) {
    3574:         PetscInt fDof;
    3575:         PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
    3576:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
    3577:       }
    3578:       PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
    3579:     } else {
    3580:       PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
    3581:       rowOffsets[1] = offsetsCopy[0];
    3582:     }

    3584:     PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
    3585:     PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
    3586:     leafEnd = leafStart + numLeaves;
    3587:     for (l = leafStart; l < leafEnd; l++) {
    3588:       PetscInt        numIndices, childId, offset;
    3589:       const PetscInt *childIndices;

    3591:       PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
    3592:       PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
    3593:       childId      = rootIndices[offset++];
    3594:       childIndices = &rootIndices[offset];
    3595:       numIndices--;

    3597:       if (childId == -1) { /* equivalent points: scatter */
    3598:         PetscInt i;

    3600:         for (i = 0; i < numIndices; i++) {
    3601:           PetscInt colIndex = childIndices[i];
    3602:           PetscInt rowIndex = parentIndices[i];
    3603:           if (rowIndex < 0) continue;
    3604:           PetscCheck(colIndex >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unconstrained fine and constrained coarse");
    3605:           if (colIndex >= colStart && colIndex < colEnd) {
    3606:             nnzD[rowIndex - rowStart] = 1;
    3607:           } else {
    3608:             nnzO[rowIndex - rowStart] = 1;
    3609:           }
    3610:         }
    3611:       } else {
    3612:         PetscInt parentId, f, lim;

    3614:         PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));

    3616:         lim        = PetscMax(1, numFields);
    3617:         offsets[0] = 0;
    3618:         if (numFields) {
    3619:           PetscInt f;

    3621:           for (f = 0; f < numFields; f++) {
    3622:             PetscInt fDof;
    3623:             PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));

    3625:             offsets[f + 1] = fDof + offsets[f];
    3626:           }
    3627:         } else {
    3628:           PetscInt cDof;

    3630:           PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
    3631:           offsets[1] = cDof;
    3632:         }
    3633:         for (f = 0; f < lim; f++) {
    3634:           PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
    3635:           PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
    3636:           PetscInt i, numD = 0, numO = 0;

    3638:           for (i = childStart; i < childEnd; i++) {
    3639:             PetscInt colIndex = childIndices[i];

    3641:             if (colIndex < 0) continue;
    3642:             if (colIndex >= colStart && colIndex < colEnd) {
    3643:               numD++;
    3644:             } else {
    3645:               numO++;
    3646:             }
    3647:           }
    3648:           for (i = parentStart; i < parentEnd; i++) {
    3649:             PetscInt rowIndex = parentIndices[i];

    3651:             if (rowIndex < 0) continue;
    3652:             nnzD[rowIndex - rowStart] += numD;
    3653:             nnzO[rowIndex - rowStart] += numO;
    3654:           }
    3655:         }
    3656:       }
    3657:     }
    3658:   }
    3659:   /* preallocate */
    3660:   PetscCall(MatXAIJSetPreallocation(mat, 1, nnzD, nnzO, NULL, NULL));
    3661:   PetscCall(PetscFree2(nnzD, nnzO));
    3662:   /* insert values */
    3663:   PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
    3664:   for (p = pStartC; p < pEndC; p++) {
    3665:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

    3667:     PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
    3668:     PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
    3669:     if ((dof - cdof) <= 0) continue;
    3670:     PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));

    3672:     rowOffsets[0]  = 0;
    3673:     offsetsCopy[0] = 0;
    3674:     if (numFields) {
    3675:       PetscInt f;

    3677:       for (f = 0; f < numFields; f++) {
    3678:         PetscInt fDof;
    3679:         PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
    3680:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
    3681:       }
    3682:       PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
    3683:     } else {
    3684:       PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
    3685:       rowOffsets[1] = offsetsCopy[0];
    3686:     }

    3688:     PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
    3689:     PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
    3690:     leafEnd = leafStart + numLeaves;
    3691:     for (l = leafStart; l < leafEnd; l++) {
    3692:       PetscInt        numIndices, childId, offset;
    3693:       const PetscInt *childIndices;

    3695:       PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
    3696:       PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
    3697:       childId      = rootIndices[offset++];
    3698:       childIndices = &rootIndices[offset];
    3699:       numIndices--;

    3701:       if (childId == -1) { /* equivalent points: scatter */
    3702:         PetscInt i;

    3704:         for (i = 0; i < numIndices; i++) PetscCall(MatSetValue(mat, parentIndices[i], childIndices[i], 1., INSERT_VALUES));
    3705:       } else {
    3706:         PetscInt parentId, f, lim;

    3708:         PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));

    3710:         lim        = PetscMax(1, numFields);
    3711:         offsets[0] = 0;
    3712:         if (numFields) {
    3713:           PetscInt f;

    3715:           for (f = 0; f < numFields; f++) {
    3716:             PetscInt fDof;
    3717:             PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));

    3719:             offsets[f + 1] = fDof + offsets[f];
    3720:           }
    3721:         } else {
    3722:           PetscInt cDof;

    3724:           PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
    3725:           offsets[1] = cDof;
    3726:         }
    3727:         for (f = 0; f < lim; f++) {
    3728:           PetscScalar    *childMat   = &childrenMats[childId - pRefStart][f][0];
    3729:           PetscInt       *rowIndices = &parentIndices[rowOffsets[f]];
    3730:           const PetscInt *colIndices = &childIndices[offsets[f]];

    3732:           PetscCall(MatSetValues(mat, rowOffsets[f + 1] - rowOffsets[f], rowIndices, offsets[f + 1] - offsets[f], colIndices, childMat, INSERT_VALUES));
    3733:         }
    3734:       }
    3735:     }
    3736:   }
    3737:   PetscCall(PetscSectionDestroy(&multiRootSec));
    3738:   PetscCall(PetscSectionDestroy(&rootIndicesSec));
    3739:   PetscCall(PetscFree(parentIndices));
    3740:   PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
    3741:   PetscCall(PetscFree(rootIndices));
    3742:   PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));

    3744:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
    3745:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
    3746:   PetscFunctionReturn(PETSC_SUCCESS);
    3747: }

    3749: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
    3750: {
    3751:   PetscSF            coarseToFineEmbedded;
    3752:   PetscSection       globalCoarse, globalFine;
    3753:   PetscSection       localCoarse, localFine;
    3754:   PetscSection       aSec, cSec;
    3755:   PetscSection       rootValuesSec;
    3756:   PetscSection       leafValuesSec;
    3757:   PetscScalar       *rootValues, *leafValues;
    3758:   IS                 aIS;
    3759:   const PetscInt    *anchors;
    3760:   Mat                cMat;
    3761:   PetscInt           numFields;
    3762:   PetscInt           pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd;
    3763:   PetscInt           aStart, aEnd, cStart, cEnd;
    3764:   PetscInt          *maxChildIds;
    3765:   PetscInt          *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
    3766:   PetscFV            fv = NULL;
    3767:   PetscInt           dim, numFVcomps = -1, fvField = -1;
    3768:   DM                 cellDM = NULL, gradDM = NULL;
    3769:   const PetscScalar *cellGeomArray = NULL;
    3770:   const PetscScalar *gradArray     = NULL;

    3772:   PetscFunctionBegin;
    3773:   PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
    3774:   PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
    3775:   PetscCall(DMPlexGetSimplexOrBoxCells(coarse, 0, &cellStart, &cellEnd));
    3776:   PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
    3777:   PetscCall(DMGetGlobalSection(fine, &globalFine));
    3778:   PetscCall(DMGetCoordinateDim(coarse, &dim));
    3779:   { /* winnow fine points that don't have global dofs out of the sf */
    3780:     PetscInt        nleaves, l;
    3781:     const PetscInt *leaves;
    3782:     PetscInt        dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;

    3784:     PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));

    3786:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
    3787:       PetscInt p = leaves ? leaves[l] : l;

    3789:       PetscCall(PetscSectionGetDof(globalFine, p, &dof));
    3790:       PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
    3791:       if ((dof - cdof) > 0) numPointsWithDofs++;
    3792:     }
    3793:     PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
    3794:     for (l = 0, offset = 0; l < nleaves; l++) {
    3795:       PetscInt p = leaves ? leaves[l] : l;

    3797:       PetscCall(PetscSectionGetDof(globalFine, p, &dof));
    3798:       PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
    3799:       if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
    3800:     }
    3801:     PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
    3802:     PetscCall(PetscFree(pointsWithDofs));
    3803:   }
    3804:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
    3805:   PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
    3806:   for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
    3807:   PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));
    3808:   PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));

    3810:   PetscCall(DMGetLocalSection(coarse, &localCoarse));
    3811:   PetscCall(DMGetGlobalSection(coarse, &globalCoarse));

    3813:   PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
    3814:   PetscCall(ISGetIndices(aIS, &anchors));
    3815:   PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));

    3817:   PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
    3818:   PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));

    3820:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
    3821:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootValuesSec));
    3822:   PetscCall(PetscSectionSetChart(rootValuesSec, pStartC, pEndC));
    3823:   PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
    3824:   {
    3825:     PetscInt maxFields = PetscMax(1, numFields) + 1;
    3826:     PetscCall(PetscMalloc7(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &newOffsets, maxFields, &newOffsetsCopy, maxFields, &rowOffsets, maxFields, &numD, maxFields, &numO));
    3827:   }
    3828:   if (grad) {
    3829:     PetscInt i;

    3831:     PetscCall(VecGetDM(cellGeom, &cellDM));
    3832:     PetscCall(VecGetArrayRead(cellGeom, &cellGeomArray));
    3833:     PetscCall(VecGetDM(grad, &gradDM));
    3834:     PetscCall(VecGetArrayRead(grad, &gradArray));
    3835:     for (i = 0; i < PetscMax(1, numFields); i++) {
    3836:       PetscObject  obj;
    3837:       PetscClassId id;

    3839:       PetscCall(DMGetField(coarse, i, NULL, &obj));
    3840:       PetscCall(PetscObjectGetClassId(obj, &id));
    3841:       if (id == PETSCFV_CLASSID) {
    3842:         fv = (PetscFV)obj;
    3843:         PetscCall(PetscFVGetNumComponents(fv, &numFVcomps));
    3844:         fvField = i;
    3845:         break;
    3846:       }
    3847:     }
    3848:   }

    3850:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
    3851:     PetscInt dof;
    3852:     PetscInt maxChildId = maxChildIds[p - pStartC];
    3853:     PetscInt numValues  = 0;

    3855:     PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
    3856:     if (dof < 0) dof = -(dof + 1);
    3857:     offsets[0]    = 0;
    3858:     newOffsets[0] = 0;
    3859:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
    3860:       PetscInt *closure = NULL, closureSize, cl;

    3862:       PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
    3863:       for (cl = 0; cl < closureSize; cl++) { /* get the closure */
    3864:         PetscInt c = closure[2 * cl], clDof;

    3866:         PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
    3867:         numValues += clDof;
    3868:       }
    3869:       PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
    3870:     } else if (maxChildId == -1) {
    3871:       PetscCall(PetscSectionGetDof(localCoarse, p, &numValues));
    3872:     }
    3873:     /* we will pack the column indices with the field offsets */
    3874:     if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
    3875:       /* also send the centroid, and the gradient */
    3876:       numValues += dim * (1 + numFVcomps);
    3877:     }
    3878:     PetscCall(PetscSectionSetDof(rootValuesSec, p, numValues));
    3879:   }
    3880:   PetscCall(PetscSectionSetUp(rootValuesSec));
    3881:   {
    3882:     PetscInt           numRootValues;
    3883:     const PetscScalar *coarseArray;

    3885:     PetscCall(PetscSectionGetStorageSize(rootValuesSec, &numRootValues));
    3886:     PetscCall(PetscMalloc1(numRootValues, &rootValues));
    3887:     PetscCall(VecGetArrayRead(vecCoarseLocal, &coarseArray));
    3888:     for (p = pStartC; p < pEndC; p++) {
    3889:       PetscInt     numValues;
    3890:       PetscInt     pValOff;
    3891:       PetscScalar *pVal;
    3892:       PetscInt     maxChildId = maxChildIds[p - pStartC];

    3894:       PetscCall(PetscSectionGetDof(rootValuesSec, p, &numValues));
    3895:       if (!numValues) continue;
    3896:       PetscCall(PetscSectionGetOffset(rootValuesSec, p, &pValOff));
    3897:       pVal = &rootValues[pValOff];
    3898:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
    3899:         PetscInt closureSize = numValues;
    3900:         PetscCall(DMPlexVecGetClosure(coarse, NULL, vecCoarseLocal, p, &closureSize, &pVal));
    3901:         if (grad && p >= cellStart && p < cellEnd) {
    3902:           PetscFVCellGeom *cg;
    3903:           PetscScalar     *gradVals = NULL;
    3904:           PetscInt         i;

    3906:           pVal += (numValues - dim * (1 + numFVcomps));

    3908:           PetscCall(DMPlexPointLocalRead(cellDM, p, cellGeomArray, (void *)&cg));
    3909:           for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
    3910:           pVal += dim;
    3911:           PetscCall(DMPlexPointGlobalRead(gradDM, p, gradArray, (void *)&gradVals));
    3912:           for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
    3913:         }
    3914:       } else if (maxChildId == -1) {
    3915:         PetscInt lDof, lOff, i;

    3917:         PetscCall(PetscSectionGetDof(localCoarse, p, &lDof));
    3918:         PetscCall(PetscSectionGetOffset(localCoarse, p, &lOff));
    3919:         for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
    3920:       }
    3921:     }
    3922:     PetscCall(VecRestoreArrayRead(vecCoarseLocal, &coarseArray));
    3923:     PetscCall(PetscFree(maxChildIds));
    3924:   }
    3925:   {
    3926:     PetscSF   valuesSF;
    3927:     PetscInt *remoteOffsetsValues, numLeafValues;

    3929:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafValuesSec));
    3930:     PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootValuesSec, &remoteOffsetsValues, leafValuesSec));
    3931:     PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootValuesSec, remoteOffsetsValues, leafValuesSec, &valuesSF));
    3932:     PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
    3933:     PetscCall(PetscFree(remoteOffsetsValues));
    3934:     PetscCall(PetscSectionGetStorageSize(leafValuesSec, &numLeafValues));
    3935:     PetscCall(PetscMalloc1(numLeafValues, &leafValues));
    3936:     PetscCall(PetscSFBcastBegin(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
    3937:     PetscCall(PetscSFBcastEnd(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
    3938:     PetscCall(PetscSFDestroy(&valuesSF));
    3939:     PetscCall(PetscFree(rootValues));
    3940:     PetscCall(PetscSectionDestroy(&rootValuesSec));
    3941:   }
    3942:   PetscCall(DMGetLocalSection(fine, &localFine));
    3943:   {
    3944:     PetscInt       maxDof;
    3945:     PetscInt      *rowIndices;
    3946:     DM             refTree;
    3947:     PetscInt     **refPointFieldN;
    3948:     PetscScalar ***refPointFieldMats;
    3949:     PetscSection   refConSec, refAnSec;
    3950:     PetscInt       pRefStart, pRefEnd, leafStart, leafEnd;
    3951:     PetscScalar   *pointWork;

    3953:     PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
    3954:     PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
    3955:     PetscCall(DMGetWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
    3956:     PetscCall(DMPlexGetReferenceTree(fine, &refTree));
    3957:     PetscCall(DMCopyDisc(fine, refTree));
    3958:     PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
    3959:     PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
    3960:     PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL));
    3961:     PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
    3962:     PetscCall(PetscSectionGetChart(leafValuesSec, &leafStart, &leafEnd));
    3963:     PetscCall(DMPlexGetSimplexOrBoxCells(fine, 0, &cellStart, &cellEnd));
    3964:     for (p = leafStart; p < leafEnd; p++) {
    3965:       PetscInt           gDof, gcDof, gOff, lDof;
    3966:       PetscInt           numValues, pValOff;
    3967:       PetscInt           childId;
    3968:       const PetscScalar *pVal;
    3969:       const PetscScalar *fvGradData = NULL;

    3971:       PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
    3972:       PetscCall(PetscSectionGetDof(localFine, p, &lDof));
    3973:       PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
    3974:       if ((gDof - gcDof) <= 0) continue;
    3975:       PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
    3976:       PetscCall(PetscSectionGetDof(leafValuesSec, p, &numValues));
    3977:       if (!numValues) continue;
    3978:       PetscCall(PetscSectionGetOffset(leafValuesSec, p, &pValOff));
    3979:       pVal              = &leafValues[pValOff];
    3980:       offsets[0]        = 0;
    3981:       offsetsCopy[0]    = 0;
    3982:       newOffsets[0]     = 0;
    3983:       newOffsetsCopy[0] = 0;
    3984:       childId           = cids[p - pStartF];
    3985:       if (numFields) {
    3986:         PetscInt f;
    3987:         for (f = 0; f < numFields; f++) {
    3988:           PetscInt rowDof;

    3990:           PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
    3991:           offsets[f + 1]     = offsets[f] + rowDof;
    3992:           offsetsCopy[f + 1] = offsets[f + 1];
    3993:           /* TODO: closure indices */
    3994:           newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
    3995:         }
    3996:         PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
    3997:       } else {
    3998:         offsets[0]    = 0;
    3999:         offsets[1]    = lDof;
    4000:         newOffsets[0] = 0;
    4001:         newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
    4002:         PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
    4003:       }
    4004:       if (childId == -1) { /* no child interpolation: one nnz per */
    4005:         PetscCall(VecSetValues(vecFine, numValues, rowIndices, pVal, INSERT_VALUES));
    4006:       } else {
    4007:         PetscInt f;

    4009:         if (grad && p >= cellStart && p < cellEnd) {
    4010:           numValues -= (dim * (1 + numFVcomps));
    4011:           fvGradData = &pVal[numValues];
    4012:         }
    4013:         for (f = 0; f < PetscMax(1, numFields); f++) {
    4014:           const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
    4015:           PetscInt           numRows  = offsets[f + 1] - offsets[f];
    4016:           PetscInt           numCols  = newOffsets[f + 1] - newOffsets[f];
    4017:           const PetscScalar *cVal     = &pVal[newOffsets[f]];
    4018:           PetscScalar       *rVal     = &pointWork[offsets[f]];
    4019:           PetscInt           i, j;

    4021: #if 0
    4022:           PetscCall(PetscInfo(coarse,"childId %" PetscInt_FMT ", numRows %" PetscInt_FMT ", numCols %" PetscInt_FMT ", refPointFieldN %" PetscInt_FMT " maxDof %" PetscInt_FMT "\n",childId,numRows,numCols,refPointFieldN[childId - pRefStart][f], maxDof));
    4023: #endif
    4024:           for (i = 0; i < numRows; i++) {
    4025:             PetscScalar val = 0.;
    4026:             for (j = 0; j < numCols; j++) val += childMat[i * numCols + j] * cVal[j];
    4027:             rVal[i] = val;
    4028:           }
    4029:           if (f == fvField && p >= cellStart && p < cellEnd) {
    4030:             PetscReal          centroid[3];
    4031:             PetscScalar        diff[3];
    4032:             const PetscScalar *parentCentroid = &fvGradData[0];
    4033:             const PetscScalar *gradient       = &fvGradData[dim];

    4035:             PetscCall(DMPlexComputeCellGeometryFVM(fine, p, NULL, centroid, NULL));
    4036:             for (i = 0; i < dim; i++) diff[i] = centroid[i] - parentCentroid[i];
    4037:             for (i = 0; i < numFVcomps; i++) {
    4038:               PetscScalar val = 0.;

    4040:               for (j = 0; j < dim; j++) val += gradient[dim * i + j] * diff[j];
    4041:               rVal[i] += val;
    4042:             }
    4043:           }
    4044:           PetscCall(VecSetValues(vecFine, numRows, &rowIndices[offsets[f]], rVal, INSERT_VALUES));
    4045:         }
    4046:       }
    4047:     }
    4048:     PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
    4049:     PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
    4050:     PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
    4051:   }
    4052:   PetscCall(PetscFree(leafValues));
    4053:   PetscCall(PetscSectionDestroy(&leafValuesSec));
    4054:   PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
    4055:   PetscCall(ISRestoreIndices(aIS, &anchors));
    4056:   PetscFunctionReturn(PETSC_SUCCESS);
    4057: }

    4059: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
    4060: {
    4061:   DM             refTree;
    4062:   PetscSection   multiRootSec, rootIndicesSec;
    4063:   PetscSection   globalCoarse, globalFine;
    4064:   PetscSection   localCoarse, localFine;
    4065:   PetscSection   cSecRef;
    4066:   PetscInt      *parentIndices, pRefStart, pRefEnd;
    4067:   PetscScalar   *rootValues, *parentValues;
    4068:   Mat            injRef;
    4069:   PetscInt       numFields, maxDof;
    4070:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
    4071:   PetscInt      *offsets, *offsetsCopy, *rowOffsets;
    4072:   PetscLayout    rowMap, colMap;
    4073:   PetscInt       rowStart, rowEnd, colStart, colEnd;
    4074:   PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */

    4076:   PetscFunctionBegin;
    4077:   /* get the templates for the fine-to-coarse injection from the reference tree */
    4078:   PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
    4079:   PetscCall(VecSetOption(vecCoarse, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
    4080:   PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
    4081:   PetscCall(DMCopyDisc(coarse, refTree));
    4082:   PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
    4083:   PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
    4084:   PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));

    4086:   PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
    4087:   PetscCall(DMGetLocalSection(fine, &localFine));
    4088:   PetscCall(DMGetGlobalSection(fine, &globalFine));
    4089:   PetscCall(PetscSectionGetNumFields(localFine, &numFields));
    4090:   PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
    4091:   PetscCall(DMGetLocalSection(coarse, &localCoarse));
    4092:   PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
    4093:   PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
    4094:   {
    4095:     PetscInt maxFields = PetscMax(1, numFields) + 1;
    4096:     PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
    4097:   }

    4099:   PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, cids, vecFine, numFields, offsets, &multiRootSec, &rootIndicesSec, NULL, &rootValues));

    4101:   PetscCall(PetscMalloc2(maxDof, &parentIndices, maxDof, &parentValues));

    4103:   /* count indices */
    4104:   PetscCall(VecGetLayout(vecFine, &colMap));
    4105:   PetscCall(VecGetLayout(vecCoarse, &rowMap));
    4106:   PetscCall(PetscLayoutSetUp(rowMap));
    4107:   PetscCall(PetscLayoutSetUp(colMap));
    4108:   PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
    4109:   PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
    4110:   /* insert values */
    4111:   PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
    4112:   for (p = pStartC; p < pEndC; p++) {
    4113:     PetscInt  numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
    4114:     PetscBool contribute = PETSC_FALSE;

    4116:     PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
    4117:     PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
    4118:     if ((dof - cdof) <= 0) continue;
    4119:     PetscCall(PetscSectionGetDof(localCoarse, p, &dof));
    4120:     PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));

    4122:     rowOffsets[0]  = 0;
    4123:     offsetsCopy[0] = 0;
    4124:     if (numFields) {
    4125:       PetscInt f;

    4127:       for (f = 0; f < numFields; f++) {
    4128:         PetscInt fDof;
    4129:         PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
    4130:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
    4131:       }
    4132:       PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
    4133:     } else {
    4134:       PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
    4135:       rowOffsets[1] = offsetsCopy[0];
    4136:     }

    4138:     PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
    4139:     PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
    4140:     leafEnd = leafStart + numLeaves;
    4141:     for (l = 0; l < dof; l++) parentValues[l] = 0.;
    4142:     for (l = leafStart; l < leafEnd; l++) {
    4143:       PetscInt           numIndices, childId, offset;
    4144:       const PetscScalar *childValues;

    4146:       PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
    4147:       PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
    4148:       childId     = (PetscInt)PetscRealPart(rootValues[offset++]);
    4149:       childValues = &rootValues[offset];
    4150:       numIndices--;

    4152:       if (childId == -2) { /* skip */
    4153:         continue;
    4154:       } else if (childId == -1) { /* equivalent points: scatter */
    4155:         PetscInt m;

    4157:         contribute = PETSC_TRUE;
    4158:         for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m];
    4159:       } else { /* contributions from children: sum with injectors from reference tree */
    4160:         PetscInt parentId, f, lim;

    4162:         contribute = PETSC_TRUE;
    4163:         PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));

    4165:         lim        = PetscMax(1, numFields);
    4166:         offsets[0] = 0;
    4167:         if (numFields) {
    4168:           PetscInt f;

    4170:           for (f = 0; f < numFields; f++) {
    4171:             PetscInt fDof;
    4172:             PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));

    4174:             offsets[f + 1] = fDof + offsets[f];
    4175:           }
    4176:         } else {
    4177:           PetscInt cDof;

    4179:           PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
    4180:           offsets[1] = cDof;
    4181:         }
    4182:         for (f = 0; f < lim; f++) {
    4183:           PetscScalar       *childMat = &childrenMats[childId - pRefStart][f][0];
    4184:           PetscInt           n        = offsets[f + 1] - offsets[f];
    4185:           PetscInt           m        = rowOffsets[f + 1] - rowOffsets[f];
    4186:           PetscInt           i, j;
    4187:           const PetscScalar *colValues = &childValues[offsets[f]];

    4189:           for (i = 0; i < m; i++) {
    4190:             PetscScalar val = 0.;
    4191:             for (j = 0; j < n; j++) val += childMat[n * i + j] * colValues[j];
    4192:             parentValues[rowOffsets[f] + i] += val;
    4193:           }
    4194:         }
    4195:       }
    4196:     }
    4197:     if (contribute) PetscCall(VecSetValues(vecCoarse, dof, parentIndices, parentValues, INSERT_VALUES));
    4198:   }
    4199:   PetscCall(PetscSectionDestroy(&multiRootSec));
    4200:   PetscCall(PetscSectionDestroy(&rootIndicesSec));
    4201:   PetscCall(PetscFree2(parentIndices, parentValues));
    4202:   PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
    4203:   PetscCall(PetscFree(rootValues));
    4204:   PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));
    4205:   PetscFunctionReturn(PETSC_SUCCESS);
    4206: }

    4208: /*@
    4209:   DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening
    4210:   that can be represented by a common reference tree used by both.  This routine can be used for a combination of
    4211:   coarsening and refinement at the same time.

    4213:   Collective

    4215:   Input Parameters:
    4216: + dmIn        - The `DMPLEX` mesh for the input vector
    4217: . dmOut       - The second `DMPLEX` mesh
    4218: . vecIn       - The input vector
    4219: . sfRefine    - A star forest indicating points in the mesh `dmIn` (roots in the star forest) that are parents to points in
    4220:                 the mesh `dmOut` (leaves in the star forest), i.e. where `dmOut` is more refined than `dmIn`
    4221: . sfCoarsen   - A star forest indicating points in the mesh `dmOut` (roots in the star forest) that are parents to points in
    4222:                 the mesh `dmIn` (leaves in the star forest), i.e. where `dmOut` is more coarsened than `dmIn`
    4223: . cidsRefine  - The childIds of the points in `dmOut`.  These childIds relate back to the reference tree: childid[j] = k implies
    4224:                 that mesh point j of `dmOut` was refined from a point in `dmIn` just as the mesh point k in the reference
    4225:                 tree was refined from its parent.  childid[j] = -1 indicates that the point j in `dmOut` is exactly
    4226:                 equivalent to its root in `dmIn`, so no interpolation is necessary.  childid[j] = -2 indicates that this
    4227:                 point j in `dmOut` is not a leaf of `sfRefine`.
    4228: . cidsCoarsen - The childIds of the points in `dmIn`.  These childIds relate back to the reference tree: childid[j] = k implies
    4229:                 that mesh point j of dmIn coarsens to a point in `dmOut` just as the mesh point k in the reference
    4230:                 tree coarsens to its parent.  childid[j] = -2 indicates that point j in `dmOut` is not a leaf in `sfCoarsen`.
    4231: . useBCs      - `PETSC_TRUE` indicates that boundary values should be inserted into `vecIn` before transfer.
    4232: - time        - Used if boundary values are time dependent.

    4234:   Output Parameter:
    4235: . vecOut - Using interpolation and injection operators calculated on the reference tree, the transferred
    4236:                 projection of `vecIn` from `dmIn` to `dmOut`.  Note that any field discretized with a `PetscFV` finite volume
    4237:                 method that uses gradient reconstruction will use reconstructed gradients when interpolating from
    4238:                 coarse points to fine points.

    4240:   Level: developer

    4242: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `Vec`, `PetscFV`, `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()`
    4243: @*/
    4244: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
    4245: {
    4246:   PetscFunctionBegin;
    4247:   PetscCall(VecSet(vecOut, 0.0));
    4248:   if (sfRefine) {
    4249:     Vec vecInLocal;
    4250:     DM  dmGrad   = NULL;
    4251:     Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;

    4253:     PetscCall(DMGetLocalVector(dmIn, &vecInLocal));
    4254:     PetscCall(VecSet(vecInLocal, 0.0));
    4255:     {
    4256:       PetscInt numFields, i;

    4258:       PetscCall(DMGetNumFields(dmIn, &numFields));
    4259:       for (i = 0; i < numFields; i++) {
    4260:         PetscObject  obj;
    4261:         PetscClassId classid;

    4263:         PetscCall(DMGetField(dmIn, i, NULL, &obj));
    4264:         PetscCall(PetscObjectGetClassId(obj, &classid));
    4265:         if (classid == PETSCFV_CLASSID) {
    4266:           PetscCall(DMPlexGetDataFVM(dmIn, (PetscFV)obj, &cellGeom, &faceGeom, &dmGrad));
    4267:           break;
    4268:         }
    4269:       }
    4270:     }
    4271:     if (useBCs) PetscCall(DMPlexInsertBoundaryValues(dmIn, PETSC_TRUE, vecInLocal, time, faceGeom, cellGeom, NULL));
    4272:     PetscCall(DMGlobalToLocalBegin(dmIn, vecIn, INSERT_VALUES, vecInLocal));
    4273:     PetscCall(DMGlobalToLocalEnd(dmIn, vecIn, INSERT_VALUES, vecInLocal));
    4274:     if (dmGrad) {
    4275:       PetscCall(DMGetGlobalVector(dmGrad, &grad));
    4276:       PetscCall(DMPlexReconstructGradientsFVM(dmIn, vecInLocal, grad));
    4277:     }
    4278:     PetscCall(DMPlexTransferVecTree_Interpolate(dmIn, vecInLocal, dmOut, vecOut, sfRefine, cidsRefine, grad, cellGeom));
    4279:     PetscCall(DMRestoreLocalVector(dmIn, &vecInLocal));
    4280:     if (dmGrad) PetscCall(DMRestoreGlobalVector(dmGrad, &grad));
    4281:   }
    4282:   if (sfCoarsen) PetscCall(DMPlexTransferVecTree_Inject(dmIn, vecIn, dmOut, vecOut, sfCoarsen, cidsCoarsen));
    4283:   PetscCall(VecAssemblyBegin(vecOut));
    4284:   PetscCall(VecAssemblyEnd(vecOut));
    4285:   PetscFunctionReturn(PETSC_SUCCESS);
    4286: }