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