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: `DMPlexGetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`
 21: @*/
 22: PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref)
 23: {
 24:   DM_Plex *mesh = (DM_Plex *)dm->data;

 28:   PetscObjectReference((PetscObject)ref);
 29:   DMDestroy(&mesh->referenceTree);
 30:   mesh->referenceTree = ref;
 31:   return 0;
 32: }

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

 37:   Not collective

 39:   Input Parameters:
 40: . dm - The DMPlex object

 42:   Output Parameters:
 43: . ref - The reference tree DMPlex object

 45:   Level: intermediate

 47: .seealso: `DMPlexSetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`
 48: @*/
 49: PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref)
 50: {
 51:   DM_Plex *mesh = (DM_Plex *)dm->data;

 55:   *ref = mesh->referenceTree;
 56:   return 0;
 57: }

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

 63:   if (parentOrientA == parentOrientB) {
 64:     if (childOrientB) *childOrientB = childOrientA;
 65:     if (childB) *childB = childA;
 66:     return 0;
 67:   }
 68:   for (dim = 0; dim < 3; dim++) {
 69:     DMPlexGetDepthStratum(dm, dim, &dStart, &dEnd);
 70:     if (parent >= dStart && parent <= dEnd) break;
 71:   }
 74:   if (childA < dStart || childA >= dEnd) {
 75:     /* this is a lower-dimensional child: bootstrap */
 76:     PetscInt        size, i, sA = -1, sB, sOrientB, sConeSize;
 77:     const PetscInt *supp, *coneA, *coneB, *oA, *oB;

 79:     DMPlexGetSupportSize(dm, childA, &size);
 80:     DMPlexGetSupport(dm, childA, &supp);

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

 86:       sA = supp[i];
 87:       if (sA == parent) continue;
 88:       DMPlexGetTreeParent(dm, sA, &sParent, NULL);
 89:       if (sParent == parent) break;
 90:     }
 92:     /* find out which point sB is in an equivalent position to sA under
 93:      * parentOrientB */
 94:     DMPlexReferenceTreeGetChildSymmetry_Default(dm, parent, parentOrientA, 0, sA, parentOrientB, &sOrientB, &sB);
 95:     DMPlexGetConeSize(dm, sA, &sConeSize);
 96:     DMPlexGetCone(dm, sA, &coneA);
 97:     DMPlexGetCone(dm, sB, &coneB);
 98:     DMPlexGetConeOrientation(dm, sA, &oA);
 99:     DMPlexGetConeOrientation(dm, sB, &oB);
100:     /* step through the cone of sA in natural order */
101:     for (i = 0; i < sConeSize; i++) {
102:       if (coneA[i] == childA) {
103:         /* if childA is at position i in coneA,
104:          * then we want the point that is at sOrientB*i in coneB */
105:         PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize - (sOrientB + 1) - i) % sConeSize);
106:         if (childB) *childB = coneB[j];
107:         if (childOrientB) {
108:           DMPolytopeType ct;
109:           PetscInt       oBtrue;

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

145:     /* count which position the child is in */
146:     DMPlexGetTreeChildren(dm, parent, &numChildren, &children);
147:     for (i = 0; i < numChildren; i++) {
148:       p = children[i];
149:       if (p == childA) {
150:         posA = i;
151:         break;
152:       }
153:     }
154:     if (posA >= coneSize) {
155:       /* this is the triangle in the middle of a uniformly refined triangle: it
156:        * is invariant */
158:       *childB = childA;
159:     } else {
160:       /* figure out position B by applying ABswapVert */
161:       PetscInt posB;

163:       posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize - (ABswapVert + 1) - posA) % coneSize);
164:       if (childB) *childB = children[posB];
165:     }
166:   }
167:   if (childOrientB) *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
168:   return 0;
169: }

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

174:   Input Parameters:
175: + dm - the reference tree DMPlex object
176: . parent - the parent point
177: . parentOrientA - the reference orientation for describing the parent
178: . childOrientA - the reference orientation for describing the child
179: . childA - the reference childID for describing the child
180: - parentOrientB - the new orientation for describing the parent

182:   Output Parameters:
183: + childOrientB - if not NULL, set to the new oreintation for describing the child
184: - childB - if not NULL, the new childID for describing the child

186:   Level: developer

188: .seealso: `DMPlexGetReferenceTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetTree()`
189: @*/
190: PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
191: {
192:   DM_Plex *mesh = (DM_Plex *)dm->data;

196:   mesh->getchildsymmetry(dm, parent, parentOrientA, childOrientA, childA, parentOrientB, childOrientB, childB);
197:   return 0;
198: }

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

202: PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
203: {
204:   DMPlexSetTree_Internal(dm, parentSection, parents, childIDs, PETSC_TRUE, PETSC_FALSE);
205:   return 0;
206: }

208: PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref)
209: {
210:   MPI_Comm     comm;
211:   PetscInt     dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs;
212:   PetscInt    *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts;
213:   DMLabel      identity, identityRef;
214:   PetscSection unionSection, unionConeSection, parentSection;
215:   PetscScalar *unionCoords;
216:   IS           perm;

218:   comm = PetscObjectComm((PetscObject)K);
219:   DMGetDimension(K, &dim);
220:   DMPlexGetChart(K, &pStart, &pEnd);
221:   DMGetLabel(K, labelName, &identity);
222:   DMGetLabel(Kref, labelName, &identityRef);
223:   DMPlexGetChart(Kref, &pRefStart, &pRefEnd);
224:   PetscSectionCreate(comm, &unionSection);
225:   PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));
226:   /* count points that will go in the union */
227:   for (p = pStart; p < pEnd; p++) PetscSectionSetDof(unionSection, p - pStart, 1);
228:   for (p = pRefStart; p < pRefEnd; p++) {
229:     PetscInt q, qSize;
230:     DMLabelGetValue(identityRef, p, &q);
231:     DMLabelGetStratumSize(identityRef, q, &qSize);
232:     if (qSize > 1) PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);
233:   }
234:   PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart, &permvals);
235:   offset = 0;
236:   /* stratify points in the union by topological dimension */
237:   for (d = 0; d <= dim; d++) {
238:     PetscInt cStart, cEnd, c;

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

243:     DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);
244:     for (c = cStart; c < cEnd; c++) permvals[offset++] = c + (pEnd - pStart);
245:   }
246:   ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);
247:   PetscSectionSetPermutation(unionSection, perm);
248:   PetscSectionSetUp(unionSection);
249:   PetscSectionGetStorageSize(unionSection, &numUnionPoints);
250:   PetscMalloc2(numUnionPoints, &coneSizes, dim + 1, &numDimPoints);
251:   /* count dimension points */
252:   for (d = 0; d <= dim; d++) {
253:     PetscInt cStart, cOff, cOff2;
254:     DMPlexGetHeightStratum(K, d, &cStart, NULL);
255:     PetscSectionGetOffset(unionSection, cStart - pStart, &cOff);
256:     if (d < dim) {
257:       DMPlexGetHeightStratum(K, d + 1, &cStart, NULL);
258:       PetscSectionGetOffset(unionSection, cStart - pStart, &cOff2);
259:     } else {
260:       cOff2 = numUnionPoints;
261:     }
262:     numDimPoints[dim - d] = cOff2 - cOff;
263:   }
264:   PetscSectionCreate(comm, &unionConeSection);
265:   PetscSectionSetChart(unionConeSection, 0, numUnionPoints);
266:   /* count the cones in the union */
267:   for (p = pStart; p < pEnd; p++) {
268:     PetscInt dof, uOff;

270:     DMPlexGetConeSize(K, p, &dof);
271:     PetscSectionGetOffset(unionSection, p - pStart, &uOff);
272:     PetscSectionSetDof(unionConeSection, uOff, dof);
273:     coneSizes[uOff] = dof;
274:   }
275:   for (p = pRefStart; p < pRefEnd; p++) {
276:     PetscInt dof, uDof, uOff;

278:     DMPlexGetConeSize(Kref, p, &dof);
279:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof);
280:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff);
281:     if (uDof) {
282:       PetscSectionSetDof(unionConeSection, uOff, dof);
283:       coneSizes[uOff] = dof;
284:     }
285:   }
286:   PetscSectionSetUp(unionConeSection);
287:   PetscSectionGetStorageSize(unionConeSection, &numCones);
288:   PetscMalloc2(numCones, &unionCones, numCones, &unionOrientations);
289:   /* write the cones in the union */
290:   for (p = pStart; p < pEnd; p++) {
291:     PetscInt        dof, uOff, c, cOff;
292:     const PetscInt *cone, *orientation;

294:     DMPlexGetConeSize(K, p, &dof);
295:     DMPlexGetCone(K, p, &cone);
296:     DMPlexGetConeOrientation(K, p, &orientation);
297:     PetscSectionGetOffset(unionSection, p - pStart, &uOff);
298:     PetscSectionGetOffset(unionConeSection, uOff, &cOff);
299:     for (c = 0; c < dof; c++) {
300:       PetscInt e, eOff;
301:       e = cone[c];
302:       PetscSectionGetOffset(unionSection, e - pStart, &eOff);
303:       unionCones[cOff + c]        = eOff;
304:       unionOrientations[cOff + c] = orientation[c];
305:     }
306:   }
307:   for (p = pRefStart; p < pRefEnd; p++) {
308:     PetscInt        dof, uDof, uOff, c, cOff;
309:     const PetscInt *cone, *orientation;

311:     DMPlexGetConeSize(Kref, p, &dof);
312:     DMPlexGetCone(Kref, p, &cone);
313:     DMPlexGetConeOrientation(Kref, p, &orientation);
314:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof);
315:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff);
316:     if (uDof) {
317:       PetscSectionGetOffset(unionConeSection, uOff, &cOff);
318:       for (c = 0; c < dof; c++) {
319:         PetscInt e, eOff, eDof;

321:         e = cone[c];
322:         PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart), &eDof);
323:         if (eDof) {
324:           PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);
325:         } else {
326:           DMLabelGetValue(identityRef, e, &e);
327:           PetscSectionGetOffset(unionSection, e - pStart, &eOff);
328:         }
329:         unionCones[cOff + c]        = eOff;
330:         unionOrientations[cOff + c] = orientation[c];
331:       }
332:     }
333:   }
334:   /* get the coordinates */
335:   {
336:     PetscInt     vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff;
337:     PetscSection KcoordsSec, KrefCoordsSec;
338:     Vec          KcoordsVec, KrefCoordsVec;
339:     PetscScalar *Kcoords;

341:     DMGetCoordinateSection(K, &KcoordsSec);
342:     DMGetCoordinatesLocal(K, &KcoordsVec);
343:     DMGetCoordinateSection(Kref, &KrefCoordsSec);
344:     DMGetCoordinatesLocal(Kref, &KrefCoordsVec);

346:     numVerts = numDimPoints[0];
347:     PetscMalloc1(numVerts * dim, &unionCoords);
348:     DMPlexGetDepthStratum(K, 0, &vStart, &vEnd);

350:     offset = 0;
351:     for (v = vStart; v < vEnd; v++) {
352:       PetscSectionGetOffset(unionSection, v - pStart, &vOff);
353:       VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);
354:       for (d = 0; d < dim; d++) unionCoords[offset * dim + d] = Kcoords[d];
355:       offset++;
356:     }
357:     DMPlexGetDepthStratum(Kref, 0, &vRefStart, &vRefEnd);
358:     for (v = vRefStart; v < vRefEnd; v++) {
359:       PetscSectionGetDof(unionSection, v - pRefStart + (pEnd - pStart), &vDof);
360:       PetscSectionGetOffset(unionSection, v - pRefStart + (pEnd - pStart), &vOff);
361:       VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);
362:       if (vDof) {
363:         for (d = 0; d < dim; d++) unionCoords[offset * dim + d] = Kcoords[d];
364:         offset++;
365:       }
366:     }
367:   }
368:   DMCreate(comm, ref);
369:   DMSetType(*ref, DMPLEX);
370:   DMSetDimension(*ref, dim);
371:   DMPlexCreateFromDAG(*ref, dim, numDimPoints, coneSizes, unionCones, unionOrientations, unionCoords);
372:   /* set the tree */
373:   PetscSectionCreate(comm, &parentSection);
374:   PetscSectionSetChart(parentSection, 0, numUnionPoints);
375:   for (p = pRefStart; p < pRefEnd; p++) {
376:     PetscInt uDof, uOff;

378:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof);
379:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff);
380:     if (uDof) PetscSectionSetDof(parentSection, uOff, 1);
381:   }
382:   PetscSectionSetUp(parentSection);
383:   PetscSectionGetStorageSize(parentSection, &parentSize);
384:   PetscMalloc2(parentSize, &parents, parentSize, &childIDs);
385:   for (p = pRefStart; p < pRefEnd; p++) {
386:     PetscInt uDof, uOff;

388:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof);
389:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff);
390:     if (uDof) {
391:       PetscInt pOff, parent, parentU;
392:       PetscSectionGetOffset(parentSection, uOff, &pOff);
393:       DMLabelGetValue(identityRef, p, &parent);
394:       PetscSectionGetOffset(unionSection, parent - pStart, &parentU);
395:       parents[pOff]  = parentU;
396:       childIDs[pOff] = uOff;
397:     }
398:   }
399:   DMPlexCreateReferenceTree_SetTree(*ref, parentSection, parents, childIDs);
400:   PetscSectionDestroy(&parentSection);
401:   PetscFree2(parents, childIDs);

403:   /* clean up */
404:   PetscSectionDestroy(&unionSection);
405:   PetscSectionDestroy(&unionConeSection);
406:   ISDestroy(&perm);
407:   PetscFree(unionCoords);
408:   PetscFree2(unionCones, unionOrientations);
409:   PetscFree2(coneSizes, numDimPoints);
410:   return 0;
411: }

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

416:   Collective

418:   Input Parameters:
419: + comm    - the MPI communicator
420: . dim     - the spatial dimension
421: - simplex - Flag for simplex, otherwise use a tensor-product cell

423:   Output Parameters:
424: . ref     - the reference tree DMPlex object

426:   Level: intermediate

428: .seealso: `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`
429: @*/
430: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
431: {
432:   DM_Plex *mesh;
433:   DM       K, Kref;
434:   PetscInt p, pStart, pEnd;
435:   DMLabel  identity;

437: #if 1
438:   comm = PETSC_COMM_SELF;
439: #endif
440:   /* create a reference element */
441:   DMPlexCreateReferenceCell(comm, DMPolytopeTypeSimpleShape(dim, simplex), &K);
442:   DMCreateLabel(K, "identity");
443:   DMGetLabel(K, "identity", &identity);
444:   DMPlexGetChart(K, &pStart, &pEnd);
445:   for (p = pStart; p < pEnd; p++) DMLabelSetValue(identity, p, p);
446:   /* refine it */
447:   DMRefine(K, comm, &Kref);

449:   /* the reference tree is the union of these two, without duplicating
450:    * points that appear in both */
451:   DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref);
452:   mesh                   = (DM_Plex *)(*ref)->data;
453:   mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
454:   DMDestroy(&K);
455:   DMDestroy(&Kref);
456:   return 0;
457: }

459: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
460: {
461:   DM_Plex     *mesh = (DM_Plex *)dm->data;
462:   PetscSection childSec, pSec;
463:   PetscInt     p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT;
464:   PetscInt    *offsets, *children, pStart, pEnd;

467:   PetscSectionDestroy(&mesh->childSection);
468:   PetscFree(mesh->children);
469:   pSec = mesh->parentSection;
470:   if (!pSec) return 0;
471:   PetscSectionGetStorageSize(pSec, &pSize);
472:   for (p = 0; p < pSize; p++) {
473:     PetscInt par = mesh->parents[p];

475:     parMax = PetscMax(parMax, par + 1);
476:     parMin = PetscMin(parMin, par);
477:   }
478:   if (parMin > parMax) {
479:     parMin = -1;
480:     parMax = -1;
481:   }
482:   PetscSectionCreate(PetscObjectComm((PetscObject)pSec), &childSec);
483:   PetscSectionSetChart(childSec, parMin, parMax);
484:   for (p = 0; p < pSize; p++) {
485:     PetscInt par = mesh->parents[p];

487:     PetscSectionAddDof(childSec, par, 1);
488:   }
489:   PetscSectionSetUp(childSec);
490:   PetscSectionGetStorageSize(childSec, &cSize);
491:   PetscMalloc1(cSize, &children);
492:   PetscCalloc1(parMax - parMin, &offsets);
493:   PetscSectionGetChart(pSec, &pStart, &pEnd);
494:   for (p = pStart; p < pEnd; p++) {
495:     PetscInt dof, off, i;

497:     PetscSectionGetDof(pSec, p, &dof);
498:     PetscSectionGetOffset(pSec, p, &off);
499:     for (i = 0; i < dof; i++) {
500:       PetscInt par = mesh->parents[off + i], cOff;

502:       PetscSectionGetOffset(childSec, par, &cOff);
503:       children[cOff + offsets[par - parMin]++] = p;
504:     }
505:   }
506:   mesh->childSection = childSec;
507:   mesh->children     = children;
508:   PetscFree(offsets);
509:   return 0;
510: }

512: static PetscErrorCode AnchorsFlatten(PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
513: {
514:   PetscInt        pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
515:   const PetscInt *vals;
516:   PetscSection    secNew;
517:   PetscBool       anyNew, globalAnyNew;
518:   PetscBool       compress;

520:   PetscSectionGetChart(section, &pStart, &pEnd);
521:   ISGetLocalSize(is, &size);
522:   ISGetIndices(is, &vals);
523:   PetscSectionCreate(PetscObjectComm((PetscObject)section), &secNew);
524:   PetscSectionSetChart(secNew, pStart, pEnd);
525:   for (i = 0; i < size; i++) {
526:     PetscInt dof;

528:     p = vals[i];
529:     if (p < pStart || p >= pEnd) continue;
530:     PetscSectionGetDof(section, p, &dof);
531:     if (dof) break;
532:   }
533:   if (i == size) {
534:     PetscSectionSetUp(secNew);
535:     anyNew   = PETSC_FALSE;
536:     compress = PETSC_FALSE;
537:     sizeNew  = 0;
538:   } else {
539:     anyNew = PETSC_TRUE;
540:     for (p = pStart; p < pEnd; p++) {
541:       PetscInt dof, off;

543:       PetscSectionGetDof(section, p, &dof);
544:       PetscSectionGetOffset(section, p, &off);
545:       for (i = 0; i < dof; i++) {
546:         PetscInt q = vals[off + i], qDof = 0;

548:         if (q >= pStart && q < pEnd) PetscSectionGetDof(section, q, &qDof);
549:         if (qDof) PetscSectionAddDof(secNew, p, qDof);
550:         else PetscSectionAddDof(secNew, p, 1);
551:       }
552:     }
553:     PetscSectionSetUp(secNew);
554:     PetscSectionGetStorageSize(secNew, &sizeNew);
555:     PetscMalloc1(sizeNew, &valsNew);
556:     compress = PETSC_FALSE;
557:     for (p = pStart; p < pEnd; p++) {
558:       PetscInt dof, off, count, offNew, dofNew;

560:       PetscSectionGetDof(section, p, &dof);
561:       PetscSectionGetOffset(section, p, &off);
562:       PetscSectionGetDof(secNew, p, &dofNew);
563:       PetscSectionGetOffset(secNew, p, &offNew);
564:       count = 0;
565:       for (i = 0; i < dof; i++) {
566:         PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;

568:         if (q >= pStart && q < pEnd) {
569:           PetscSectionGetDof(section, q, &qDof);
570:           PetscSectionGetOffset(section, q, &qOff);
571:         }
572:         if (qDof) {
573:           PetscInt oldCount = count;

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

578:             for (k = 0; k < oldCount; k++) {
579:               if (valsNew[offNew + k] == r) break;
580:             }
581:             if (k == oldCount) valsNew[offNew + count++] = r;
582:           }
583:         } else {
584:           PetscInt k, oldCount = count;

586:           for (k = 0; k < oldCount; k++) {
587:             if (valsNew[offNew + k] == q) break;
588:           }
589:           if (k == oldCount) valsNew[offNew + count++] = q;
590:         }
591:       }
592:       if (count < dofNew) {
593:         PetscSectionSetDof(secNew, p, count);
594:         compress = PETSC_TRUE;
595:       }
596:     }
597:   }
598:   ISRestoreIndices(is, &vals);
599:   MPIU_Allreduce(&anyNew, &globalAnyNew, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)secNew));
600:   if (!globalAnyNew) {
601:     PetscSectionDestroy(&secNew);
602:     *sectionNew = NULL;
603:     *isNew      = NULL;
604:   } else {
605:     PetscBool globalCompress;

607:     MPIU_Allreduce(&compress, &globalCompress, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)secNew));
608:     if (compress) {
609:       PetscSection secComp;
610:       PetscInt    *valsComp = NULL;

612:       PetscSectionCreate(PetscObjectComm((PetscObject)section), &secComp);
613:       PetscSectionSetChart(secComp, pStart, pEnd);
614:       for (p = pStart; p < pEnd; p++) {
615:         PetscInt dof;

617:         PetscSectionGetDof(secNew, p, &dof);
618:         PetscSectionSetDof(secComp, p, dof);
619:       }
620:       PetscSectionSetUp(secComp);
621:       PetscSectionGetStorageSize(secComp, &sizeNew);
622:       PetscMalloc1(sizeNew, &valsComp);
623:       for (p = pStart; p < pEnd; p++) {
624:         PetscInt dof, off, offNew, j;

626:         PetscSectionGetDof(secNew, p, &dof);
627:         PetscSectionGetOffset(secNew, p, &off);
628:         PetscSectionGetOffset(secComp, p, &offNew);
629:         for (j = 0; j < dof; j++) valsComp[offNew + j] = valsNew[off + j];
630:       }
631:       PetscSectionDestroy(&secNew);
632:       secNew = secComp;
633:       PetscFree(valsNew);
634:       valsNew = valsComp;
635:     }
636:     ISCreateGeneral(PetscObjectComm((PetscObject)is), sizeNew, valsNew, PETSC_OWN_POINTER, isNew);
637:   }
638:   return 0;
639: }

641: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
642: {
643:   PetscInt     p, pStart, pEnd, *anchors, size;
644:   PetscInt     aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
645:   PetscSection aSec;
646:   DMLabel      canonLabel;
647:   IS           aIS;

650:   DMPlexGetChart(dm, &pStart, &pEnd);
651:   DMGetLabel(dm, "canonical", &canonLabel);
652:   for (p = pStart; p < pEnd; p++) {
653:     PetscInt parent;

655:     if (canonLabel) {
656:       PetscInt canon;

658:       DMLabelGetValue(canonLabel, p, &canon);
659:       if (p != canon) continue;
660:     }
661:     DMPlexGetTreeParent(dm, p, &parent, NULL);
662:     if (parent != p) {
663:       aMin = PetscMin(aMin, p);
664:       aMax = PetscMax(aMax, p + 1);
665:     }
666:   }
667:   if (aMin > aMax) {
668:     aMin = -1;
669:     aMax = -1;
670:   }
671:   PetscSectionCreate(PETSC_COMM_SELF, &aSec);
672:   PetscSectionSetChart(aSec, aMin, aMax);
673:   for (p = aMin; p < aMax; p++) {
674:     PetscInt parent, ancestor = p;

676:     if (canonLabel) {
677:       PetscInt canon;

679:       DMLabelGetValue(canonLabel, p, &canon);
680:       if (p != canon) continue;
681:     }
682:     DMPlexGetTreeParent(dm, p, &parent, NULL);
683:     while (parent != ancestor) {
684:       ancestor = parent;
685:       DMPlexGetTreeParent(dm, ancestor, &parent, NULL);
686:     }
687:     if (ancestor != p) {
688:       PetscInt closureSize, *closure = NULL;

690:       DMPlexGetTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure);
691:       PetscSectionSetDof(aSec, p, closureSize);
692:       DMPlexRestoreTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure);
693:     }
694:   }
695:   PetscSectionSetUp(aSec);
696:   PetscSectionGetStorageSize(aSec, &size);
697:   PetscMalloc1(size, &anchors);
698:   for (p = aMin; p < aMax; p++) {
699:     PetscInt parent, ancestor = p;

701:     if (canonLabel) {
702:       PetscInt canon;

704:       DMLabelGetValue(canonLabel, p, &canon);
705:       if (p != canon) continue;
706:     }
707:     DMPlexGetTreeParent(dm, p, &parent, NULL);
708:     while (parent != ancestor) {
709:       ancestor = parent;
710:       DMPlexGetTreeParent(dm, ancestor, &parent, NULL);
711:     }
712:     if (ancestor != p) {
713:       PetscInt j, closureSize, *closure = NULL, aOff;

715:       PetscSectionGetOffset(aSec, p, &aOff);

717:       DMPlexGetTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure);
718:       for (j = 0; j < closureSize; j++) anchors[aOff + j] = closure[2 * j];
719:       DMPlexRestoreTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure);
720:     }
721:   }
722:   ISCreateGeneral(PETSC_COMM_SELF, size, anchors, PETSC_OWN_POINTER, &aIS);
723:   {
724:     PetscSection aSecNew = aSec;
725:     IS           aISNew  = aIS;

727:     PetscObjectReference((PetscObject)aSec);
728:     PetscObjectReference((PetscObject)aIS);
729:     while (aSecNew) {
730:       PetscSectionDestroy(&aSec);
731:       ISDestroy(&aIS);
732:       aSec    = aSecNew;
733:       aIS     = aISNew;
734:       aSecNew = NULL;
735:       aISNew  = NULL;
736:       AnchorsFlatten(aSec, aIS, &aSecNew, &aISNew);
737:     }
738:   }
739:   DMPlexSetAnchors(dm, aSec, aIS);
740:   PetscSectionDestroy(&aSec);
741:   ISDestroy(&aIS);
742:   return 0;
743: }

745: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm, PetscInt p, PetscInt *dof, PetscInt *numTrueSupp)
746: {
747:   if (numTrueSupp[p] == -1) {
748:     PetscInt        i, alldof;
749:     const PetscInt *supp;
750:     PetscInt        count = 0;

752:     DMPlexGetSupportSize(dm, p, &alldof);
753:     DMPlexGetSupport(dm, p, &supp);
754:     for (i = 0; i < alldof; i++) {
755:       PetscInt        q = supp[i], numCones, j;
756:       const PetscInt *cone;

758:       DMPlexGetConeSize(dm, q, &numCones);
759:       DMPlexGetCone(dm, q, &cone);
760:       for (j = 0; j < numCones; j++) {
761:         if (cone[j] == p) break;
762:       }
763:       if (j < numCones) count++;
764:     }
765:     numTrueSupp[p] = count;
766:   }
767:   *dof = numTrueSupp[p];
768:   return 0;
769: }

771: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
772: {
773:   DM_Plex     *mesh = (DM_Plex *)dm->data;
774:   PetscSection newSupportSection;
775:   PetscInt     newSize, *newSupports, pStart, pEnd, p, d, depth;
776:   PetscInt    *numTrueSupp;
777:   PetscInt    *offsets;

780:   /* symmetrize the hierarchy */
781:   DMPlexGetDepth(dm, &depth);
782:   PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)), &newSupportSection);
783:   DMPlexGetChart(dm, &pStart, &pEnd);
784:   PetscSectionSetChart(newSupportSection, pStart, pEnd);
785:   PetscCalloc1(pEnd, &offsets);
786:   PetscMalloc1(pEnd, &numTrueSupp);
787:   for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
788:   /* if a point is in the (true) support of q, it should be in the support of
789:    * parent(q) */
790:   for (d = 0; d <= depth; d++) {
791:     DMPlexGetHeightStratum(dm, d, &pStart, &pEnd);
792:     for (p = pStart; p < pEnd; ++p) {
793:       PetscInt dof, q, qdof, parent;

795:       DMPlexGetTrueSupportSize(dm, p, &dof, numTrueSupp);
796:       PetscSectionAddDof(newSupportSection, p, dof);
797:       q = p;
798:       DMPlexGetTreeParent(dm, q, &parent, NULL);
799:       while (parent != q && parent >= pStart && parent < pEnd) {
800:         q = parent;

802:         DMPlexGetTrueSupportSize(dm, q, &qdof, numTrueSupp);
803:         PetscSectionAddDof(newSupportSection, p, qdof);
804:         PetscSectionAddDof(newSupportSection, q, dof);
805:         DMPlexGetTreeParent(dm, q, &parent, NULL);
806:       }
807:     }
808:   }
809:   PetscSectionSetUp(newSupportSection);
810:   PetscSectionGetStorageSize(newSupportSection, &newSize);
811:   PetscMalloc1(newSize, &newSupports);
812:   for (d = 0; d <= depth; d++) {
813:     DMPlexGetHeightStratum(dm, d, &pStart, &pEnd);
814:     for (p = pStart; p < pEnd; p++) {
815:       PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;

817:       PetscSectionGetDof(mesh->supportSection, p, &dof);
818:       PetscSectionGetOffset(mesh->supportSection, p, &off);
819:       PetscSectionGetDof(newSupportSection, p, &newDof);
820:       PetscSectionGetOffset(newSupportSection, p, &newOff);
821:       for (i = 0; i < dof; i++) {
822:         PetscInt        numCones, j;
823:         const PetscInt *cone;
824:         PetscInt        q = mesh->supports[off + i];

826:         DMPlexGetConeSize(dm, q, &numCones);
827:         DMPlexGetCone(dm, q, &cone);
828:         for (j = 0; j < numCones; j++) {
829:           if (cone[j] == p) break;
830:         }
831:         if (j < numCones) newSupports[newOff + offsets[p]++] = q;
832:       }

834:       q = p;
835:       DMPlexGetTreeParent(dm, q, &parent, NULL);
836:       while (parent != q && parent >= pStart && parent < pEnd) {
837:         q = parent;
838:         PetscSectionGetDof(mesh->supportSection, q, &qdof);
839:         PetscSectionGetOffset(mesh->supportSection, q, &qoff);
840:         PetscSectionGetOffset(newSupportSection, q, &newqOff);
841:         for (i = 0; i < qdof; i++) {
842:           PetscInt        numCones, j;
843:           const PetscInt *cone;
844:           PetscInt        r = mesh->supports[qoff + i];

846:           DMPlexGetConeSize(dm, r, &numCones);
847:           DMPlexGetCone(dm, r, &cone);
848:           for (j = 0; j < numCones; j++) {
849:             if (cone[j] == q) break;
850:           }
851:           if (j < numCones) newSupports[newOff + offsets[p]++] = r;
852:         }
853:         for (i = 0; i < dof; i++) {
854:           PetscInt        numCones, j;
855:           const PetscInt *cone;
856:           PetscInt        r = mesh->supports[off + i];

858:           DMPlexGetConeSize(dm, r, &numCones);
859:           DMPlexGetCone(dm, r, &cone);
860:           for (j = 0; j < numCones; j++) {
861:             if (cone[j] == p) break;
862:           }
863:           if (j < numCones) newSupports[newqOff + offsets[q]++] = r;
864:         }
865:         DMPlexGetTreeParent(dm, q, &parent, NULL);
866:       }
867:     }
868:   }
869:   PetscSectionDestroy(&mesh->supportSection);
870:   mesh->supportSection = newSupportSection;
871:   PetscFree(mesh->supports);
872:   mesh->supports = newSupports;
873:   PetscFree(offsets);
874:   PetscFree(numTrueSupp);

876:   return 0;
877: }

879: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM, PetscSection, PetscSection, Mat);
880: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM, PetscSection, PetscSection, Mat);

882: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
883: {
884:   DM_Plex *mesh = (DM_Plex *)dm->data;
885:   DM       refTree;
886:   PetscInt size;

890:   PetscObjectReference((PetscObject)parentSection);
891:   PetscSectionDestroy(&mesh->parentSection);
892:   mesh->parentSection = parentSection;
893:   PetscSectionGetStorageSize(parentSection, &size);
894:   if (parents != mesh->parents) {
895:     PetscFree(mesh->parents);
896:     PetscMalloc1(size, &mesh->parents);
897:     PetscArraycpy(mesh->parents, parents, size);
898:   }
899:   if (childIDs != mesh->childIDs) {
900:     PetscFree(mesh->childIDs);
901:     PetscMalloc1(size, &mesh->childIDs);
902:     PetscArraycpy(mesh->childIDs, childIDs, size);
903:   }
904:   DMPlexGetReferenceTree(dm, &refTree);
905:   if (refTree) {
906:     DMLabel canonLabel;

908:     DMGetLabel(refTree, "canonical", &canonLabel);
909:     if (canonLabel) {
910:       PetscInt i;

912:       for (i = 0; i < size; i++) {
913:         PetscInt canon;
914:         DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);
915:         if (canon >= 0) mesh->childIDs[i] = canon;
916:       }
917:     }
918:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
919:   } else {
920:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
921:   }
922:   DMPlexTreeSymmetrize(dm);
923:   if (computeCanonical) {
924:     PetscInt d, dim;

926:     /* add the canonical label */
927:     DMGetDimension(dm, &dim);
928:     DMCreateLabel(dm, "canonical");
929:     for (d = 0; d <= dim; d++) {
930:       PetscInt        p, dStart, dEnd, canon = -1, cNumChildren;
931:       const PetscInt *cChildren;

933:       DMPlexGetDepthStratum(dm, d, &dStart, &dEnd);
934:       for (p = dStart; p < dEnd; p++) {
935:         DMPlexGetTreeChildren(dm, p, &cNumChildren, &cChildren);
936:         if (cNumChildren) {
937:           canon = p;
938:           break;
939:         }
940:       }
941:       if (canon == -1) continue;
942:       for (p = dStart; p < dEnd; p++) {
943:         PetscInt        numChildren, i;
944:         const PetscInt *children;

946:         DMPlexGetTreeChildren(dm, p, &numChildren, &children);
947:         if (numChildren) {
949:           DMSetLabelValue(dm, "canonical", p, canon);
950:           for (i = 0; i < numChildren; i++) DMSetLabelValue(dm, "canonical", children[i], cChildren[i]);
951:         }
952:       }
953:     }
954:   }
955:   if (exchangeSupports) DMPlexTreeExchangeSupports(dm);
956:   mesh->createanchors = DMPlexCreateAnchors_Tree;
957:   /* reset anchors */
958:   DMPlexSetAnchors(dm, NULL, NULL);
959:   return 0;
960: }

962: /*@
963:   DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points.  This routine also creates
964:   the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its
965:   tree root.

967:   Collective on dm

969:   Input Parameters:
970: + dm - the DMPlex object
971: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
972:                   offset indexes the parent and childID list; the reference count of parentSection is incremented
973: . parents - a list of the point parents; copied, can be destroyed
974: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
975:              the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed

977:   Level: intermediate

979: .seealso: `DMPlexGetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()`
980: @*/
981: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
982: {
983:   DMPlexSetTree_Internal(dm, parentSection, parents, childIDs, PETSC_FALSE, PETSC_TRUE);
984:   return 0;
985: }

987: /*@
988:   DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
989:   Collective on dm

991:   Input Parameter:
992: . dm - the DMPlex object

994:   Output Parameters:
995: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
996:                   offset indexes the parent and childID list
997: . parents - a list of the point parents
998: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
999:              the child corresponds to the point in the reference tree with index childID
1000: . childSection - the inverse of the parent section
1001: - children - a list of the point children

1003:   Level: intermediate

1005: .seealso: `DMPlexSetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()`
1006: @*/
1007: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1008: {
1009:   DM_Plex *mesh = (DM_Plex *)dm->data;

1012:   if (parentSection) *parentSection = mesh->parentSection;
1013:   if (parents) *parents = mesh->parents;
1014:   if (childIDs) *childIDs = mesh->childIDs;
1015:   if (childSection) *childSection = mesh->childSection;
1016:   if (children) *children = mesh->children;
1017:   return 0;
1018: }

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

1023:   Input Parameters:
1024: + dm - the DMPlex object
1025: - point - the query point

1027:   Output Parameters:
1028: + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent
1029: - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point
1030:             does not have a parent

1032:   Level: intermediate

1034: .seealso: `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeChildren()`
1035: @*/
1036: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1037: {
1038:   DM_Plex     *mesh = (DM_Plex *)dm->data;
1039:   PetscSection pSec;

1042:   pSec = mesh->parentSection;
1043:   if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1044:     PetscInt dof;

1046:     PetscSectionGetDof(pSec, point, &dof);
1047:     if (dof) {
1048:       PetscInt off;

1050:       PetscSectionGetOffset(pSec, point, &off);
1051:       if (parent) *parent = mesh->parents[off];
1052:       if (childID) *childID = mesh->childIDs[off];
1053:       return 0;
1054:     }
1055:   }
1056:   if (parent) *parent = point;
1057:   if (childID) *childID = 0;
1058:   return 0;
1059: }

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

1064:   Input Parameters:
1065: + dm - the DMPlex object
1066: - point - the query point

1068:   Output Parameters:
1069: + numChildren - if not NULL, set to the number of children
1070: - children - if not NULL, set to a list children, or set to NULL if the point has no children

1072:   Level: intermediate

1074:   Fortran Notes:
1075:   Since it returns an array, this routine is only available in Fortran 90, and you must
1076:   include petsc.h90 in your code.

1078: .seealso: `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeParent()`
1079: @*/
1080: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1081: {
1082:   DM_Plex     *mesh = (DM_Plex *)dm->data;
1083:   PetscSection childSec;
1084:   PetscInt     dof = 0;

1087:   childSec = mesh->childSection;
1088:   if (childSec && point >= childSec->pStart && point < childSec->pEnd) PetscSectionGetDof(childSec, point, &dof);
1089:   if (numChildren) *numChildren = dof;
1090:   if (children) {
1091:     if (dof) {
1092:       PetscInt off;

1094:       PetscSectionGetOffset(childSec, point, &off);
1095:       *children = &mesh->children[off];
1096:     } else {
1097:       *children = NULL;
1098:     }
1099:   }
1100:   return 0;
1101: }

1103: 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)
1104: {
1105:   PetscInt f, b, p, c, offset, qPoints;

1107:   PetscSpaceEvaluate(space, nPoints, points, work, NULL, NULL);
1108:   for (f = 0, offset = 0; f < nFunctionals; f++) {
1109:     qPoints = pointsPerFn[f];
1110:     for (b = 0; b < nBasis; b++) {
1111:       PetscScalar val = 0.;

1113:       for (p = 0; p < qPoints; p++) {
1114:         for (c = 0; c < nComps; c++) val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c];
1115:       }
1116:       MatSetValue(basisAtPoints, b, f, val, INSERT_VALUES);
1117:     }
1118:     offset += qPoints;
1119:   }
1120:   MatAssemblyBegin(basisAtPoints, MAT_FINAL_ASSEMBLY);
1121:   MatAssemblyEnd(basisAtPoints, MAT_FINAL_ASSEMBLY);
1122:   return 0;
1123: }

1125: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
1126: {
1127:   PetscDS         ds;
1128:   PetscInt        spdim;
1129:   PetscInt        numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
1130:   const PetscInt *anchors;
1131:   PetscSection    aSec;
1132:   PetscReal      *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
1133:   IS              aIS;

1135:   DMPlexGetChart(dm, &pStart, &pEnd);
1136:   DMGetDS(dm, &ds);
1137:   PetscDSGetNumFields(ds, &numFields);
1138:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1139:   DMPlexGetAnchors(dm, &aSec, &aIS);
1140:   ISGetIndices(aIS, &anchors);
1141:   PetscSectionGetChart(cSec, &conStart, &conEnd);
1142:   DMGetDimension(dm, &spdim);
1143:   PetscMalloc6(spdim, &v0, spdim, &v0parent, spdim, &vtmp, spdim * spdim, &J, spdim * spdim, &Jparent, spdim * spdim, &invJparent);

1145:   for (f = 0; f < numFields; f++) {
1146:     PetscObject          disc;
1147:     PetscClassId         id;
1148:     PetscSpace           bspace;
1149:     PetscDualSpace       dspace;
1150:     PetscInt             i, j, k, nPoints, Nc, offset;
1151:     PetscInt             fSize, maxDof;
1152:     PetscReal           *weights, *pointsRef, *pointsReal, *work;
1153:     PetscScalar         *scwork;
1154:     const PetscScalar   *X;
1155:     PetscInt            *sizes, *workIndRow, *workIndCol;
1156:     Mat                  Amat, Bmat, Xmat;
1157:     const PetscInt      *numDof = NULL;
1158:     const PetscInt    ***perms  = NULL;
1159:     const PetscScalar ***flips  = NULL;

1161:     PetscDSGetDiscretization(ds, f, &disc);
1162:     PetscObjectGetClassId(disc, &id);
1163:     if (id == PETSCFE_CLASSID) {
1164:       PetscFE fe = (PetscFE)disc;

1166:       PetscFEGetBasisSpace(fe, &bspace);
1167:       PetscFEGetDualSpace(fe, &dspace);
1168:       PetscDualSpaceGetDimension(dspace, &fSize);
1169:       PetscFEGetNumComponents(fe, &Nc);
1170:     } else if (id == PETSCFV_CLASSID) {
1171:       PetscFV fv = (PetscFV)disc;

1173:       PetscFVGetNumComponents(fv, &Nc);
1174:       PetscSpaceCreate(PetscObjectComm((PetscObject)fv), &bspace);
1175:       PetscSpaceSetType(bspace, PETSCSPACEPOLYNOMIAL);
1176:       PetscSpaceSetDegree(bspace, 0, PETSC_DETERMINE);
1177:       PetscSpaceSetNumComponents(bspace, Nc);
1178:       PetscSpaceSetNumVariables(bspace, spdim);
1179:       PetscSpaceSetUp(bspace);
1180:       PetscFVGetDualSpace(fv, &dspace);
1181:       PetscDualSpaceGetDimension(dspace, &fSize);
1182:     } else SETERRQ(PetscObjectComm(disc), PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1183:     PetscDualSpaceGetNumDof(dspace, &numDof);
1184:     for (i = 0, maxDof = 0; i <= spdim; i++) maxDof = PetscMax(maxDof, numDof[i]);
1185:     PetscDualSpaceGetSymmetries(dspace, &perms, &flips);

1187:     MatCreate(PETSC_COMM_SELF, &Amat);
1188:     MatSetSizes(Amat, fSize, fSize, fSize, fSize);
1189:     MatSetType(Amat, MATSEQDENSE);
1190:     MatSetUp(Amat);
1191:     MatDuplicate(Amat, MAT_DO_NOT_COPY_VALUES, &Bmat);
1192:     MatDuplicate(Amat, MAT_DO_NOT_COPY_VALUES, &Xmat);
1193:     nPoints = 0;
1194:     for (i = 0; i < fSize; i++) {
1195:       PetscInt        qPoints, thisNc;
1196:       PetscQuadrature quad;

1198:       PetscDualSpaceGetFunctional(dspace, i, &quad);
1199:       PetscQuadratureGetData(quad, NULL, &thisNc, &qPoints, NULL, NULL);
1201:       nPoints += qPoints;
1202:     }
1203:     PetscMalloc7(fSize, &sizes, nPoints * Nc, &weights, spdim * nPoints, &pointsRef, spdim * nPoints, &pointsReal, nPoints * fSize * Nc, &work, maxDof, &workIndRow, maxDof, &workIndCol);
1204:     PetscMalloc1(maxDof * maxDof, &scwork);
1205:     offset = 0;
1206:     for (i = 0; i < fSize; i++) {
1207:       PetscInt         qPoints;
1208:       const PetscReal *p, *w;
1209:       PetscQuadrature  quad;

1211:       PetscDualSpaceGetFunctional(dspace, i, &quad);
1212:       PetscQuadratureGetData(quad, NULL, NULL, &qPoints, &p, &w);
1213:       PetscArraycpy(weights + Nc * offset, w, Nc * qPoints);
1214:       PetscArraycpy(pointsRef + spdim * offset, p, spdim * qPoints);
1215:       sizes[i] = qPoints;
1216:       offset += qPoints;
1217:     }
1218:     EvaluateBasis(bspace, fSize, fSize, Nc, nPoints, sizes, pointsRef, weights, work, Amat);
1219:     MatLUFactor(Amat, NULL, NULL, NULL);
1220:     for (c = cStart; c < cEnd; c++) {
1221:       PetscInt  parent;
1222:       PetscInt  closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1223:       PetscInt *childOffsets, *parentOffsets;

1225:       DMPlexGetTreeParent(dm, c, &parent, NULL);
1226:       if (parent == c) continue;
1227:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1228:       for (i = 0; i < closureSize; i++) {
1229:         PetscInt p = closure[2 * i];
1230:         PetscInt conDof;

1232:         if (p < conStart || p >= conEnd) continue;
1233:         if (numFields) {
1234:           PetscSectionGetFieldDof(cSec, p, f, &conDof);
1235:         } else {
1236:           PetscSectionGetDof(cSec, p, &conDof);
1237:         }
1238:         if (conDof) break;
1239:       }
1240:       if (i == closureSize) {
1241:         DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1242:         continue;
1243:       }

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

1250:         CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i * spdim], vtmp);
1251:         CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i * spdim]);
1252:       }
1253:       EvaluateBasis(bspace, fSize, fSize, Nc, nPoints, sizes, pointsReal, weights, work, Bmat);
1254:       MatMatSolve(Amat, Bmat, Xmat);
1255:       MatDenseGetArrayRead(Xmat, &X);
1256:       DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSizeP, &closureP);
1257:       PetscMalloc2(closureSize + 1, &childOffsets, closureSizeP + 1, &parentOffsets);
1258:       childOffsets[0] = 0;
1259:       for (i = 0; i < closureSize; i++) {
1260:         PetscInt p = closure[2 * i];
1261:         PetscInt dof;

1263:         if (numFields) {
1264:           PetscSectionGetFieldDof(section, p, f, &dof);
1265:         } else {
1266:           PetscSectionGetDof(section, p, &dof);
1267:         }
1268:         childOffsets[i + 1] = childOffsets[i] + dof;
1269:       }
1270:       parentOffsets[0] = 0;
1271:       for (i = 0; i < closureSizeP; i++) {
1272:         PetscInt p = closureP[2 * i];
1273:         PetscInt dof;

1275:         if (numFields) {
1276:           PetscSectionGetFieldDof(section, p, f, &dof);
1277:         } else {
1278:           PetscSectionGetDof(section, p, &dof);
1279:         }
1280:         parentOffsets[i + 1] = parentOffsets[i] + dof;
1281:       }
1282:       for (i = 0; i < closureSize; i++) {
1283:         PetscInt           conDof, conOff, aDof, aOff, nWork;
1284:         PetscInt           p = closure[2 * i];
1285:         PetscInt           o = closure[2 * i + 1];
1286:         const PetscInt    *perm;
1287:         const PetscScalar *flip;

1289:         if (p < conStart || p >= conEnd) continue;
1290:         if (numFields) {
1291:           PetscSectionGetFieldDof(cSec, p, f, &conDof);
1292:           PetscSectionGetFieldOffset(cSec, p, f, &conOff);
1293:         } else {
1294:           PetscSectionGetDof(cSec, p, &conDof);
1295:           PetscSectionGetOffset(cSec, p, &conOff);
1296:         }
1297:         if (!conDof) continue;
1298:         perm = (perms && perms[i]) ? perms[i][o] : NULL;
1299:         flip = (flips && flips[i]) ? flips[i][o] : NULL;
1300:         PetscSectionGetDof(aSec, p, &aDof);
1301:         PetscSectionGetOffset(aSec, p, &aOff);
1302:         nWork = childOffsets[i + 1] - childOffsets[i];
1303:         for (k = 0; k < aDof; k++) {
1304:           PetscInt a = anchors[aOff + k];
1305:           PetscInt aSecDof, aSecOff;

1307:           if (numFields) {
1308:             PetscSectionGetFieldDof(section, a, f, &aSecDof);
1309:             PetscSectionGetFieldOffset(section, a, f, &aSecOff);
1310:           } else {
1311:             PetscSectionGetDof(section, a, &aSecDof);
1312:             PetscSectionGetOffset(section, a, &aSecOff);
1313:           }
1314:           if (!aSecDof) continue;

1316:           for (j = 0; j < closureSizeP; j++) {
1317:             PetscInt q  = closureP[2 * j];
1318:             PetscInt oq = closureP[2 * j + 1];

1320:             if (q == a) {
1321:               PetscInt           r, s, nWorkP;
1322:               const PetscInt    *permP;
1323:               const PetscScalar *flipP;

1325:               permP  = (perms && perms[j]) ? perms[j][oq] : NULL;
1326:               flipP  = (flips && flips[j]) ? flips[j][oq] : NULL;
1327:               nWorkP = parentOffsets[j + 1] - parentOffsets[j];
1328:               /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the
1329:                * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArrayRead() is
1330:                * column-major, so transpose-transpose = do nothing */
1331:               for (r = 0; r < nWork; r++) {
1332:                 for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])];
1333:               }
1334:               for (r = 0; r < nWork; r++) workIndRow[perm ? perm[r] : r] = conOff + r;
1335:               for (s = 0; s < nWorkP; s++) workIndCol[permP ? permP[s] : s] = aSecOff + s;
1336:               if (flip) {
1337:                 for (r = 0; r < nWork; r++) {
1338:                   for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] *= flip[r];
1339:                 }
1340:               }
1341:               if (flipP) {
1342:                 for (r = 0; r < nWork; r++) {
1343:                   for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] *= flipP[s];
1344:                 }
1345:               }
1346:               MatSetValues(cMat, nWork, workIndRow, nWorkP, workIndCol, scwork, INSERT_VALUES);
1347:               break;
1348:             }
1349:           }
1350:         }
1351:       }
1352:       MatDenseRestoreArrayRead(Xmat, &X);
1353:       PetscFree2(childOffsets, parentOffsets);
1354:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
1355:       DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSizeP, &closureP);
1356:     }
1357:     MatDestroy(&Amat);
1358:     MatDestroy(&Bmat);
1359:     MatDestroy(&Xmat);
1360:     PetscFree(scwork);
1361:     PetscFree7(sizes, weights, pointsRef, pointsReal, work, workIndRow, workIndCol);
1362:     if (id == PETSCFV_CLASSID) PetscSpaceDestroy(&bspace);
1363:   }
1364:   MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY);
1365:   MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY);
1366:   PetscFree6(v0, v0parent, vtmp, J, Jparent, invJparent);
1367:   ISRestoreIndices(aIS, &anchors);

1369:   return 0;
1370: }

1372: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1373: {
1374:   Mat                 refCmat;
1375:   PetscDS             ds;
1376:   PetscInt            numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1377:   PetscScalar      ***refPointFieldMats;
1378:   PetscSection        refConSec, refAnSec, refSection;
1379:   IS                  refAnIS;
1380:   const PetscInt     *refAnchors;
1381:   const PetscInt    **perms;
1382:   const PetscScalar **flips;

1384:   DMGetDS(refTree, &ds);
1385:   PetscDSGetNumFields(ds, &numFields);
1386:   maxFields = PetscMax(1, numFields);
1387:   DMGetDefaultConstraints(refTree, &refConSec, &refCmat, NULL);
1388:   DMPlexGetAnchors(refTree, &refAnSec, &refAnIS);
1389:   ISGetIndices(refAnIS, &refAnchors);
1390:   DMGetLocalSection(refTree, &refSection);
1391:   PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd);
1392:   PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats);
1393:   PetscMalloc1(pRefEnd - pRefStart, &refPointFieldN);
1394:   PetscSectionGetMaxDof(refConSec, &maxDof);
1395:   PetscSectionGetMaxDof(refAnSec, &maxAnDof);
1396:   PetscMalloc1(maxDof, &rows);
1397:   PetscMalloc1(maxDof * maxAnDof, &cols);
1398:   for (p = pRefStart; p < pRefEnd; p++) {
1399:     PetscInt parent, closureSize, *closure = NULL, pDof;

1401:     DMPlexGetTreeParent(refTree, p, &parent, NULL);
1402:     PetscSectionGetDof(refConSec, p, &pDof);
1403:     if (!pDof || parent == p) continue;

1405:     PetscMalloc1(maxFields, &refPointFieldMats[p - pRefStart]);
1406:     PetscCalloc1(maxFields, &refPointFieldN[p - pRefStart]);
1407:     DMPlexGetTransitiveClosure(refTree, parent, PETSC_TRUE, &closureSize, &closure);
1408:     for (f = 0; f < maxFields; f++) {
1409:       PetscInt cDof, cOff, numCols, r, i;

1411:       if (f < numFields) {
1412:         PetscSectionGetFieldDof(refConSec, p, f, &cDof);
1413:         PetscSectionGetFieldOffset(refConSec, p, f, &cOff);
1414:         PetscSectionGetFieldPointSyms(refSection, f, closureSize, closure, &perms, &flips);
1415:       } else {
1416:         PetscSectionGetDof(refConSec, p, &cDof);
1417:         PetscSectionGetOffset(refConSec, p, &cOff);
1418:         PetscSectionGetPointSyms(refSection, closureSize, closure, &perms, &flips);
1419:       }

1421:       for (r = 0; r < cDof; r++) rows[r] = cOff + r;
1422:       numCols = 0;
1423:       for (i = 0; i < closureSize; i++) {
1424:         PetscInt        q = closure[2 * i];
1425:         PetscInt        aDof, aOff, j;
1426:         const PetscInt *perm = perms ? perms[i] : NULL;

1428:         if (numFields) {
1429:           PetscSectionGetFieldDof(refSection, q, f, &aDof);
1430:           PetscSectionGetFieldOffset(refSection, q, f, &aOff);
1431:         } else {
1432:           PetscSectionGetDof(refSection, q, &aDof);
1433:           PetscSectionGetOffset(refSection, q, &aOff);
1434:         }

1436:         for (j = 0; j < aDof; j++) cols[numCols++] = aOff + (perm ? perm[j] : j);
1437:       }
1438:       refPointFieldN[p - pRefStart][f] = numCols;
1439:       PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]);
1440:       MatGetValues(refCmat, cDof, rows, numCols, cols, refPointFieldMats[p - pRefStart][f]);
1441:       if (flips) {
1442:         PetscInt colOff = 0;

1444:         for (i = 0; i < closureSize; i++) {
1445:           PetscInt           q = closure[2 * i];
1446:           PetscInt           aDof, aOff, j;
1447:           const PetscScalar *flip = flips ? flips[i] : NULL;

1449:           if (numFields) {
1450:             PetscSectionGetFieldDof(refSection, q, f, &aDof);
1451:             PetscSectionGetFieldOffset(refSection, q, f, &aOff);
1452:           } else {
1453:             PetscSectionGetDof(refSection, q, &aDof);
1454:             PetscSectionGetOffset(refSection, q, &aOff);
1455:           }
1456:           if (flip) {
1457:             PetscInt k;
1458:             for (k = 0; k < cDof; k++) {
1459:               for (j = 0; j < aDof; j++) refPointFieldMats[p - pRefStart][f][k * numCols + colOff + j] *= flip[j];
1460:             }
1461:           }
1462:           colOff += aDof;
1463:         }
1464:       }
1465:       if (numFields) {
1466:         PetscSectionRestoreFieldPointSyms(refSection, f, closureSize, closure, &perms, &flips);
1467:       } else {
1468:         PetscSectionRestorePointSyms(refSection, closureSize, closure, &perms, &flips);
1469:       }
1470:     }
1471:     DMPlexRestoreTransitiveClosure(refTree, parent, PETSC_TRUE, &closureSize, &closure);
1472:   }
1473:   *childrenMats = refPointFieldMats;
1474:   *childrenN    = refPointFieldN;
1475:   ISRestoreIndices(refAnIS, &refAnchors);
1476:   PetscFree(rows);
1477:   PetscFree(cols);
1478:   return 0;
1479: }

1481: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1482: {
1483:   PetscDS        ds;
1484:   PetscInt     **refPointFieldN;
1485:   PetscScalar ***refPointFieldMats;
1486:   PetscInt       numFields, maxFields, pRefStart, pRefEnd, p, f;
1487:   PetscSection   refConSec;

1489:   refPointFieldN    = *childrenN;
1490:   *childrenN        = NULL;
1491:   refPointFieldMats = *childrenMats;
1492:   *childrenMats     = NULL;
1493:   DMGetDS(refTree, &ds);
1494:   PetscDSGetNumFields(ds, &numFields);
1495:   maxFields = PetscMax(1, numFields);
1496:   DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL);
1497:   PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd);
1498:   for (p = pRefStart; p < pRefEnd; p++) {
1499:     PetscInt parent, pDof;

1501:     DMPlexGetTreeParent(refTree, p, &parent, NULL);
1502:     PetscSectionGetDof(refConSec, p, &pDof);
1503:     if (!pDof || parent == p) continue;

1505:     for (f = 0; f < maxFields; f++) {
1506:       PetscInt cDof;

1508:       if (numFields) {
1509:         PetscSectionGetFieldDof(refConSec, p, f, &cDof);
1510:       } else {
1511:         PetscSectionGetDof(refConSec, p, &cDof);
1512:       }

1514:       PetscFree(refPointFieldMats[p - pRefStart][f]);
1515:     }
1516:     PetscFree(refPointFieldMats[p - pRefStart]);
1517:     PetscFree(refPointFieldN[p - pRefStart]);
1518:   }
1519:   PetscFree(refPointFieldMats);
1520:   PetscFree(refPointFieldN);
1521:   return 0;
1522: }

1524: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1525: {
1526:   DM              refTree;
1527:   PetscDS         ds;
1528:   Mat             refCmat;
1529:   PetscInt        numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1530:   PetscScalar  ***refPointFieldMats, *pointWork;
1531:   PetscSection    refConSec, refAnSec, anSec;
1532:   IS              refAnIS, anIS;
1533:   const PetscInt *anchors;

1536:   DMGetDS(dm, &ds);
1537:   PetscDSGetNumFields(ds, &numFields);
1538:   maxFields = PetscMax(1, numFields);
1539:   DMPlexGetReferenceTree(dm, &refTree);
1540:   DMCopyDisc(dm, refTree);
1541:   DMGetDefaultConstraints(refTree, &refConSec, &refCmat, NULL);
1542:   DMPlexGetAnchors(refTree, &refAnSec, &refAnIS);
1543:   DMPlexGetAnchors(dm, &anSec, &anIS);
1544:   ISGetIndices(anIS, &anchors);
1545:   PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd);
1546:   PetscSectionGetChart(conSec, &conStart, &conEnd);
1547:   PetscSectionGetMaxDof(refConSec, &maxDof);
1548:   PetscSectionGetMaxDof(refAnSec, &maxAnDof);
1549:   PetscMalloc1(maxDof * maxDof * maxAnDof, &pointWork);

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

1554:   /* step 2: compute the preorder */
1555:   DMPlexGetChart(dm, &pStart, &pEnd);
1556:   PetscMalloc2(pEnd - pStart, &perm, pEnd - pStart, &iperm);
1557:   for (p = pStart; p < pEnd; p++) {
1558:     perm[p - pStart]  = p;
1559:     iperm[p - pStart] = p - pStart;
1560:   }
1561:   for (p = 0; p < pEnd - pStart;) {
1562:     PetscInt point = perm[p];
1563:     PetscInt parent;

1565:     DMPlexGetTreeParent(dm, point, &parent, NULL);
1566:     if (parent == point) {
1567:       p++;
1568:     } else {
1569:       PetscInt size, closureSize, *closure = NULL, i;

1571:       DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
1572:       for (i = 0; i < closureSize; i++) {
1573:         PetscInt q = closure[2 * i];
1574:         if (iperm[q - pStart] > iperm[point - pStart]) {
1575:           /* swap */
1576:           perm[p]                 = q;
1577:           perm[iperm[q - pStart]] = point;
1578:           iperm[point - pStart]   = iperm[q - pStart];
1579:           iperm[q - pStart]       = p;
1580:           break;
1581:         }
1582:       }
1583:       size = closureSize;
1584:       DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
1585:       if (i == size) p++;
1586:     }
1587:   }

1589:   /* step 3: fill the constraint matrix */
1590:   /* we are going to use a preorder progressive fill strategy.  Mat doesn't
1591:    * allow progressive fill without assembly, so we are going to set up the
1592:    * values outside of the Mat first.
1593:    */
1594:   {
1595:     PetscInt        nRows, row, nnz;
1596:     PetscBool       done;
1597:     PetscInt        secStart, secEnd;
1598:     const PetscInt *ia, *ja;
1599:     PetscScalar    *vals;

1601:     PetscSectionGetChart(section, &secStart, &secEnd);
1602:     MatGetRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done);
1604:     nnz = ia[nRows];
1605:     /* malloc and then zero rows right before we fill them: this way valgrind
1606:      * can tell if we are doing progressive fill in the wrong order */
1607:     PetscMalloc1(nnz, &vals);
1608:     for (p = 0; p < pEnd - pStart; p++) {
1609:       PetscInt parent, childid, closureSize, *closure = NULL;
1610:       PetscInt point = perm[p], pointDof;

1612:       DMPlexGetTreeParent(dm, point, &parent, &childid);
1613:       if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1614:       PetscSectionGetDof(conSec, point, &pointDof);
1615:       if (!pointDof) continue;
1616:       DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
1617:       for (f = 0; f < maxFields; f++) {
1618:         PetscInt            cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
1619:         PetscScalar        *pointMat;
1620:         const PetscInt    **perms;
1621:         const PetscScalar **flips;

1623:         if (numFields) {
1624:           PetscSectionGetFieldDof(conSec, point, f, &cDof);
1625:           PetscSectionGetFieldOffset(conSec, point, f, &cOff);
1626:         } else {
1627:           PetscSectionGetDof(conSec, point, &cDof);
1628:           PetscSectionGetOffset(conSec, point, &cOff);
1629:         }
1630:         if (!cDof) continue;
1631:         if (numFields) PetscSectionGetFieldPointSyms(section, f, closureSize, closure, &perms, &flips);
1632:         else PetscSectionGetPointSyms(section, closureSize, closure, &perms, &flips);

1634:         /* make sure that every row for this point is the same size */
1635:         if (PetscDefined(USE_DEBUG)) {
1636:           for (r = 0; r < cDof; r++) {
1637:             if (cDof > 1 && r) {
1639:             }
1640:           }
1641:         }
1642:         /* zero rows */
1643:         for (i = ia[cOff]; i < ia[cOff + cDof]; i++) vals[i] = 0.;
1644:         matOffset   = ia[cOff];
1645:         numFillCols = ia[cOff + 1] - matOffset;
1646:         pointMat    = refPointFieldMats[childid - pRefStart][f];
1647:         numCols     = refPointFieldN[childid - pRefStart][f];
1648:         offset      = 0;
1649:         for (i = 0; i < closureSize; i++) {
1650:           PetscInt           q = closure[2 * i];
1651:           PetscInt           aDof, aOff, j, k, qConDof, qConOff;
1652:           const PetscInt    *perm = perms ? perms[i] : NULL;
1653:           const PetscScalar *flip = flips ? flips[i] : NULL;

1655:           qConDof = qConOff = 0;
1656:           if (q < secStart || q >= secEnd) continue;
1657:           if (numFields) {
1658:             PetscSectionGetFieldDof(section, q, f, &aDof);
1659:             PetscSectionGetFieldOffset(section, q, f, &aOff);
1660:             if (q >= conStart && q < conEnd) {
1661:               PetscSectionGetFieldDof(conSec, q, f, &qConDof);
1662:               PetscSectionGetFieldOffset(conSec, q, f, &qConOff);
1663:             }
1664:           } else {
1665:             PetscSectionGetDof(section, q, &aDof);
1666:             PetscSectionGetOffset(section, q, &aOff);
1667:             if (q >= conStart && q < conEnd) {
1668:               PetscSectionGetDof(conSec, q, &qConDof);
1669:               PetscSectionGetOffset(conSec, q, &qConOff);
1670:             }
1671:           }
1672:           if (!aDof) continue;
1673:           if (qConDof) {
1674:             /* this point has anchors: its rows of the matrix should already
1675:              * be filled, thanks to preordering */
1676:             /* first multiply into pointWork, then set in matrix */
1677:             PetscInt aMatOffset   = ia[qConOff];
1678:             PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1679:             for (r = 0; r < cDof; r++) {
1680:               for (j = 0; j < aNumFillCols; j++) {
1681:                 PetscScalar inVal = 0;
1682:                 for (k = 0; k < aDof; k++) {
1683:                   PetscInt col = perm ? perm[k] : k;

1685:                   inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.);
1686:                 }
1687:                 pointWork[r * aNumFillCols + j] = inVal;
1688:               }
1689:             }
1690:             /* assume that the columns are sorted, spend less time searching */
1691:             for (j = 0, k = 0; j < aNumFillCols; j++) {
1692:               PetscInt col = ja[aMatOffset + j];
1693:               for (; k < numFillCols; k++) {
1694:                 if (ja[matOffset + k] == col) break;
1695:               }
1697:               for (r = 0; r < cDof; r++) vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1698:             }
1699:           } else {
1700:             /* find where to put this portion of pointMat into the matrix */
1701:             for (k = 0; k < numFillCols; k++) {
1702:               if (ja[matOffset + k] == aOff) break;
1703:             }
1705:             for (r = 0; r < cDof; r++) {
1706:               for (j = 0; j < aDof; j++) {
1707:                 PetscInt col = perm ? perm[j] : j;

1709:                 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
1710:               }
1711:             }
1712:           }
1713:           offset += aDof;
1714:         }
1715:         if (numFields) {
1716:           PetscSectionRestoreFieldPointSyms(section, f, closureSize, closure, &perms, &flips);
1717:         } else {
1718:           PetscSectionRestorePointSyms(section, closureSize, closure, &perms, &flips);
1719:         }
1720:       }
1721:       DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
1722:     }
1723:     for (row = 0; row < nRows; row++) MatSetValues(cMat, 1, &row, ia[row + 1] - ia[row], &ja[ia[row]], &vals[ia[row]], INSERT_VALUES);
1724:     MatRestoreRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done);
1726:     MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY);
1727:     MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY);
1728:     PetscFree(vals);
1729:   }

1731:   /* clean up */
1732:   ISRestoreIndices(anIS, &anchors);
1733:   PetscFree2(perm, iperm);
1734:   PetscFree(pointWork);
1735:   DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN);
1736:   return 0;
1737: }

1739: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1740:  * a non-conforming mesh.  Local refinement comes later */
1741: PetscErrorCode DMPlexTreeRefineCell(DM dm, PetscInt cell, DM *ncdm)
1742: {
1743:   DM           K;
1744:   PetscMPIInt  rank;
1745:   PetscInt     dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1746:   PetscInt     numNewCones, *newConeSizes, *newCones, *newOrientations;
1747:   PetscInt    *Kembedding;
1748:   PetscInt    *cellClosure = NULL, nc;
1749:   PetscScalar *newVertexCoords;
1750:   PetscInt     numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1751:   PetscSection parentSection;

1753:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1754:   DMGetDimension(dm, &dim);
1755:   DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);
1756:   DMSetDimension(*ncdm, dim);

1758:   DMPlexGetChart(dm, &pStart, &pEnd);
1759:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &parentSection);
1760:   DMPlexGetReferenceTree(dm, &K);
1761:   DMGetCoordinatesLocalSetUp(dm);
1762:   if (rank == 0) {
1763:     /* compute the new charts */
1764:     PetscMalloc5(dim + 1, &pNewCount, dim + 1, &pNewStart, dim + 1, &pNewEnd, dim + 1, &pOldStart, dim + 1, &pOldEnd);
1765:     offset = 0;
1766:     for (d = 0; d <= dim; d++) {
1767:       PetscInt pOldCount, kStart, kEnd, k;

1769:       pNewStart[d] = offset;
1770:       DMPlexGetHeightStratum(dm, d, &pOldStart[d], &pOldEnd[d]);
1771:       DMPlexGetHeightStratum(K, d, &kStart, &kEnd);
1772:       pOldCount = pOldEnd[d] - pOldStart[d];
1773:       /* adding the new points */
1774:       pNewCount[d] = pOldCount + kEnd - kStart;
1775:       if (!d) {
1776:         /* removing the cell */
1777:         pNewCount[d]--;
1778:       }
1779:       for (k = kStart; k < kEnd; k++) {
1780:         PetscInt parent;
1781:         DMPlexGetTreeParent(K, k, &parent, NULL);
1782:         if (parent == k) {
1783:           /* avoid double counting points that won't actually be new */
1784:           pNewCount[d]--;
1785:         }
1786:       }
1787:       pNewEnd[d] = pNewStart[d] + pNewCount[d];
1788:       offset     = pNewEnd[d];
1789:     }
1791:     /* get the current closure of the cell that we are removing */
1792:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure);

1794:     PetscMalloc1(pNewEnd[dim], &newConeSizes);
1795:     {
1796:       DMPolytopeType pct, qct;
1797:       PetscInt       kStart, kEnd, k, closureSizeK, *closureK = NULL, j;

1799:       DMPlexGetChart(K, &kStart, &kEnd);
1800:       PetscMalloc4(kEnd - kStart, &Kembedding, kEnd - kStart, &perm, kEnd - kStart, &iperm, kEnd - kStart, &preOrient);

1802:       for (k = kStart; k < kEnd; k++) {
1803:         perm[k - kStart]      = k;
1804:         iperm[k - kStart]     = k - kStart;
1805:         preOrient[k - kStart] = 0;
1806:       }

1808:       DMPlexGetTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK);
1809:       for (j = 1; j < closureSizeK; j++) {
1810:         PetscInt parentOrientA = closureK[2 * j + 1];
1811:         PetscInt parentOrientB = cellClosure[2 * j + 1];
1812:         PetscInt p, q;

1814:         p = closureK[2 * j];
1815:         q = cellClosure[2 * j];
1816:         DMPlexGetCellType(K, p, &pct);
1817:         DMPlexGetCellType(dm, q, &qct);
1818:         for (d = 0; d <= dim; d++) {
1819:           if (q >= pOldStart[d] && q < pOldEnd[d]) Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1820:         }
1821:         parentOrientA = DMPolytopeConvertNewOrientation_Internal(pct, parentOrientA);
1822:         parentOrientB = DMPolytopeConvertNewOrientation_Internal(qct, parentOrientB);
1823:         if (parentOrientA != parentOrientB) {
1824:           PetscInt        numChildren, i;
1825:           const PetscInt *children;

1827:           DMPlexGetTreeChildren(K, p, &numChildren, &children);
1828:           for (i = 0; i < numChildren; i++) {
1829:             PetscInt kPerm, oPerm;

1831:             k = children[i];
1832:             DMPlexReferenceTreeGetChildSymmetry(K, p, parentOrientA, 0, k, parentOrientB, &oPerm, &kPerm);
1833:             /* perm = what refTree position I'm in */
1834:             perm[kPerm - kStart] = k;
1835:             /* iperm = who is at this position */
1836:             iperm[k - kStart]         = kPerm - kStart;
1837:             preOrient[kPerm - kStart] = oPerm;
1838:           }
1839:         }
1840:       }
1841:       DMPlexRestoreTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK);
1842:     }
1843:     PetscSectionSetChart(parentSection, 0, pNewEnd[dim]);
1844:     offset      = 0;
1845:     numNewCones = 0;
1846:     for (d = 0; d <= dim; d++) {
1847:       PetscInt kStart, kEnd, k;
1848:       PetscInt p;
1849:       PetscInt size;

1851:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1852:         /* skip cell 0 */
1853:         if (p == cell) continue;
1854:         /* old cones to new cones */
1855:         DMPlexGetConeSize(dm, p, &size);
1856:         newConeSizes[offset++] = size;
1857:         numNewCones += size;
1858:       }

1860:       DMPlexGetHeightStratum(K, d, &kStart, &kEnd);
1861:       for (k = kStart; k < kEnd; k++) {
1862:         PetscInt kParent;

1864:         DMPlexGetTreeParent(K, k, &kParent, NULL);
1865:         if (kParent != k) {
1866:           Kembedding[k] = offset;
1867:           DMPlexGetConeSize(K, k, &size);
1868:           newConeSizes[offset++] = size;
1869:           numNewCones += size;
1870:           if (kParent != 0) PetscSectionSetDof(parentSection, Kembedding[k], 1);
1871:         }
1872:       }
1873:     }

1875:     PetscSectionSetUp(parentSection);
1876:     PetscSectionGetStorageSize(parentSection, &numPointsWithParents);
1877:     PetscMalloc2(numNewCones, &newCones, numNewCones, &newOrientations);
1878:     PetscMalloc2(numPointsWithParents, &parents, numPointsWithParents, &childIDs);

1880:     /* fill new cones */
1881:     offset = 0;
1882:     for (d = 0; d <= dim; d++) {
1883:       PetscInt        kStart, kEnd, k, l;
1884:       PetscInt        p;
1885:       PetscInt        size;
1886:       const PetscInt *cone, *orientation;

1888:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1889:         /* skip cell 0 */
1890:         if (p == cell) continue;
1891:         /* old cones to new cones */
1892:         DMPlexGetConeSize(dm, p, &size);
1893:         DMPlexGetCone(dm, p, &cone);
1894:         DMPlexGetConeOrientation(dm, p, &orientation);
1895:         for (l = 0; l < size; l++) {
1896:           newCones[offset]          = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
1897:           newOrientations[offset++] = orientation[l];
1898:         }
1899:       }

1901:       DMPlexGetHeightStratum(K, d, &kStart, &kEnd);
1902:       for (k = kStart; k < kEnd; k++) {
1903:         PetscInt kPerm = perm[k], kParent;
1904:         PetscInt preO  = preOrient[k];

1906:         DMPlexGetTreeParent(K, k, &kParent, NULL);
1907:         if (kParent != k) {
1908:           /* embed new cones */
1909:           DMPlexGetConeSize(K, k, &size);
1910:           DMPlexGetCone(K, kPerm, &cone);
1911:           DMPlexGetConeOrientation(K, kPerm, &orientation);
1912:           for (l = 0; l < size; l++) {
1913:             PetscInt       q, m = (preO >= 0) ? ((preO + l) % size) : ((size - (preO + 1) - l) % size);
1914:             PetscInt       newO, lSize, oTrue;
1915:             DMPolytopeType ct = DM_NUM_POLYTOPES;

1917:             q                = iperm[cone[m]];
1918:             newCones[offset] = Kembedding[q];
1919:             DMPlexGetConeSize(K, q, &lSize);
1920:             if (lSize == 2) ct = DM_POLYTOPE_SEGMENT;
1921:             else if (lSize == 4) ct = DM_POLYTOPE_QUADRILATERAL;
1922:             oTrue                     = DMPolytopeConvertNewOrientation_Internal(ct, orientation[m]);
1923:             oTrue                     = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
1924:             newO                      = DihedralCompose(lSize, oTrue, preOrient[q]);
1925:             newOrientations[offset++] = DMPolytopeConvertOldOrientation_Internal(ct, newO);
1926:           }
1927:           if (kParent != 0) {
1928:             PetscInt newPoint = Kembedding[kParent];
1929:             PetscSectionGetOffset(parentSection, Kembedding[k], &pOffset);
1930:             parents[pOffset]  = newPoint;
1931:             childIDs[pOffset] = k;
1932:           }
1933:         }
1934:       }
1935:     }

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

1939:     /* fill coordinates */
1940:     offset = 0;
1941:     {
1942:       PetscInt     kStart, kEnd, l;
1943:       PetscSection vSection;
1944:       PetscInt     v;
1945:       Vec          coords;
1946:       PetscScalar *coordvals;
1947:       PetscInt     dof, off;
1948:       PetscReal    v0[3], J[9], detJ;

1950:       if (PetscDefined(USE_DEBUG)) {
1951:         PetscInt k;
1952:         DMPlexGetHeightStratum(K, 0, &kStart, &kEnd);
1953:         for (k = kStart; k < kEnd; k++) {
1954:           DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);
1956:         }
1957:       }
1958:       DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);
1959:       DMGetCoordinateSection(dm, &vSection);
1960:       DMGetCoordinatesLocal(dm, &coords);
1961:       VecGetArray(coords, &coordvals);
1962:       for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
1963:         PetscSectionGetDof(vSection, v, &dof);
1964:         PetscSectionGetOffset(vSection, v, &off);
1965:         for (l = 0; l < dof; l++) newVertexCoords[offset++] = coordvals[off + l];
1966:       }
1967:       VecRestoreArray(coords, &coordvals);

1969:       DMGetCoordinateSection(K, &vSection);
1970:       DMGetCoordinatesLocal(K, &coords);
1971:       VecGetArray(coords, &coordvals);
1972:       DMPlexGetDepthStratum(K, 0, &kStart, &kEnd);
1973:       for (v = kStart; v < kEnd; v++) {
1974:         PetscReal       coord[3], newCoord[3];
1975:         PetscInt        vPerm = perm[v];
1976:         PetscInt        kParent;
1977:         const PetscReal xi0[3] = {-1., -1., -1.};

1979:         DMPlexGetTreeParent(K, v, &kParent, NULL);
1980:         if (kParent != v) {
1981:           /* this is a new vertex */
1982:           PetscSectionGetOffset(vSection, vPerm, &off);
1983:           for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off + l]);
1984:           CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
1985:           for (l = 0; l < dim; ++l) newVertexCoords[offset + l] = newCoord[l];
1986:           offset += dim;
1987:         }
1988:       }
1989:       VecRestoreArray(coords, &coordvals);
1990:     }

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

1996:       tmp                = pNewCount[d];
1997:       pNewCount[d]       = pNewCount[dim - d];
1998:       pNewCount[dim - d] = tmp;
1999:     }

2001:     DMPlexCreateFromDAG(*ncdm, dim, pNewCount, newConeSizes, newCones, newOrientations, newVertexCoords);
2002:     DMPlexSetReferenceTree(*ncdm, K);
2003:     DMPlexSetTree(*ncdm, parentSection, parents, childIDs);

2005:     /* clean up */
2006:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure);
2007:     PetscFree5(pNewCount, pNewStart, pNewEnd, pOldStart, pOldEnd);
2008:     PetscFree(newConeSizes);
2009:     PetscFree2(newCones, newOrientations);
2010:     PetscFree(newVertexCoords);
2011:     PetscFree2(parents, childIDs);
2012:     PetscFree4(Kembedding, perm, iperm, preOrient);
2013:   } else {
2014:     PetscInt     p, counts[4];
2015:     PetscInt    *coneSizes, *cones, *orientations;
2016:     Vec          coordVec;
2017:     PetscScalar *coords;

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

2022:       DMPlexGetDepthStratum(dm, d, &dStart, &dEnd);
2023:       counts[d] = dEnd - dStart;
2024:     }
2025:     PetscMalloc1(pEnd - pStart, &coneSizes);
2026:     for (p = pStart; p < pEnd; p++) DMPlexGetConeSize(dm, p, &coneSizes[p - pStart]);
2027:     DMPlexGetCones(dm, &cones);
2028:     DMPlexGetConeOrientations(dm, &orientations);
2029:     DMGetCoordinatesLocal(dm, &coordVec);
2030:     VecGetArray(coordVec, &coords);

2032:     PetscSectionSetChart(parentSection, pStart, pEnd);
2033:     PetscSectionSetUp(parentSection);
2034:     DMPlexCreateFromDAG(*ncdm, dim, counts, coneSizes, cones, orientations, NULL);
2035:     DMPlexSetReferenceTree(*ncdm, K);
2036:     DMPlexSetTree(*ncdm, parentSection, NULL, NULL);
2037:     VecRestoreArray(coordVec, &coords);
2038:   }
2039:   PetscSectionDestroy(&parentSection);

2041:   return 0;
2042: }

2044: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2045: {
2046:   PetscSF              coarseToFineEmbedded;
2047:   PetscSection         globalCoarse, globalFine;
2048:   PetscSection         localCoarse, localFine;
2049:   PetscSection         aSec, cSec;
2050:   PetscSection         rootIndicesSec, rootMatricesSec;
2051:   PetscSection         leafIndicesSec, leafMatricesSec;
2052:   PetscInt            *rootIndices, *leafIndices;
2053:   PetscScalar         *rootMatrices, *leafMatrices;
2054:   IS                   aIS;
2055:   const PetscInt      *anchors;
2056:   Mat                  cMat;
2057:   PetscInt             numFields, maxFields;
2058:   PetscInt             pStartC, pEndC, pStartF, pEndF, p;
2059:   PetscInt             aStart, aEnd, cStart, cEnd;
2060:   PetscInt            *maxChildIds;
2061:   PetscInt            *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2062:   const PetscInt    ***perms;
2063:   const PetscScalar ***flips;

2065:   DMPlexGetChart(coarse, &pStartC, &pEndC);
2066:   DMPlexGetChart(fine, &pStartF, &pEndF);
2067:   DMGetGlobalSection(fine, &globalFine);
2068:   { /* winnow fine points that don't have global dofs out of the sf */
2069:     PetscInt        dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2070:     const PetscInt *leaves;

2072:     PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL);
2073:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2074:       p = leaves ? leaves[l] : l;
2075:       PetscSectionGetDof(globalFine, p, &dof);
2076:       PetscSectionGetConstraintDof(globalFine, p, &cdof);
2077:       if ((dof - cdof) > 0) numPointsWithDofs++;
2078:     }
2079:     PetscMalloc1(numPointsWithDofs, &pointsWithDofs);
2080:     for (l = 0, offset = 0; l < nleaves; l++) {
2081:       p = leaves ? leaves[l] : l;
2082:       PetscSectionGetDof(globalFine, p, &dof);
2083:       PetscSectionGetConstraintDof(globalFine, p, &cdof);
2084:       if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
2085:     }
2086:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
2087:     PetscFree(pointsWithDofs);
2088:   }
2089:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2090:   PetscMalloc1(pEndC - pStartC, &maxChildIds);
2091:   for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
2092:   PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX);
2093:   PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX);

2095:   DMGetLocalSection(coarse, &localCoarse);
2096:   DMGetGlobalSection(coarse, &globalCoarse);

2098:   DMPlexGetAnchors(coarse, &aSec, &aIS);
2099:   ISGetIndices(aIS, &anchors);
2100:   PetscSectionGetChart(aSec, &aStart, &aEnd);

2102:   DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL);
2103:   PetscSectionGetChart(cSec, &cStart, &cEnd);

2105:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2106:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec);
2107:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootMatricesSec);
2108:   PetscSectionSetChart(rootIndicesSec, pStartC, pEndC);
2109:   PetscSectionSetChart(rootMatricesSec, pStartC, pEndC);
2110:   PetscSectionGetNumFields(localCoarse, &numFields);
2111:   maxFields = PetscMax(1, numFields);
2112:   PetscMalloc7(maxFields + 1, &offsets, maxFields + 1, &offsetsCopy, maxFields + 1, &newOffsets, maxFields + 1, &newOffsetsCopy, maxFields + 1, &rowOffsets, maxFields + 1, &numD, maxFields + 1, &numO);
2113:   PetscMalloc2(maxFields + 1, (PetscInt ****)&perms, maxFields + 1, (PetscScalar ****)&flips);
2114:   PetscMemzero((void *)perms, (maxFields + 1) * sizeof(const PetscInt **));
2115:   PetscMemzero((void *)flips, (maxFields + 1) * sizeof(const PetscScalar **));

2117:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2118:     PetscInt dof, matSize = 0;
2119:     PetscInt aDof          = 0;
2120:     PetscInt cDof          = 0;
2121:     PetscInt maxChildId    = maxChildIds[p - pStartC];
2122:     PetscInt numRowIndices = 0;
2123:     PetscInt numColIndices = 0;
2124:     PetscInt f;

2126:     PetscSectionGetDof(globalCoarse, p, &dof);
2127:     if (dof < 0) dof = -(dof + 1);
2128:     if (p >= aStart && p < aEnd) PetscSectionGetDof(aSec, p, &aDof);
2129:     if (p >= cStart && p < cEnd) PetscSectionGetDof(cSec, p, &cDof);
2130:     for (f = 0; f <= numFields; f++) offsets[f] = 0;
2131:     for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2132:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2133:       PetscInt *closure = NULL, closureSize, cl;

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

2139:         PetscSectionGetDof(localCoarse, c, &clDof);
2140:         numRowIndices += clDof;
2141:         for (f = 0; f < numFields; f++) {
2142:           PetscSectionGetFieldDof(localCoarse, c, f, &clDof);
2143:           offsets[f + 1] += clDof;
2144:         }
2145:       }
2146:       for (f = 0; f < numFields; f++) {
2147:         offsets[f + 1] += offsets[f];
2148:         newOffsets[f + 1] = offsets[f + 1];
2149:       }
2150:       /* get the number of indices needed and their field offsets */
2151:       DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, NULL, NULL, NULL, &numColIndices, NULL, NULL, newOffsets, PETSC_FALSE);
2152:       DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2153:       if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2154:         numColIndices = numRowIndices;
2155:         matSize       = 0;
2156:       } else if (numFields) { /* we send one submat for each field: sum their sizes */
2157:         matSize = 0;
2158:         for (f = 0; f < numFields; f++) {
2159:           PetscInt numRow, numCol;

2161:           numRow = offsets[f + 1] - offsets[f];
2162:           numCol = newOffsets[f + 1] - newOffsets[f];
2163:           matSize += numRow * numCol;
2164:         }
2165:       } else {
2166:         matSize = numRowIndices * numColIndices;
2167:       }
2168:     } else if (maxChildId == -1) {
2169:       if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2170:         PetscInt aOff, a;

2172:         PetscSectionGetOffset(aSec, p, &aOff);
2173:         for (f = 0; f < numFields; f++) {
2174:           PetscInt fDof;

2176:           PetscSectionGetFieldDof(localCoarse, p, f, &fDof);
2177:           offsets[f + 1] = fDof;
2178:         }
2179:         for (a = 0; a < aDof; a++) {
2180:           PetscInt anchor = anchors[a + aOff], aLocalDof;

2182:           PetscSectionGetDof(localCoarse, anchor, &aLocalDof);
2183:           numColIndices += aLocalDof;
2184:           for (f = 0; f < numFields; f++) {
2185:             PetscInt fDof;

2187:             PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof);
2188:             newOffsets[f + 1] += fDof;
2189:           }
2190:         }
2191:         if (numFields) {
2192:           matSize = 0;
2193:           for (f = 0; f < numFields; f++) matSize += offsets[f + 1] * newOffsets[f + 1];
2194:         } else {
2195:           matSize = numColIndices * dof;
2196:         }
2197:       } else { /* no children, and no constraints on dofs: just get the global indices */
2198:         numColIndices = dof;
2199:         matSize       = 0;
2200:       }
2201:     }
2202:     /* we will pack the column indices with the field offsets */
2203:     PetscSectionSetDof(rootIndicesSec, p, numColIndices ? numColIndices + 2 * numFields : 0);
2204:     PetscSectionSetDof(rootMatricesSec, p, matSize);
2205:   }
2206:   PetscSectionSetUp(rootIndicesSec);
2207:   PetscSectionSetUp(rootMatricesSec);
2208:   {
2209:     PetscInt numRootIndices, numRootMatrices;

2211:     PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices);
2212:     PetscSectionGetStorageSize(rootMatricesSec, &numRootMatrices);
2213:     PetscMalloc2(numRootIndices, &rootIndices, numRootMatrices, &rootMatrices);
2214:     for (p = pStartC; p < pEndC; p++) {
2215:       PetscInt     numRowIndices, numColIndices, matSize, dof;
2216:       PetscInt     pIndOff, pMatOff, f;
2217:       PetscInt    *pInd;
2218:       PetscInt     maxChildId = maxChildIds[p - pStartC];
2219:       PetscScalar *pMat       = NULL;

2221:       PetscSectionGetDof(rootIndicesSec, p, &numColIndices);
2222:       if (!numColIndices) continue;
2223:       for (f = 0; f <= numFields; f++) {
2224:         offsets[f]        = 0;
2225:         newOffsets[f]     = 0;
2226:         offsetsCopy[f]    = 0;
2227:         newOffsetsCopy[f] = 0;
2228:       }
2229:       numColIndices -= 2 * numFields;
2230:       PetscSectionGetOffset(rootIndicesSec, p, &pIndOff);
2231:       pInd = &(rootIndices[pIndOff]);
2232:       PetscSectionGetDof(rootMatricesSec, p, &matSize);
2233:       if (matSize) {
2234:         PetscSectionGetOffset(rootMatricesSec, p, &pMatOff);
2235:         pMat = &rootMatrices[pMatOff];
2236:       }
2237:       PetscSectionGetDof(globalCoarse, p, &dof);
2238:       if (dof < 0) dof = -(dof + 1);
2239:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2240:         PetscInt i, j;
2241:         PetscInt numRowIndices = matSize / numColIndices;

2243:         if (!numRowIndices) { /* don't need to calculate the mat, just the indices */
2244:           PetscInt numIndices, *indices;
2245:           DMPlexGetClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL);
2247:           for (i = 0; i < numColIndices; i++) pInd[i] = indices[i];
2248:           for (i = 0; i < numFields; i++) {
2249:             pInd[numColIndices + i]             = offsets[i + 1];
2250:             pInd[numColIndices + numFields + i] = offsets[i + 1];
2251:           }
2252:           DMPlexRestoreClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL);
2253:         } else {
2254:           PetscInt     closureSize, *closure = NULL, cl;
2255:           PetscScalar *pMatIn, *pMatModified;
2256:           PetscInt     numPoints, *points;

2258:           DMGetWorkArray(coarse, numRowIndices * numRowIndices, MPIU_SCALAR, &pMatIn);
2259:           for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2260:             for (j = 0; j < numRowIndices; j++) pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2261:           }
2262:           DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2263:           for (f = 0; f < maxFields; f++) {
2264:             if (numFields) PetscSectionGetFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]);
2265:             else PetscSectionGetPointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]);
2266:           }
2267:           if (numFields) {
2268:             for (cl = 0; cl < closureSize; cl++) {
2269:               PetscInt c = closure[2 * cl];

2271:               for (f = 0; f < numFields; f++) {
2272:                 PetscInt fDof;

2274:                 PetscSectionGetFieldDof(localCoarse, c, f, &fDof);
2275:                 offsets[f + 1] += fDof;
2276:               }
2277:             }
2278:             for (f = 0; f < numFields; f++) {
2279:               offsets[f + 1] += offsets[f];
2280:               newOffsets[f + 1] = offsets[f + 1];
2281:             }
2282:           }
2283:           /* TODO : flips here ? */
2284:           /* apply hanging node constraints on the right, get the new points and the new offsets */
2285:           DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, perms, pMatIn, &numPoints, NULL, &points, &pMatModified, newOffsets, PETSC_FALSE);
2286:           for (f = 0; f < maxFields; f++) {
2287:             if (numFields) PetscSectionRestoreFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]);
2288:             else PetscSectionRestorePointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]);
2289:           }
2290:           for (f = 0; f < maxFields; f++) {
2291:             if (numFields) PetscSectionGetFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]);
2292:             else PetscSectionGetPointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]);
2293:           }
2294:           if (!numFields) {
2295:             for (i = 0; i < numRowIndices * numColIndices; i++) pMat[i] = pMatModified[i];
2296:           } else {
2297:             PetscInt i, j, count;
2298:             for (f = 0, count = 0; f < numFields; f++) {
2299:               for (i = offsets[f]; i < offsets[f + 1]; i++) {
2300:                 for (j = newOffsets[f]; j < newOffsets[f + 1]; j++, count++) pMat[count] = pMatModified[i * numColIndices + j];
2301:               }
2302:             }
2303:           }
2304:           DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatModified);
2305:           DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2306:           DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatIn);
2307:           if (numFields) {
2308:             for (f = 0; f < numFields; f++) {
2309:               pInd[numColIndices + f]             = offsets[f + 1];
2310:               pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2311:             }
2312:             for (cl = 0; cl < numPoints; cl++) {
2313:               PetscInt globalOff, c = points[2 * cl];
2314:               PetscSectionGetOffset(globalCoarse, c, &globalOff);
2315:               DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd);
2316:             }
2317:           } else {
2318:             for (cl = 0; cl < numPoints; cl++) {
2319:               PetscInt        c    = points[2 * cl], globalOff;
2320:               const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;

2322:               PetscSectionGetOffset(globalCoarse, c, &globalOff);
2323:               DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd);
2324:             }
2325:           }
2326:           for (f = 0; f < maxFields; f++) {
2327:             if (numFields) PetscSectionRestoreFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]);
2328:             else PetscSectionRestorePointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]);
2329:           }
2330:           DMRestoreWorkArray(coarse, numPoints, MPIU_SCALAR, &points);
2331:         }
2332:       } else if (matSize) {
2333:         PetscInt  cOff;
2334:         PetscInt *rowIndices, *colIndices, a, aDof, aOff;

2336:         numRowIndices = matSize / numColIndices;
2338:         DMGetWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices);
2339:         DMGetWorkArray(coarse, numColIndices, MPIU_INT, &colIndices);
2340:         PetscSectionGetOffset(cSec, p, &cOff);
2341:         PetscSectionGetDof(aSec, p, &aDof);
2342:         PetscSectionGetOffset(aSec, p, &aOff);
2343:         if (numFields) {
2344:           for (f = 0; f < numFields; f++) {
2345:             PetscInt fDof;

2347:             PetscSectionGetFieldDof(cSec, p, f, &fDof);
2348:             offsets[f + 1] = fDof;
2349:             for (a = 0; a < aDof; a++) {
2350:               PetscInt anchor = anchors[a + aOff];
2351:               PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof);
2352:               newOffsets[f + 1] += fDof;
2353:             }
2354:           }
2355:           for (f = 0; f < numFields; f++) {
2356:             offsets[f + 1] += offsets[f];
2357:             offsetsCopy[f + 1] = offsets[f + 1];
2358:             newOffsets[f + 1] += newOffsets[f];
2359:             newOffsetsCopy[f + 1] = newOffsets[f + 1];
2360:           }
2361:           DMPlexGetIndicesPointFields_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, -1, NULL, rowIndices);
2362:           for (a = 0; a < aDof; a++) {
2363:             PetscInt anchor = anchors[a + aOff], lOff;
2364:             PetscSectionGetOffset(localCoarse, anchor, &lOff);
2365:             DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, -1, NULL, colIndices);
2366:           }
2367:         } else {
2368:           DMPlexGetIndicesPoint_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, NULL, rowIndices);
2369:           for (a = 0; a < aDof; a++) {
2370:             PetscInt anchor = anchors[a + aOff], lOff;
2371:             PetscSectionGetOffset(localCoarse, anchor, &lOff);
2372:             DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, NULL, colIndices);
2373:           }
2374:         }
2375:         if (numFields) {
2376:           PetscInt count, a;

2378:           for (f = 0, count = 0; f < numFields; f++) {
2379:             PetscInt iSize = offsets[f + 1] - offsets[f];
2380:             PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2381:             MatGetValues(cMat, iSize, &rowIndices[offsets[f]], jSize, &colIndices[newOffsets[f]], &pMat[count]);
2382:             count += iSize * jSize;
2383:             pInd[numColIndices + f]             = offsets[f + 1];
2384:             pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2385:           }
2386:           for (a = 0; a < aDof; a++) {
2387:             PetscInt anchor = anchors[a + aOff];
2388:             PetscInt gOff;
2389:             PetscSectionGetOffset(globalCoarse, anchor, &gOff);
2390:             DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, -1, NULL, pInd);
2391:           }
2392:         } else {
2393:           PetscInt a;
2394:           MatGetValues(cMat, numRowIndices, rowIndices, numColIndices, colIndices, pMat);
2395:           for (a = 0; a < aDof; a++) {
2396:             PetscInt anchor = anchors[a + aOff];
2397:             PetscInt gOff;
2398:             PetscSectionGetOffset(globalCoarse, anchor, &gOff);
2399:             DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, NULL, pInd);
2400:           }
2401:         }
2402:         DMRestoreWorkArray(coarse, numColIndices, MPIU_INT, &colIndices);
2403:         DMRestoreWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices);
2404:       } else {
2405:         PetscInt gOff;

2407:         PetscSectionGetOffset(globalCoarse, p, &gOff);
2408:         if (numFields) {
2409:           for (f = 0; f < numFields; f++) {
2410:             PetscInt fDof;
2411:             PetscSectionGetFieldDof(localCoarse, p, f, &fDof);
2412:             offsets[f + 1] = fDof + offsets[f];
2413:           }
2414:           for (f = 0; f < numFields; f++) {
2415:             pInd[numColIndices + f]             = offsets[f + 1];
2416:             pInd[numColIndices + numFields + f] = offsets[f + 1];
2417:           }
2418:           DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd);
2419:         } else {
2420:           DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd);
2421:         }
2422:       }
2423:     }
2424:     PetscFree(maxChildIds);
2425:   }
2426:   {
2427:     PetscSF   indicesSF, matricesSF;
2428:     PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;

2430:     PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec);
2431:     PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafMatricesSec);
2432:     PetscSFDistributeSection(coarseToFineEmbedded, rootIndicesSec, &remoteOffsetsIndices, leafIndicesSec);
2433:     PetscSFDistributeSection(coarseToFineEmbedded, rootMatricesSec, &remoteOffsetsMatrices, leafMatricesSec);
2434:     PetscSFCreateSectionSF(coarseToFineEmbedded, rootIndicesSec, remoteOffsetsIndices, leafIndicesSec, &indicesSF);
2435:     PetscSFCreateSectionSF(coarseToFineEmbedded, rootMatricesSec, remoteOffsetsMatrices, leafMatricesSec, &matricesSF);
2436:     PetscSFDestroy(&coarseToFineEmbedded);
2437:     PetscFree(remoteOffsetsIndices);
2438:     PetscFree(remoteOffsetsMatrices);
2439:     PetscSectionGetStorageSize(leafIndicesSec, &numLeafIndices);
2440:     PetscSectionGetStorageSize(leafMatricesSec, &numLeafMatrices);
2441:     PetscMalloc2(numLeafIndices, &leafIndices, numLeafMatrices, &leafMatrices);
2442:     PetscSFBcastBegin(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE);
2443:     PetscSFBcastBegin(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE);
2444:     PetscSFBcastEnd(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE);
2445:     PetscSFBcastEnd(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE);
2446:     PetscSFDestroy(&matricesSF);
2447:     PetscSFDestroy(&indicesSF);
2448:     PetscFree2(rootIndices, rootMatrices);
2449:     PetscSectionDestroy(&rootIndicesSec);
2450:     PetscSectionDestroy(&rootMatricesSec);
2451:   }
2452:   /* count to preallocate */
2453:   DMGetLocalSection(fine, &localFine);
2454:   {
2455:     PetscInt       nGlobal;
2456:     PetscInt      *dnnz, *onnz;
2457:     PetscLayout    rowMap, colMap;
2458:     PetscInt       rowStart, rowEnd, colStart, colEnd;
2459:     PetscInt       maxDof;
2460:     PetscInt      *rowIndices;
2461:     DM             refTree;
2462:     PetscInt     **refPointFieldN;
2463:     PetscScalar ***refPointFieldMats;
2464:     PetscSection   refConSec, refAnSec;
2465:     PetscInt       pRefStart, pRefEnd, maxConDof, maxColumns, leafStart, leafEnd;
2466:     PetscScalar   *pointWork;

2468:     PetscSectionGetConstrainedStorageSize(globalFine, &nGlobal);
2469:     PetscCalloc2(nGlobal, &dnnz, nGlobal, &onnz);
2470:     MatGetLayouts(mat, &rowMap, &colMap);
2471:     PetscLayoutSetUp(rowMap);
2472:     PetscLayoutSetUp(colMap);
2473:     PetscLayoutGetRange(rowMap, &rowStart, &rowEnd);
2474:     PetscLayoutGetRange(colMap, &colStart, &colEnd);
2475:     PetscSectionGetMaxDof(localFine, &maxDof);
2476:     PetscSectionGetChart(leafIndicesSec, &leafStart, &leafEnd);
2477:     DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices);
2478:     for (p = leafStart; p < leafEnd; p++) {
2479:       PetscInt gDof, gcDof, gOff;
2480:       PetscInt numColIndices, pIndOff, *pInd;
2481:       PetscInt matSize;
2482:       PetscInt i;

2484:       PetscSectionGetDof(globalFine, p, &gDof);
2485:       PetscSectionGetConstraintDof(globalFine, p, &gcDof);
2486:       if ((gDof - gcDof) <= 0) continue;
2487:       PetscSectionGetOffset(globalFine, p, &gOff);
2490:       PetscSectionGetDof(leafIndicesSec, p, &numColIndices);
2491:       PetscSectionGetOffset(leafIndicesSec, p, &pIndOff);
2492:       numColIndices -= 2 * numFields;
2494:       pInd              = &leafIndices[pIndOff];
2495:       offsets[0]        = 0;
2496:       offsetsCopy[0]    = 0;
2497:       newOffsets[0]     = 0;
2498:       newOffsetsCopy[0] = 0;
2499:       if (numFields) {
2500:         PetscInt f;
2501:         for (f = 0; f < numFields; f++) {
2502:           PetscInt rowDof;

2504:           PetscSectionGetFieldDof(localFine, p, f, &rowDof);
2505:           offsets[f + 1]     = offsets[f] + rowDof;
2506:           offsetsCopy[f + 1] = offsets[f + 1];
2507:           newOffsets[f + 1]  = pInd[numColIndices + numFields + f];
2508:           numD[f]            = 0;
2509:           numO[f]            = 0;
2510:         }
2511:         DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices);
2512:         for (f = 0; f < numFields; f++) {
2513:           PetscInt colOffset    = newOffsets[f];
2514:           PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];

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

2519:             if (gInd >= colStart && gInd < colEnd) {
2520:               numD[f]++;
2521:             } else if (gInd >= 0) { /* negative means non-entry */
2522:               numO[f]++;
2523:             }
2524:           }
2525:         }
2526:       } else {
2527:         DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices);
2528:         numD[0] = 0;
2529:         numO[0] = 0;
2530:         for (i = 0; i < numColIndices; i++) {
2531:           PetscInt gInd = pInd[i];

2533:           if (gInd >= colStart && gInd < colEnd) {
2534:             numD[0]++;
2535:           } else if (gInd >= 0) { /* negative means non-entry */
2536:             numO[0]++;
2537:           }
2538:         }
2539:       }
2540:       PetscSectionGetDof(leafMatricesSec, p, &matSize);
2541:       if (!matSize) { /* incoming matrix is identity */
2542:         PetscInt childId;

2544:         childId = childIds[p - pStartF];
2545:         if (childId < 0) { /* no child interpolation: one nnz per */
2546:           if (numFields) {
2547:             PetscInt f;
2548:             for (f = 0; f < numFields; f++) {
2549:               PetscInt numRows = offsets[f + 1] - offsets[f], row;
2550:               for (row = 0; row < numRows; row++) {
2551:                 PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2552:                 PetscInt gIndFine   = rowIndices[offsets[f] + row];
2553:                 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2555:                   dnnz[gIndFine - rowStart] = 1;
2556:                 } else if (gIndCoarse >= 0) { /* remote */
2558:                   onnz[gIndFine - rowStart] = 1;
2559:                 } else { /* constrained */
2561:                 }
2562:               }
2563:             }
2564:           } else {
2565:             PetscInt i;
2566:             for (i = 0; i < gDof; i++) {
2567:               PetscInt gIndCoarse = pInd[i];
2568:               PetscInt gIndFine   = rowIndices[i];
2569:               if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2571:                 dnnz[gIndFine - rowStart] = 1;
2572:               } else if (gIndCoarse >= 0) { /* remote */
2574:                 onnz[gIndFine - rowStart] = 1;
2575:               } else { /* constrained */
2577:               }
2578:             }
2579:           }
2580:         } else { /* interpolate from all */
2581:           if (numFields) {
2582:             PetscInt f;
2583:             for (f = 0; f < numFields; f++) {
2584:               PetscInt numRows = offsets[f + 1] - offsets[f], row;
2585:               for (row = 0; row < numRows; row++) {
2586:                 PetscInt gIndFine = rowIndices[offsets[f] + row];
2587:                 if (gIndFine >= 0) {
2589:                   dnnz[gIndFine - rowStart] = numD[f];
2590:                   onnz[gIndFine - rowStart] = numO[f];
2591:                 }
2592:               }
2593:             }
2594:           } else {
2595:             PetscInt i;
2596:             for (i = 0; i < gDof; i++) {
2597:               PetscInt gIndFine = rowIndices[i];
2598:               if (gIndFine >= 0) {
2600:                 dnnz[gIndFine - rowStart] = numD[0];
2601:                 onnz[gIndFine - rowStart] = numO[0];
2602:               }
2603:             }
2604:           }
2605:         }
2606:       } else { /* interpolate from all */
2607:         if (numFields) {
2608:           PetscInt f;
2609:           for (f = 0; f < numFields; f++) {
2610:             PetscInt numRows = offsets[f + 1] - offsets[f], row;
2611:             for (row = 0; row < numRows; row++) {
2612:               PetscInt gIndFine = rowIndices[offsets[f] + row];
2613:               if (gIndFine >= 0) {
2615:                 dnnz[gIndFine - rowStart] = numD[f];
2616:                 onnz[gIndFine - rowStart] = numO[f];
2617:               }
2618:             }
2619:           }
2620:         } else { /* every dof get a full row */
2621:           PetscInt i;
2622:           for (i = 0; i < gDof; i++) {
2623:             PetscInt gIndFine = rowIndices[i];
2624:             if (gIndFine >= 0) {
2626:               dnnz[gIndFine - rowStart] = numD[0];
2627:               onnz[gIndFine - rowStart] = numO[0];
2628:             }
2629:           }
2630:         }
2631:       }
2632:     }
2633:     MatXAIJSetPreallocation(mat, 1, dnnz, onnz, NULL, NULL);
2634:     PetscFree2(dnnz, onnz);

2636:     DMPlexGetReferenceTree(fine, &refTree);
2637:     DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN);
2638:     DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL);
2639:     DMPlexGetAnchors(refTree, &refAnSec, NULL);
2640:     PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd);
2641:     PetscSectionGetMaxDof(refConSec, &maxConDof);
2642:     PetscSectionGetMaxDof(leafIndicesSec, &maxColumns);
2643:     PetscMalloc1(maxConDof * maxColumns, &pointWork);
2644:     for (p = leafStart; p < leafEnd; p++) {
2645:       PetscInt gDof, gcDof, gOff;
2646:       PetscInt numColIndices, pIndOff, *pInd;
2647:       PetscInt matSize;
2648:       PetscInt childId;

2650:       PetscSectionGetDof(globalFine, p, &gDof);
2651:       PetscSectionGetConstraintDof(globalFine, p, &gcDof);
2652:       if ((gDof - gcDof) <= 0) continue;
2653:       childId = childIds[p - pStartF];
2654:       PetscSectionGetOffset(globalFine, p, &gOff);
2655:       PetscSectionGetDof(leafIndicesSec, p, &numColIndices);
2656:       PetscSectionGetOffset(leafIndicesSec, p, &pIndOff);
2657:       numColIndices -= 2 * numFields;
2658:       pInd              = &leafIndices[pIndOff];
2659:       offsets[0]        = 0;
2660:       offsetsCopy[0]    = 0;
2661:       newOffsets[0]     = 0;
2662:       newOffsetsCopy[0] = 0;
2663:       rowOffsets[0]     = 0;
2664:       if (numFields) {
2665:         PetscInt f;
2666:         for (f = 0; f < numFields; f++) {
2667:           PetscInt rowDof;

2669:           PetscSectionGetFieldDof(localFine, p, f, &rowDof);
2670:           offsets[f + 1]     = offsets[f] + rowDof;
2671:           offsetsCopy[f + 1] = offsets[f + 1];
2672:           rowOffsets[f + 1]  = pInd[numColIndices + f];
2673:           newOffsets[f + 1]  = pInd[numColIndices + numFields + f];
2674:         }
2675:         DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices);
2676:       } else {
2677:         DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices);
2678:       }
2679:       PetscSectionGetDof(leafMatricesSec, p, &matSize);
2680:       if (!matSize) {      /* incoming matrix is identity */
2681:         if (childId < 0) { /* no child interpolation: scatter */
2682:           if (numFields) {
2683:             PetscInt f;
2684:             for (f = 0; f < numFields; f++) {
2685:               PetscInt numRows = offsets[f + 1] - offsets[f], row;
2686:               for (row = 0; row < numRows; row++) MatSetValue(mat, rowIndices[offsets[f] + row], pInd[newOffsets[f] + row], 1., INSERT_VALUES);
2687:             }
2688:           } else {
2689:             PetscInt numRows = gDof, row;
2690:             for (row = 0; row < numRows; row++) MatSetValue(mat, rowIndices[row], pInd[row], 1., INSERT_VALUES);
2691:           }
2692:         } else { /* interpolate from all */
2693:           if (numFields) {
2694:             PetscInt f;
2695:             for (f = 0; f < numFields; f++) {
2696:               PetscInt numRows = offsets[f + 1] - offsets[f];
2697:               PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2698:               MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], refPointFieldMats[childId - pRefStart][f], INSERT_VALUES);
2699:             }
2700:           } else {
2701:             MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, refPointFieldMats[childId - pRefStart][0], INSERT_VALUES);
2702:           }
2703:         }
2704:       } else { /* interpolate from all */
2705:         PetscInt     pMatOff;
2706:         PetscScalar *pMat;

2708:         PetscSectionGetOffset(leafMatricesSec, p, &pMatOff);
2709:         pMat = &leafMatrices[pMatOff];
2710:         if (childId < 0) { /* copy the incoming matrix */
2711:           if (numFields) {
2712:             PetscInt f, count;
2713:             for (f = 0, count = 0; f < numFields; f++) {
2714:               PetscInt     numRows   = offsets[f + 1] - offsets[f];
2715:               PetscInt     numCols   = newOffsets[f + 1] - newOffsets[f];
2716:               PetscInt     numInRows = rowOffsets[f + 1] - rowOffsets[f];
2717:               PetscScalar *inMat     = &pMat[count];

2719:               MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], inMat, INSERT_VALUES);
2720:               count += numCols * numInRows;
2721:             }
2722:           } else {
2723:             MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, pMat, INSERT_VALUES);
2724:           }
2725:         } else { /* multiply the incoming matrix by the child interpolation */
2726:           if (numFields) {
2727:             PetscInt f, count;
2728:             for (f = 0, count = 0; f < numFields; f++) {
2729:               PetscInt     numRows   = offsets[f + 1] - offsets[f];
2730:               PetscInt     numCols   = newOffsets[f + 1] - newOffsets[f];
2731:               PetscInt     numInRows = rowOffsets[f + 1] - rowOffsets[f];
2732:               PetscScalar *inMat     = &pMat[count];
2733:               PetscInt     i, j, k;
2735:               for (i = 0; i < numRows; i++) {
2736:                 for (j = 0; j < numCols; j++) {
2737:                   PetscScalar val = 0.;
2738:                   for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2739:                   pointWork[i * numCols + j] = val;
2740:                 }
2741:               }
2742:               MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], pointWork, INSERT_VALUES);
2743:               count += numCols * numInRows;
2744:             }
2745:           } else { /* every dof gets a full row */
2746:             PetscInt numRows   = gDof;
2747:             PetscInt numCols   = numColIndices;
2748:             PetscInt numInRows = matSize / numColIndices;
2749:             PetscInt i, j, k;
2751:             for (i = 0; i < numRows; i++) {
2752:               for (j = 0; j < numCols; j++) {
2753:                 PetscScalar val = 0.;
2754:                 for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2755:                 pointWork[i * numCols + j] = val;
2756:               }
2757:             }
2758:             MatSetValues(mat, numRows, rowIndices, numCols, pInd, pointWork, INSERT_VALUES);
2759:           }
2760:         }
2761:       }
2762:     }
2763:     DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN);
2764:     DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices);
2765:     PetscFree(pointWork);
2766:   }
2767:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2768:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
2769:   PetscSectionDestroy(&leafIndicesSec);
2770:   PetscSectionDestroy(&leafMatricesSec);
2771:   PetscFree2(leafIndices, leafMatrices);
2772:   PetscFree2(*(PetscInt ****)&perms, *(PetscScalar ****)&flips);
2773:   PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO);
2774:   ISRestoreIndices(aIS, &anchors);
2775:   return 0;
2776: }

2778: /*
2779:  * Assuming a nodal basis (w.r.t. the dual basis) basis:
2780:  *
2781:  * for each coarse dof \phi^c_i:
2782:  *   for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
2783:  *     for each fine dof \phi^f_j;
2784:  *       a_{i,j} = 0;
2785:  *       for each fine dof \phi^f_k:
2786:  *         a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
2787:  *                    [^^^ this is = \phi^c_i ^^^]
2788:  */
2789: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
2790: {
2791:   PetscDS      ds;
2792:   PetscSection section, cSection;
2793:   DMLabel      canonical, depth;
2794:   Mat          cMat, mat;
2795:   PetscInt    *nnz;
2796:   PetscInt     f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
2797:   PetscInt     m, n;
2798:   PetscScalar *pointScalar;
2799:   PetscReal   *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;

2801:   DMGetLocalSection(refTree, &section);
2802:   DMGetDimension(refTree, &dim);
2803:   PetscMalloc6(dim, &v0, dim, &v0parent, dim, &vtmp, dim * dim, &J, dim * dim, &Jparent, dim * dim, &invJ);
2804:   PetscMalloc2(dim, &pointScalar, dim, &pointRef);
2805:   DMGetDS(refTree, &ds);
2806:   PetscDSGetNumFields(ds, &numFields);
2807:   PetscSectionGetNumFields(section, &numSecFields);
2808:   DMGetLabel(refTree, "canonical", &canonical);
2809:   DMGetLabel(refTree, "depth", &depth);
2810:   DMGetDefaultConstraints(refTree, &cSection, &cMat, NULL);
2811:   DMPlexGetChart(refTree, &pStart, &pEnd);
2812:   DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);
2813:   MatGetSize(cMat, &n, &m); /* the injector has transpose sizes from the constraint matrix */
2814:   /* Step 1: compute non-zero pattern.  A proper subset of constraint matrix non-zero */
2815:   PetscCalloc1(m, &nnz);
2816:   for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
2817:     const PetscInt *children;
2818:     PetscInt        numChildren;
2819:     PetscInt        i, numChildDof, numSelfDof;

2821:     if (canonical) {
2822:       PetscInt pCanonical;
2823:       DMLabelGetValue(canonical, p, &pCanonical);
2824:       if (p != pCanonical) continue;
2825:     }
2826:     DMPlexGetTreeChildren(refTree, p, &numChildren, &children);
2827:     if (!numChildren) continue;
2828:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
2829:       PetscInt child = children[i];
2830:       PetscInt dof;

2832:       PetscSectionGetDof(section, child, &dof);
2833:       numChildDof += dof;
2834:     }
2835:     PetscSectionGetDof(section, p, &numSelfDof);
2836:     if (!numChildDof || !numSelfDof) continue;
2837:     for (f = 0; f < numFields; f++) {
2838:       PetscInt selfOff;

2840:       if (numSecFields) { /* count the dofs for just this field */
2841:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
2842:           PetscInt child = children[i];
2843:           PetscInt dof;

2845:           PetscSectionGetFieldDof(section, child, f, &dof);
2846:           numChildDof += dof;
2847:         }
2848:         PetscSectionGetFieldDof(section, p, f, &numSelfDof);
2849:         PetscSectionGetFieldOffset(section, p, f, &selfOff);
2850:       } else {
2851:         PetscSectionGetOffset(section, p, &selfOff);
2852:       }
2853:       for (i = 0; i < numSelfDof; i++) nnz[selfOff + i] = numChildDof;
2854:     }
2855:   }
2856:   MatCreateAIJ(PETSC_COMM_SELF, m, n, m, n, -1, nnz, -1, NULL, &mat);
2857:   PetscFree(nnz);
2858:   /* Setp 2: compute entries */
2859:   for (p = pStart; p < pEnd; p++) {
2860:     const PetscInt *children;
2861:     PetscInt        numChildren;
2862:     PetscInt        i, numChildDof, numSelfDof;

2864:     /* same conditions about when entries occur */
2865:     if (canonical) {
2866:       PetscInt pCanonical;
2867:       DMLabelGetValue(canonical, p, &pCanonical);
2868:       if (p != pCanonical) continue;
2869:     }
2870:     DMPlexGetTreeChildren(refTree, p, &numChildren, &children);
2871:     if (!numChildren) continue;
2872:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
2873:       PetscInt child = children[i];
2874:       PetscInt dof;

2876:       PetscSectionGetDof(section, child, &dof);
2877:       numChildDof += dof;
2878:     }
2879:     PetscSectionGetDof(section, p, &numSelfDof);
2880:     if (!numChildDof || !numSelfDof) continue;

2882:     for (f = 0; f < numFields; f++) {
2883:       PetscInt        pI = -1, cI = -1;
2884:       PetscInt        selfOff, Nc, parentCell;
2885:       PetscInt        cellShapeOff;
2886:       PetscObject     disc;
2887:       PetscDualSpace  dsp;
2888:       PetscClassId    classId;
2889:       PetscScalar    *pointMat;
2890:       PetscInt       *matRows, *matCols;
2891:       PetscInt        pO = PETSC_MIN_INT;
2892:       const PetscInt *depthNumDof;

2894:       if (numSecFields) {
2895:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
2896:           PetscInt child = children[i];
2897:           PetscInt dof;

2899:           PetscSectionGetFieldDof(section, child, f, &dof);
2900:           numChildDof += dof;
2901:         }
2902:         PetscSectionGetFieldDof(section, p, f, &numSelfDof);
2903:         PetscSectionGetFieldOffset(section, p, f, &selfOff);
2904:       } else {
2905:         PetscSectionGetOffset(section, p, &selfOff);
2906:       }

2908:       /* find a cell whose closure contains p */
2909:       if (p >= cStart && p < cEnd) {
2910:         parentCell = p;
2911:       } else {
2912:         PetscInt *star = NULL;
2913:         PetscInt  numStar;

2915:         parentCell = -1;
2916:         DMPlexGetTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star);
2917:         for (i = numStar - 1; i >= 0; i--) {
2918:           PetscInt c = star[2 * i];

2920:           if (c >= cStart && c < cEnd) {
2921:             parentCell = c;
2922:             break;
2923:           }
2924:         }
2925:         DMPlexRestoreTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star);
2926:       }
2927:       /* determine the offset of p's shape functions within parentCell's shape functions */
2928:       PetscDSGetDiscretization(ds, f, &disc);
2929:       PetscObjectGetClassId(disc, &classId);
2930:       if (classId == PETSCFE_CLASSID) {
2931:         PetscFEGetDualSpace((PetscFE)disc, &dsp);
2932:       } else if (classId == PETSCFV_CLASSID) {
2933:         PetscFVGetDualSpace((PetscFV)disc, &dsp);
2934:       } else {
2935:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported discretization object");
2936:       }
2937:       PetscDualSpaceGetNumDof(dsp, &depthNumDof);
2938:       PetscDualSpaceGetNumComponents(dsp, &Nc);
2939:       {
2940:         PetscInt *closure = NULL;
2941:         PetscInt  numClosure;

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

2947:           pO = closure[2 * i + 1];
2948:           if (point == p) {
2949:             pI = i;
2950:             break;
2951:           }
2952:           DMLabelGetValue(depth, point, &pointDepth);
2953:           cellShapeOff += depthNumDof[pointDepth];
2954:         }
2955:         DMPlexRestoreTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure);
2956:       }

2958:       DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat);
2959:       DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows);
2960:       matCols = matRows + numSelfDof;
2961:       for (i = 0; i < numSelfDof; i++) matRows[i] = selfOff + i;
2962:       for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
2963:       {
2964:         PetscInt colOff = 0;

2966:         for (i = 0; i < numChildren; i++) {
2967:           PetscInt child = children[i];
2968:           PetscInt dof, off, j;

2970:           if (numSecFields) {
2971:             PetscSectionGetFieldDof(cSection, child, f, &dof);
2972:             PetscSectionGetFieldOffset(cSection, child, f, &off);
2973:           } else {
2974:             PetscSectionGetDof(cSection, child, &dof);
2975:             PetscSectionGetOffset(cSection, child, &off);
2976:           }

2978:           for (j = 0; j < dof; j++) matCols[colOff++] = off + j;
2979:         }
2980:       }
2981:       if (classId == PETSCFE_CLASSID) {
2982:         PetscFE              fe = (PetscFE)disc;
2983:         PetscInt             fSize;
2984:         const PetscInt    ***perms;
2985:         const PetscScalar ***flips;
2986:         const PetscInt      *pperms;

2988:         PetscFEGetDualSpace(fe, &dsp);
2989:         PetscDualSpaceGetDimension(dsp, &fSize);
2990:         PetscDualSpaceGetSymmetries(dsp, &perms, &flips);
2991:         pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL;
2992:         for (i = 0; i < numSelfDof; i++) { /* for every shape function */
2993:           PetscQuadrature  q;
2994:           PetscInt         dim, thisNc, numPoints, j, k;
2995:           const PetscReal *points;
2996:           const PetscReal *weights;
2997:           PetscInt        *closure = NULL;
2998:           PetscInt         numClosure;
2999:           PetscInt         iCell              = pperms ? pperms[i] : i;
3000:           PetscInt         parentCellShapeDof = cellShapeOff + iCell;
3001:           PetscTabulation  Tparent;

3003:           PetscDualSpaceGetFunctional(dsp, parentCellShapeDof, &q);
3004:           PetscQuadratureGetData(q, &dim, &thisNc, &numPoints, &points, &weights);
3006:           PetscFECreateTabulation(fe, 1, numPoints, points, 0, &Tparent); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3007:           for (j = 0; j < numPoints; j++) {
3008:             PetscInt           childCell = -1;
3009:             PetscReal         *parentValAtPoint;
3010:             const PetscReal    xi0[3]    = {-1., -1., -1.};
3011:             const PetscReal   *pointReal = &points[dim * j];
3012:             const PetscScalar *point;
3013:             PetscTabulation    Tchild;
3014:             PetscInt           childCellShapeOff, pointMatOff;
3015: #if defined(PETSC_USE_COMPLEX)
3016:             PetscInt d;

3018:             for (d = 0; d < dim; d++) pointScalar[d] = points[dim * j + d];
3019:             point = pointScalar;
3020: #else
3021:             point = pointReal;
3022: #endif

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

3026:             for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3027:               PetscInt  child = children[k];
3028:               PetscInt *star  = NULL;
3029:               PetscInt  numStar, s;

3031:               DMPlexGetTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star);
3032:               for (s = numStar - 1; s >= 0; s--) {
3033:                 PetscInt c = star[2 * s];

3035:                 if (c < cStart || c >= cEnd) continue;
3036:                 DMPlexLocatePoint_Internal(refTree, dim, point, c, &childCell);
3037:                 if (childCell >= 0) break;
3038:               }
3039:               DMPlexRestoreTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star);
3040:               if (childCell >= 0) break;
3041:             }
3043:             DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);
3044:             DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);
3045:             CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
3046:             CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);

3048:             PetscFECreateTabulation(fe, 1, 1, pointRef, 0, &Tchild);
3049:             DMPlexGetTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure);
3050:             for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3051:               PetscInt        child = children[k], childDepth, childDof, childO = PETSC_MIN_INT;
3052:               PetscInt        l;
3053:               const PetscInt *cperms;

3055:               DMLabelGetValue(depth, child, &childDepth);
3056:               childDof = depthNumDof[childDepth];
3057:               for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) {
3058:                 PetscInt point = closure[2 * l];
3059:                 PetscInt pointDepth;

3061:                 childO = closure[2 * l + 1];
3062:                 if (point == child) {
3063:                   cI = l;
3064:                   break;
3065:                 }
3066:                 DMLabelGetValue(depth, point, &pointDepth);
3067:                 childCellShapeOff += depthNumDof[pointDepth];
3068:               }
3069:               if (l == numClosure) {
3070:                 pointMatOff += childDof;
3071:                 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3072:               }
3073:               cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL;
3074:               for (l = 0; l < childDof; l++) {
3075:                 PetscInt   lCell        = cperms ? cperms[l] : l;
3076:                 PetscInt   childCellDof = childCellShapeOff + lCell;
3077:                 PetscReal *childValAtPoint;
3078:                 PetscReal  val = 0.;

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

3083:                 pointMat[i * numChildDof + pointMatOff + l] += val;
3084:               }
3085:               pointMatOff += childDof;
3086:             }
3087:             DMPlexRestoreTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure);
3088:             PetscTabulationDestroy(&Tchild);
3089:           }
3090:           PetscTabulationDestroy(&Tparent);
3091:         }
3092:       } else { /* just the volume-weighted averages of the children */
3093:         PetscReal parentVol;
3094:         PetscInt  childCell;

3096:         DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);
3097:         for (i = 0, childCell = 0; i < numChildren; i++) {
3098:           PetscInt  child = children[i], j;
3099:           PetscReal childVol;

3101:           if (child < cStart || child >= cEnd) continue;
3102:           DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);
3103:           for (j = 0; j < Nc; j++) pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3104:           childCell++;
3105:         }
3106:       }
3107:       /* Insert pointMat into mat */
3108:       MatSetValues(mat, numSelfDof, matRows, numChildDof, matCols, pointMat, INSERT_VALUES);
3109:       DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows);
3110:       DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat);
3111:     }
3112:   }
3113:   PetscFree6(v0, v0parent, vtmp, J, Jparent, invJ);
3114:   PetscFree2(pointScalar, pointRef);
3115:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
3116:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
3117:   *inj = mat;
3118:   return 0;
3119: }

3121: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3122: {
3123:   PetscDS        ds;
3124:   PetscInt       numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3125:   PetscScalar ***refPointFieldMats;
3126:   PetscSection   refConSec, refSection;

3128:   DMGetDS(refTree, &ds);
3129:   PetscDSGetNumFields(ds, &numFields);
3130:   DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL);
3131:   DMGetLocalSection(refTree, &refSection);
3132:   PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd);
3133:   PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats);
3134:   PetscSectionGetMaxDof(refConSec, &maxDof);
3135:   PetscMalloc1(maxDof, &rows);
3136:   PetscMalloc1(maxDof * maxDof, &cols);
3137:   for (p = pRefStart; p < pRefEnd; p++) {
3138:     PetscInt parent, pDof, parentDof;

3140:     DMPlexGetTreeParent(refTree, p, &parent, NULL);
3141:     PetscSectionGetDof(refConSec, p, &pDof);
3142:     PetscSectionGetDof(refSection, parent, &parentDof);
3143:     if (!pDof || !parentDof || parent == p) continue;

3145:     PetscMalloc1(numFields, &refPointFieldMats[p - pRefStart]);
3146:     for (f = 0; f < numFields; f++) {
3147:       PetscInt cDof, cOff, numCols, r;

3149:       if (numFields > 1) {
3150:         PetscSectionGetFieldDof(refConSec, p, f, &cDof);
3151:         PetscSectionGetFieldOffset(refConSec, p, f, &cOff);
3152:       } else {
3153:         PetscSectionGetDof(refConSec, p, &cDof);
3154:         PetscSectionGetOffset(refConSec, p, &cOff);
3155:       }

3157:       for (r = 0; r < cDof; r++) rows[r] = cOff + r;
3158:       numCols = 0;
3159:       {
3160:         PetscInt aDof, aOff, j;

3162:         if (numFields > 1) {
3163:           PetscSectionGetFieldDof(refSection, parent, f, &aDof);
3164:           PetscSectionGetFieldOffset(refSection, parent, f, &aOff);
3165:         } else {
3166:           PetscSectionGetDof(refSection, parent, &aDof);
3167:           PetscSectionGetOffset(refSection, parent, &aOff);
3168:         }

3170:         for (j = 0; j < aDof; j++) cols[numCols++] = aOff + j;
3171:       }
3172:       PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]);
3173:       /* transpose of constraint matrix */
3174:       MatGetValues(inj, numCols, cols, cDof, rows, refPointFieldMats[p - pRefStart][f]);
3175:     }
3176:   }
3177:   *childrenMats = refPointFieldMats;
3178:   PetscFree(rows);
3179:   PetscFree(cols);
3180:   return 0;
3181: }

3183: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3184: {
3185:   PetscDS        ds;
3186:   PetscScalar ***refPointFieldMats;
3187:   PetscInt       numFields, pRefStart, pRefEnd, p, f;
3188:   PetscSection   refConSec, refSection;

3190:   refPointFieldMats = *childrenMats;
3191:   *childrenMats     = NULL;
3192:   DMGetDS(refTree, &ds);
3193:   DMGetLocalSection(refTree, &refSection);
3194:   PetscDSGetNumFields(ds, &numFields);
3195:   DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL);
3196:   PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd);
3197:   for (p = pRefStart; p < pRefEnd; p++) {
3198:     PetscInt parent, pDof, parentDof;

3200:     DMPlexGetTreeParent(refTree, p, &parent, NULL);
3201:     PetscSectionGetDof(refConSec, p, &pDof);
3202:     PetscSectionGetDof(refSection, parent, &parentDof);
3203:     if (!pDof || !parentDof || parent == p) continue;

3205:     for (f = 0; f < numFields; f++) {
3206:       PetscInt cDof;

3208:       if (numFields > 1) {
3209:         PetscSectionGetFieldDof(refConSec, p, f, &cDof);
3210:       } else {
3211:         PetscSectionGetDof(refConSec, p, &cDof);
3212:       }

3214:       PetscFree(refPointFieldMats[p - pRefStart][f]);
3215:     }
3216:     PetscFree(refPointFieldMats[p - pRefStart]);
3217:   }
3218:   PetscFree(refPointFieldMats);
3219:   return 0;
3220: }

3222: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree, Mat *injRef)
3223: {
3224:   Mat         cMatRef;
3225:   PetscObject injRefObj;

3227:   DMGetDefaultConstraints(refTree, NULL, &cMatRef, NULL);
3228:   PetscObjectQuery((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", &injRefObj);
3229:   *injRef = (Mat)injRefObj;
3230:   if (!*injRef) {
3231:     DMPlexComputeInjectorReferenceTree(refTree, injRef);
3232:     PetscObjectCompose((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", (PetscObject)*injRef);
3233:     /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3234:     PetscObjectDereference((PetscObject)*injRef);
3235:   }
3236:   return 0;
3237: }

3239: 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)
3240: {
3241:   PetscInt        pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3242:   PetscSection    globalCoarse, globalFine;
3243:   PetscSection    localCoarse, localFine, leafIndicesSec;
3244:   PetscSection    multiRootSec, rootIndicesSec;
3245:   PetscInt       *leafInds, *rootInds = NULL;
3246:   const PetscInt *rootDegrees;
3247:   PetscScalar    *leafVals = NULL, *rootVals = NULL;
3248:   PetscSF         coarseToFineEmbedded;

3250:   DMPlexGetChart(coarse, &pStartC, &pEndC);
3251:   DMPlexGetChart(fine, &pStartF, &pEndF);
3252:   DMGetLocalSection(fine, &localFine);
3253:   DMGetGlobalSection(fine, &globalFine);
3254:   PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec);
3255:   PetscSectionSetChart(leafIndicesSec, pStartF, pEndF);
3256:   PetscSectionGetMaxDof(localFine, &maxDof);
3257:   { /* winnow fine points that don't have global dofs out of the sf */
3258:     PetscInt        l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3259:     const PetscInt *leaves;

3261:     PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL);
3262:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3263:       p = leaves ? leaves[l] : l;
3264:       PetscSectionGetDof(globalFine, p, &dof);
3265:       PetscSectionGetConstraintDof(globalFine, p, &cdof);
3266:       if ((dof - cdof) > 0) {
3267:         numPointsWithDofs++;

3269:         PetscSectionGetDof(localFine, p, &dof);
3270:         PetscSectionSetDof(leafIndicesSec, p, dof + 1);
3271:       }
3272:     }
3273:     PetscMalloc1(numPointsWithDofs, &pointsWithDofs);
3274:     PetscSectionSetUp(leafIndicesSec);
3275:     PetscSectionGetStorageSize(leafIndicesSec, &numIndices);
3276:     PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1), &leafInds);
3277:     if (gatheredValues) PetscMalloc1(numIndices, &leafVals);
3278:     for (l = 0, offset = 0; l < nleaves; l++) {
3279:       p = leaves ? leaves[l] : l;
3280:       PetscSectionGetDof(globalFine, p, &dof);
3281:       PetscSectionGetConstraintDof(globalFine, p, &cdof);
3282:       if ((dof - cdof) > 0) {
3283:         PetscInt     off, gOff;
3284:         PetscInt    *pInd;
3285:         PetscScalar *pVal = NULL;

3287:         pointsWithDofs[offset++] = l;

3289:         PetscSectionGetOffset(leafIndicesSec, p, &off);

3291:         pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3292:         if (gatheredValues) {
3293:           PetscInt i;

3295:           pVal = &leafVals[off + 1];
3296:           for (i = 0; i < dof; i++) pVal[i] = 0.;
3297:         }
3298:         PetscSectionGetOffset(globalFine, p, &gOff);

3300:         offsets[0] = 0;
3301:         if (numFields) {
3302:           PetscInt f;

3304:           for (f = 0; f < numFields; f++) {
3305:             PetscInt fDof;
3306:             PetscSectionGetFieldDof(localFine, p, f, &fDof);
3307:             offsets[f + 1] = fDof + offsets[f];
3308:           }
3309:           DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd);
3310:         } else {
3311:           DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd);
3312:         }
3313:         if (gatheredValues) VecGetValues(fineVec, dof, pInd, pVal);
3314:       }
3315:     }
3316:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
3317:     PetscFree(pointsWithDofs);
3318:   }

3320:   DMPlexGetChart(coarse, &pStartC, &pEndC);
3321:   DMGetLocalSection(coarse, &localCoarse);
3322:   DMGetGlobalSection(coarse, &globalCoarse);

3324:   { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3325:     MPI_Datatype threeInt;
3326:     PetscMPIInt  rank;
3327:     PetscInt(*parentNodeAndIdCoarse)[3];
3328:     PetscInt(*parentNodeAndIdFine)[3];
3329:     PetscInt           p, nleaves, nleavesToParents;
3330:     PetscSF            pointSF, sfToParents;
3331:     const PetscInt    *ilocal;
3332:     const PetscSFNode *iremote;
3333:     PetscSFNode       *iremoteToParents;
3334:     PetscInt          *ilocalToParents;

3336:     MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank);
3337:     MPI_Type_contiguous(3, MPIU_INT, &threeInt);
3338:     MPI_Type_commit(&threeInt);
3339:     PetscMalloc2(pEndC - pStartC, &parentNodeAndIdCoarse, pEndF - pStartF, &parentNodeAndIdFine);
3340:     DMGetPointSF(coarse, &pointSF);
3341:     PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, &iremote);
3342:     for (p = pStartC; p < pEndC; p++) {
3343:       PetscInt parent, childId;
3344:       DMPlexGetTreeParent(coarse, p, &parent, &childId);
3345:       parentNodeAndIdCoarse[p - pStartC][0] = rank;
3346:       parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3347:       parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3348:       if (nleaves > 0) {
3349:         PetscInt leaf = -1;

3351:         if (ilocal) {
3352:           PetscFindInt(parent, nleaves, ilocal, &leaf);
3353:         } else {
3354:           leaf = p - pStartC;
3355:         }
3356:         if (leaf >= 0) {
3357:           parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3358:           parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3359:         }
3360:       }
3361:     }
3362:     for (p = pStartF; p < pEndF; p++) {
3363:       parentNodeAndIdFine[p - pStartF][0] = -1;
3364:       parentNodeAndIdFine[p - pStartF][1] = -1;
3365:       parentNodeAndIdFine[p - pStartF][2] = -1;
3366:     }
3367:     PetscSFBcastBegin(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE);
3368:     PetscSFBcastEnd(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE);
3369:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3370:       PetscInt dof;

3372:       PetscSectionGetDof(leafIndicesSec, p, &dof);
3373:       if (dof) {
3374:         PetscInt off;

3376:         PetscSectionGetOffset(leafIndicesSec, p, &off);
3377:         if (gatheredIndices) {
3378:           leafInds[off] = PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3379:         } else if (gatheredValues) {
3380:           leafVals[off] = (PetscScalar)PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3381:         }
3382:       }
3383:       if (parentNodeAndIdFine[p - pStartF][0] >= 0) nleavesToParents++;
3384:     }
3385:     PetscMalloc1(nleavesToParents, &ilocalToParents);
3386:     PetscMalloc1(nleavesToParents, &iremoteToParents);
3387:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3388:       if (parentNodeAndIdFine[p - pStartF][0] >= 0) {
3389:         ilocalToParents[nleavesToParents]        = p - pStartF;
3390:         iremoteToParents[nleavesToParents].rank  = parentNodeAndIdFine[p - pStartF][0];
3391:         iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p - pStartF][1];
3392:         nleavesToParents++;
3393:       }
3394:     }
3395:     PetscSFCreate(PetscObjectComm((PetscObject)coarse), &sfToParents);
3396:     PetscSFSetGraph(sfToParents, pEndC - pStartC, nleavesToParents, ilocalToParents, PETSC_OWN_POINTER, iremoteToParents, PETSC_OWN_POINTER);
3397:     PetscSFDestroy(&coarseToFineEmbedded);

3399:     coarseToFineEmbedded = sfToParents;

3401:     PetscFree2(parentNodeAndIdCoarse, parentNodeAndIdFine);
3402:     MPI_Type_free(&threeInt);
3403:   }

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

3409:     for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3410:       PetscSectionGetDof(globalCoarse, p, &dof);
3411:       PetscSectionGetConstraintDof(globalCoarse, p, &cdof);
3412:       if ((dof - cdof) > 0) numPointsWithDofs++;
3413:     }
3414:     PetscMalloc1(numPointsWithDofs, &pointsWithDofs);
3415:     for (p = pStartC, offset = 0; p < pEndC; p++) {
3416:       PetscSectionGetDof(globalCoarse, p, &dof);
3417:       PetscSectionGetConstraintDof(globalCoarse, p, &cdof);
3418:       if ((dof - cdof) > 0) pointsWithDofs[offset++] = p - pStartC;
3419:     }
3420:     PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);
3421:     PetscSFDestroy(&coarseToFineEmbedded);
3422:     PetscFree(pointsWithDofs);
3423:     coarseToFineEmbedded = sfDofsOnly;
3424:   }

3426:   /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3427:   PetscSFComputeDegreeBegin(coarseToFineEmbedded, &rootDegrees);
3428:   PetscSFComputeDegreeEnd(coarseToFineEmbedded, &rootDegrees);
3429:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &multiRootSec);
3430:   PetscSectionSetChart(multiRootSec, pStartC, pEndC);
3431:   for (p = pStartC; p < pEndC; p++) PetscSectionSetDof(multiRootSec, p, rootDegrees[p - pStartC]);
3432:   PetscSectionSetUp(multiRootSec);
3433:   PetscSectionGetStorageSize(multiRootSec, &numMulti);
3434:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec);
3435:   { /* distribute the leaf section */
3436:     PetscSF   multi, multiInv, indicesSF;
3437:     PetscInt *remoteOffsets, numRootIndices;

3439:     PetscSFGetMultiSF(coarseToFineEmbedded, &multi);
3440:     PetscSFCreateInverseSF(multi, &multiInv);
3441:     PetscSFDistributeSection(multiInv, leafIndicesSec, &remoteOffsets, rootIndicesSec);
3442:     PetscSFCreateSectionSF(multiInv, leafIndicesSec, remoteOffsets, rootIndicesSec, &indicesSF);
3443:     PetscFree(remoteOffsets);
3444:     PetscSFDestroy(&multiInv);
3445:     PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices);
3446:     if (gatheredIndices) {
3447:       PetscMalloc1(numRootIndices, &rootInds);
3448:       PetscSFBcastBegin(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE);
3449:       PetscSFBcastEnd(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE);
3450:     }
3451:     if (gatheredValues) {
3452:       PetscMalloc1(numRootIndices, &rootVals);
3453:       PetscSFBcastBegin(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE);
3454:       PetscSFBcastEnd(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE);
3455:     }
3456:     PetscSFDestroy(&indicesSF);
3457:   }
3458:   PetscSectionDestroy(&leafIndicesSec);
3459:   PetscFree(leafInds);
3460:   PetscFree(leafVals);
3461:   PetscSFDestroy(&coarseToFineEmbedded);
3462:   *rootMultiSec = multiRootSec;
3463:   *multiLeafSec = rootIndicesSec;
3464:   if (gatheredIndices) *gatheredIndices = rootInds;
3465:   if (gatheredValues) *gatheredValues = rootVals;
3466:   return 0;
3467: }

3469: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3470: {
3471:   DM             refTree;
3472:   PetscSection   multiRootSec, rootIndicesSec;
3473:   PetscSection   globalCoarse, globalFine;
3474:   PetscSection   localCoarse, localFine;
3475:   PetscSection   cSecRef;
3476:   PetscInt      *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd;
3477:   Mat            injRef;
3478:   PetscInt       numFields, maxDof;
3479:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
3480:   PetscInt      *offsets, *offsetsCopy, *rowOffsets;
3481:   PetscLayout    rowMap, colMap;
3482:   PetscInt       rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3483:   PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */


3486:   /* get the templates for the fine-to-coarse injection from the reference tree */
3487:   DMPlexGetReferenceTree(coarse, &refTree);
3488:   DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL);
3489:   PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd);
3490:   DMPlexReferenceTreeGetInjector(refTree, &injRef);

3492:   DMPlexGetChart(fine, &pStartF, &pEndF);
3493:   DMGetLocalSection(fine, &localFine);
3494:   DMGetGlobalSection(fine, &globalFine);
3495:   PetscSectionGetNumFields(localFine, &numFields);
3496:   DMPlexGetChart(coarse, &pStartC, &pEndC);
3497:   DMGetLocalSection(coarse, &localCoarse);
3498:   DMGetGlobalSection(coarse, &globalCoarse);
3499:   PetscSectionGetMaxDof(localCoarse, &maxDof);
3500:   {
3501:     PetscInt maxFields = PetscMax(1, numFields) + 1;
3502:     PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets);
3503:   }

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

3507:   PetscMalloc1(maxDof, &parentIndices);

3509:   /* count indices */
3510:   MatGetLayouts(mat, &rowMap, &colMap);
3511:   PetscLayoutSetUp(rowMap);
3512:   PetscLayoutSetUp(colMap);
3513:   PetscLayoutGetRange(rowMap, &rowStart, &rowEnd);
3514:   PetscLayoutGetRange(colMap, &colStart, &colEnd);
3515:   PetscCalloc2(rowEnd - rowStart, &nnzD, rowEnd - rowStart, &nnzO);
3516:   for (p = pStartC; p < pEndC; p++) {
3517:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3519:     PetscSectionGetDof(globalCoarse, p, &dof);
3520:     PetscSectionGetConstraintDof(globalCoarse, p, &cdof);
3521:     if ((dof - cdof) <= 0) continue;
3522:     PetscSectionGetOffset(globalCoarse, p, &gOff);

3524:     rowOffsets[0]  = 0;
3525:     offsetsCopy[0] = 0;
3526:     if (numFields) {
3527:       PetscInt f;

3529:       for (f = 0; f < numFields; f++) {
3530:         PetscInt fDof;
3531:         PetscSectionGetFieldDof(localCoarse, p, f, &fDof);
3532:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3533:       }
3534:       DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices);
3535:     } else {
3536:       DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices);
3537:       rowOffsets[1] = offsetsCopy[0];
3538:     }

3540:     PetscSectionGetDof(multiRootSec, p, &numLeaves);
3541:     PetscSectionGetOffset(multiRootSec, p, &leafStart);
3542:     leafEnd = leafStart + numLeaves;
3543:     for (l = leafStart; l < leafEnd; l++) {
3544:       PetscInt        numIndices, childId, offset;
3545:       const PetscInt *childIndices;

3547:       PetscSectionGetDof(rootIndicesSec, l, &numIndices);
3548:       PetscSectionGetOffset(rootIndicesSec, l, &offset);
3549:       childId      = rootIndices[offset++];
3550:       childIndices = &rootIndices[offset];
3551:       numIndices--;

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

3556:         for (i = 0; i < numIndices; i++) {
3557:           PetscInt colIndex = childIndices[i];
3558:           PetscInt rowIndex = parentIndices[i];
3559:           if (rowIndex < 0) continue;
3561:           if (colIndex >= colStart && colIndex < colEnd) {
3562:             nnzD[rowIndex - rowStart] = 1;
3563:           } else {
3564:             nnzO[rowIndex - rowStart] = 1;
3565:           }
3566:         }
3567:       } else {
3568:         PetscInt parentId, f, lim;

3570:         DMPlexGetTreeParent(refTree, childId, &parentId, NULL);

3572:         lim        = PetscMax(1, numFields);
3573:         offsets[0] = 0;
3574:         if (numFields) {
3575:           PetscInt f;

3577:           for (f = 0; f < numFields; f++) {
3578:             PetscInt fDof;
3579:             PetscSectionGetFieldDof(cSecRef, childId, f, &fDof);

3581:             offsets[f + 1] = fDof + offsets[f];
3582:           }
3583:         } else {
3584:           PetscInt cDof;

3586:           PetscSectionGetDof(cSecRef, childId, &cDof);
3587:           offsets[1] = cDof;
3588:         }
3589:         for (f = 0; f < lim; f++) {
3590:           PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3591:           PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3592:           PetscInt i, numD = 0, numO = 0;

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

3597:             if (colIndex < 0) continue;
3598:             if (colIndex >= colStart && colIndex < colEnd) {
3599:               numD++;
3600:             } else {
3601:               numO++;
3602:             }
3603:           }
3604:           for (i = parentStart; i < parentEnd; i++) {
3605:             PetscInt rowIndex = parentIndices[i];

3607:             if (rowIndex < 0) continue;
3608:             nnzD[rowIndex - rowStart] += numD;
3609:             nnzO[rowIndex - rowStart] += numO;
3610:           }
3611:         }
3612:       }
3613:     }
3614:   }
3615:   /* preallocate */
3616:   MatXAIJSetPreallocation(mat, 1, nnzD, nnzO, NULL, NULL);
3617:   PetscFree2(nnzD, nnzO);
3618:   /* insert values */
3619:   DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats);
3620:   for (p = pStartC; p < pEndC; p++) {
3621:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3623:     PetscSectionGetDof(globalCoarse, p, &dof);
3624:     PetscSectionGetConstraintDof(globalCoarse, p, &cdof);
3625:     if ((dof - cdof) <= 0) continue;
3626:     PetscSectionGetOffset(globalCoarse, p, &gOff);

3628:     rowOffsets[0]  = 0;
3629:     offsetsCopy[0] = 0;
3630:     if (numFields) {
3631:       PetscInt f;

3633:       for (f = 0; f < numFields; f++) {
3634:         PetscInt fDof;
3635:         PetscSectionGetFieldDof(localCoarse, p, f, &fDof);
3636:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3637:       }
3638:       DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices);
3639:     } else {
3640:       DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices);
3641:       rowOffsets[1] = offsetsCopy[0];
3642:     }

3644:     PetscSectionGetDof(multiRootSec, p, &numLeaves);
3645:     PetscSectionGetOffset(multiRootSec, p, &leafStart);
3646:     leafEnd = leafStart + numLeaves;
3647:     for (l = leafStart; l < leafEnd; l++) {
3648:       PetscInt        numIndices, childId, offset;
3649:       const PetscInt *childIndices;

3651:       PetscSectionGetDof(rootIndicesSec, l, &numIndices);
3652:       PetscSectionGetOffset(rootIndicesSec, l, &offset);
3653:       childId      = rootIndices[offset++];
3654:       childIndices = &rootIndices[offset];
3655:       numIndices--;

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

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

3664:         DMPlexGetTreeParent(refTree, childId, &parentId, NULL);

3666:         lim        = PetscMax(1, numFields);
3667:         offsets[0] = 0;
3668:         if (numFields) {
3669:           PetscInt f;

3671:           for (f = 0; f < numFields; f++) {
3672:             PetscInt fDof;
3673:             PetscSectionGetFieldDof(cSecRef, childId, f, &fDof);

3675:             offsets[f + 1] = fDof + offsets[f];
3676:           }
3677:         } else {
3678:           PetscInt cDof;

3680:           PetscSectionGetDof(cSecRef, childId, &cDof);
3681:           offsets[1] = cDof;
3682:         }
3683:         for (f = 0; f < lim; f++) {
3684:           PetscScalar    *childMat   = &childrenMats[childId - pRefStart][f][0];
3685:           PetscInt       *rowIndices = &parentIndices[rowOffsets[f]];
3686:           const PetscInt *colIndices = &childIndices[offsets[f]];

3688:           MatSetValues(mat, rowOffsets[f + 1] - rowOffsets[f], rowIndices, offsets[f + 1] - offsets[f], colIndices, childMat, INSERT_VALUES);
3689:         }
3690:       }
3691:     }
3692:   }
3693:   PetscSectionDestroy(&multiRootSec);
3694:   PetscSectionDestroy(&rootIndicesSec);
3695:   PetscFree(parentIndices);
3696:   DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats);
3697:   PetscFree(rootIndices);
3698:   PetscFree3(offsets, offsetsCopy, rowOffsets);

3700:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
3701:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
3702:   return 0;
3703: }

3705: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3706: {
3707:   PetscSF            coarseToFineEmbedded;
3708:   PetscSection       globalCoarse, globalFine;
3709:   PetscSection       localCoarse, localFine;
3710:   PetscSection       aSec, cSec;
3711:   PetscSection       rootValuesSec;
3712:   PetscSection       leafValuesSec;
3713:   PetscScalar       *rootValues, *leafValues;
3714:   IS                 aIS;
3715:   const PetscInt    *anchors;
3716:   Mat                cMat;
3717:   PetscInt           numFields;
3718:   PetscInt           pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd;
3719:   PetscInt           aStart, aEnd, cStart, cEnd;
3720:   PetscInt          *maxChildIds;
3721:   PetscInt          *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3722:   PetscFV            fv = NULL;
3723:   PetscInt           dim, numFVcomps = -1, fvField = -1;
3724:   DM                 cellDM = NULL, gradDM = NULL;
3725:   const PetscScalar *cellGeomArray = NULL;
3726:   const PetscScalar *gradArray     = NULL;

3728:   VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
3729:   DMPlexGetChart(coarse, &pStartC, &pEndC);
3730:   DMPlexGetSimplexOrBoxCells(coarse, 0, &cellStart, &cellEnd);
3731:   DMPlexGetChart(fine, &pStartF, &pEndF);
3732:   DMGetGlobalSection(fine, &globalFine);
3733:   DMGetCoordinateDim(coarse, &dim);
3734:   { /* winnow fine points that don't have global dofs out of the sf */
3735:     PetscInt        nleaves, l;
3736:     const PetscInt *leaves;
3737:     PetscInt        dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;

3739:     PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL);

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

3744:       PetscSectionGetDof(globalFine, p, &dof);
3745:       PetscSectionGetConstraintDof(globalFine, p, &cdof);
3746:       if ((dof - cdof) > 0) numPointsWithDofs++;
3747:     }
3748:     PetscMalloc1(numPointsWithDofs, &pointsWithDofs);
3749:     for (l = 0, offset = 0; l < nleaves; l++) {
3750:       PetscInt p = leaves ? leaves[l] : l;

3752:       PetscSectionGetDof(globalFine, p, &dof);
3753:       PetscSectionGetConstraintDof(globalFine, p, &cdof);
3754:       if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
3755:     }
3756:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
3757:     PetscFree(pointsWithDofs);
3758:   }
3759:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
3760:   PetscMalloc1(pEndC - pStartC, &maxChildIds);
3761:   for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
3762:   PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX);
3763:   PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX);

3765:   DMGetLocalSection(coarse, &localCoarse);
3766:   DMGetGlobalSection(coarse, &globalCoarse);

3768:   DMPlexGetAnchors(coarse, &aSec, &aIS);
3769:   ISGetIndices(aIS, &anchors);
3770:   PetscSectionGetChart(aSec, &aStart, &aEnd);

3772:   DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL);
3773:   PetscSectionGetChart(cSec, &cStart, &cEnd);

3775:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
3776:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootValuesSec);
3777:   PetscSectionSetChart(rootValuesSec, pStartC, pEndC);
3778:   PetscSectionGetNumFields(localCoarse, &numFields);
3779:   {
3780:     PetscInt maxFields = PetscMax(1, numFields) + 1;
3781:     PetscMalloc7(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &newOffsets, maxFields, &newOffsetsCopy, maxFields, &rowOffsets, maxFields, &numD, maxFields, &numO);
3782:   }
3783:   if (grad) {
3784:     PetscInt i;

3786:     VecGetDM(cellGeom, &cellDM);
3787:     VecGetArrayRead(cellGeom, &cellGeomArray);
3788:     VecGetDM(grad, &gradDM);
3789:     VecGetArrayRead(grad, &gradArray);
3790:     for (i = 0; i < PetscMax(1, numFields); i++) {
3791:       PetscObject  obj;
3792:       PetscClassId id;

3794:       DMGetField(coarse, i, NULL, &obj);
3795:       PetscObjectGetClassId(obj, &id);
3796:       if (id == PETSCFV_CLASSID) {
3797:         fv = (PetscFV)obj;
3798:         PetscFVGetNumComponents(fv, &numFVcomps);
3799:         fvField = i;
3800:         break;
3801:       }
3802:     }
3803:   }

3805:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
3806:     PetscInt dof;
3807:     PetscInt maxChildId = maxChildIds[p - pStartC];
3808:     PetscInt numValues  = 0;

3810:     PetscSectionGetDof(globalCoarse, p, &dof);
3811:     if (dof < 0) dof = -(dof + 1);
3812:     offsets[0]    = 0;
3813:     newOffsets[0] = 0;
3814:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
3815:       PetscInt *closure = NULL, closureSize, cl;

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

3821:         PetscSectionGetDof(localCoarse, c, &clDof);
3822:         numValues += clDof;
3823:       }
3824:       DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
3825:     } else if (maxChildId == -1) {
3826:       PetscSectionGetDof(localCoarse, p, &numValues);
3827:     }
3828:     /* we will pack the column indices with the field offsets */
3829:     if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
3830:       /* also send the centroid, and the gradient */
3831:       numValues += dim * (1 + numFVcomps);
3832:     }
3833:     PetscSectionSetDof(rootValuesSec, p, numValues);
3834:   }
3835:   PetscSectionSetUp(rootValuesSec);
3836:   {
3837:     PetscInt           numRootValues;
3838:     const PetscScalar *coarseArray;

3840:     PetscSectionGetStorageSize(rootValuesSec, &numRootValues);
3841:     PetscMalloc1(numRootValues, &rootValues);
3842:     VecGetArrayRead(vecCoarseLocal, &coarseArray);
3843:     for (p = pStartC; p < pEndC; p++) {
3844:       PetscInt     numValues;
3845:       PetscInt     pValOff;
3846:       PetscScalar *pVal;
3847:       PetscInt     maxChildId = maxChildIds[p - pStartC];

3849:       PetscSectionGetDof(rootValuesSec, p, &numValues);
3850:       if (!numValues) continue;
3851:       PetscSectionGetOffset(rootValuesSec, p, &pValOff);
3852:       pVal = &(rootValues[pValOff]);
3853:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
3854:         PetscInt closureSize = numValues;
3855:         DMPlexVecGetClosure(coarse, NULL, vecCoarseLocal, p, &closureSize, &pVal);
3856:         if (grad && p >= cellStart && p < cellEnd) {
3857:           PetscFVCellGeom *cg;
3858:           PetscScalar     *gradVals = NULL;
3859:           PetscInt         i;

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

3863:           DMPlexPointLocalRead(cellDM, p, cellGeomArray, (void *)&cg);
3864:           for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
3865:           pVal += dim;
3866:           DMPlexPointGlobalRead(gradDM, p, gradArray, (void *)&gradVals);
3867:           for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
3868:         }
3869:       } else if (maxChildId == -1) {
3870:         PetscInt lDof, lOff, i;

3872:         PetscSectionGetDof(localCoarse, p, &lDof);
3873:         PetscSectionGetOffset(localCoarse, p, &lOff);
3874:         for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
3875:       }
3876:     }
3877:     VecRestoreArrayRead(vecCoarseLocal, &coarseArray);
3878:     PetscFree(maxChildIds);
3879:   }
3880:   {
3881:     PetscSF   valuesSF;
3882:     PetscInt *remoteOffsetsValues, numLeafValues;

3884:     PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafValuesSec);
3885:     PetscSFDistributeSection(coarseToFineEmbedded, rootValuesSec, &remoteOffsetsValues, leafValuesSec);
3886:     PetscSFCreateSectionSF(coarseToFineEmbedded, rootValuesSec, remoteOffsetsValues, leafValuesSec, &valuesSF);
3887:     PetscSFDestroy(&coarseToFineEmbedded);
3888:     PetscFree(remoteOffsetsValues);
3889:     PetscSectionGetStorageSize(leafValuesSec, &numLeafValues);
3890:     PetscMalloc1(numLeafValues, &leafValues);
3891:     PetscSFBcastBegin(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE);
3892:     PetscSFBcastEnd(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE);
3893:     PetscSFDestroy(&valuesSF);
3894:     PetscFree(rootValues);
3895:     PetscSectionDestroy(&rootValuesSec);
3896:   }
3897:   DMGetLocalSection(fine, &localFine);
3898:   {
3899:     PetscInt       maxDof;
3900:     PetscInt      *rowIndices;
3901:     DM             refTree;
3902:     PetscInt     **refPointFieldN;
3903:     PetscScalar ***refPointFieldMats;
3904:     PetscSection   refConSec, refAnSec;
3905:     PetscInt       pRefStart, pRefEnd, leafStart, leafEnd;
3906:     PetscScalar   *pointWork;

3908:     PetscSectionGetMaxDof(localFine, &maxDof);
3909:     DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices);
3910:     DMGetWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork);
3911:     DMPlexGetReferenceTree(fine, &refTree);
3912:     DMCopyDisc(fine, refTree);
3913:     DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN);
3914:     DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL);
3915:     DMPlexGetAnchors(refTree, &refAnSec, NULL);
3916:     PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd);
3917:     PetscSectionGetChart(leafValuesSec, &leafStart, &leafEnd);
3918:     DMPlexGetSimplexOrBoxCells(fine, 0, &cellStart, &cellEnd);
3919:     for (p = leafStart; p < leafEnd; p++) {
3920:       PetscInt           gDof, gcDof, gOff, lDof;
3921:       PetscInt           numValues, pValOff;
3922:       PetscInt           childId;
3923:       const PetscScalar *pVal;
3924:       const PetscScalar *fvGradData = NULL;

3926:       PetscSectionGetDof(globalFine, p, &gDof);
3927:       PetscSectionGetDof(localFine, p, &lDof);
3928:       PetscSectionGetConstraintDof(globalFine, p, &gcDof);
3929:       if ((gDof - gcDof) <= 0) continue;
3930:       PetscSectionGetOffset(globalFine, p, &gOff);
3931:       PetscSectionGetDof(leafValuesSec, p, &numValues);
3932:       if (!numValues) continue;
3933:       PetscSectionGetOffset(leafValuesSec, p, &pValOff);
3934:       pVal              = &leafValues[pValOff];
3935:       offsets[0]        = 0;
3936:       offsetsCopy[0]    = 0;
3937:       newOffsets[0]     = 0;
3938:       newOffsetsCopy[0] = 0;
3939:       childId           = cids[p - pStartF];
3940:       if (numFields) {
3941:         PetscInt f;
3942:         for (f = 0; f < numFields; f++) {
3943:           PetscInt rowDof;

3945:           PetscSectionGetFieldDof(localFine, p, f, &rowDof);
3946:           offsets[f + 1]     = offsets[f] + rowDof;
3947:           offsetsCopy[f + 1] = offsets[f + 1];
3948:           /* TODO: closure indices */
3949:           newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
3950:         }
3951:         DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices);
3952:       } else {
3953:         offsets[0]    = 0;
3954:         offsets[1]    = lDof;
3955:         newOffsets[0] = 0;
3956:         newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
3957:         DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices);
3958:       }
3959:       if (childId == -1) { /* no child interpolation: one nnz per */
3960:         VecSetValues(vecFine, numValues, rowIndices, pVal, INSERT_VALUES);
3961:       } else {
3962:         PetscInt f;

3964:         if (grad && p >= cellStart && p < cellEnd) {
3965:           numValues -= (dim * (1 + numFVcomps));
3966:           fvGradData = &pVal[numValues];
3967:         }
3968:         for (f = 0; f < PetscMax(1, numFields); f++) {
3969:           const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
3970:           PetscInt           numRows  = offsets[f + 1] - offsets[f];
3971:           PetscInt           numCols  = newOffsets[f + 1] - newOffsets[f];
3972:           const PetscScalar *cVal     = &pVal[newOffsets[f]];
3973:           PetscScalar       *rVal     = &pointWork[offsets[f]];
3974:           PetscInt           i, j;

3976: #if 0
3977:           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);
3978: #endif
3979:           for (i = 0; i < numRows; i++) {
3980:             PetscScalar val = 0.;
3981:             for (j = 0; j < numCols; j++) val += childMat[i * numCols + j] * cVal[j];
3982:             rVal[i] = val;
3983:           }
3984:           if (f == fvField && p >= cellStart && p < cellEnd) {
3985:             PetscReal          centroid[3];
3986:             PetscScalar        diff[3];
3987:             const PetscScalar *parentCentroid = &fvGradData[0];
3988:             const PetscScalar *gradient       = &fvGradData[dim];

3990:             DMPlexComputeCellGeometryFVM(fine, p, NULL, centroid, NULL);
3991:             for (i = 0; i < dim; i++) diff[i] = centroid[i] - parentCentroid[i];
3992:             for (i = 0; i < numFVcomps; i++) {
3993:               PetscScalar val = 0.;

3995:               for (j = 0; j < dim; j++) val += gradient[dim * i + j] * diff[j];
3996:               rVal[i] += val;
3997:             }
3998:           }
3999:           VecSetValues(vecFine, numRows, &rowIndices[offsets[f]], rVal, INSERT_VALUES);
4000:         }
4001:       }
4002:     }
4003:     DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN);
4004:     DMRestoreWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork);
4005:     DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices);
4006:   }
4007:   PetscFree(leafValues);
4008:   PetscSectionDestroy(&leafValuesSec);
4009:   PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO);
4010:   ISRestoreIndices(aIS, &anchors);
4011:   return 0;
4012: }

4014: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4015: {
4016:   DM             refTree;
4017:   PetscSection   multiRootSec, rootIndicesSec;
4018:   PetscSection   globalCoarse, globalFine;
4019:   PetscSection   localCoarse, localFine;
4020:   PetscSection   cSecRef;
4021:   PetscInt      *parentIndices, pRefStart, pRefEnd;
4022:   PetscScalar   *rootValues, *parentValues;
4023:   Mat            injRef;
4024:   PetscInt       numFields, maxDof;
4025:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
4026:   PetscInt      *offsets, *offsetsCopy, *rowOffsets;
4027:   PetscLayout    rowMap, colMap;
4028:   PetscInt       rowStart, rowEnd, colStart, colEnd;
4029:   PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */


4032:   /* get the templates for the fine-to-coarse injection from the reference tree */
4033:   VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
4034:   VecSetOption(vecCoarse, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
4035:   DMPlexGetReferenceTree(coarse, &refTree);
4036:   DMCopyDisc(coarse, refTree);
4037:   DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL);
4038:   PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd);
4039:   DMPlexReferenceTreeGetInjector(refTree, &injRef);

4041:   DMPlexGetChart(fine, &pStartF, &pEndF);
4042:   DMGetLocalSection(fine, &localFine);
4043:   DMGetGlobalSection(fine, &globalFine);
4044:   PetscSectionGetNumFields(localFine, &numFields);
4045:   DMPlexGetChart(coarse, &pStartC, &pEndC);
4046:   DMGetLocalSection(coarse, &localCoarse);
4047:   DMGetGlobalSection(coarse, &globalCoarse);
4048:   PetscSectionGetMaxDof(localCoarse, &maxDof);
4049:   {
4050:     PetscInt maxFields = PetscMax(1, numFields) + 1;
4051:     PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets);
4052:   }

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

4056:   PetscMalloc2(maxDof, &parentIndices, maxDof, &parentValues);

4058:   /* count indices */
4059:   VecGetLayout(vecFine, &colMap);
4060:   VecGetLayout(vecCoarse, &rowMap);
4061:   PetscLayoutSetUp(rowMap);
4062:   PetscLayoutSetUp(colMap);
4063:   PetscLayoutGetRange(rowMap, &rowStart, &rowEnd);
4064:   PetscLayoutGetRange(colMap, &colStart, &colEnd);
4065:   /* insert values */
4066:   DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats);
4067:   for (p = pStartC; p < pEndC; p++) {
4068:     PetscInt  numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4069:     PetscBool contribute = PETSC_FALSE;

4071:     PetscSectionGetDof(globalCoarse, p, &dof);
4072:     PetscSectionGetConstraintDof(globalCoarse, p, &cdof);
4073:     if ((dof - cdof) <= 0) continue;
4074:     PetscSectionGetDof(localCoarse, p, &dof);
4075:     PetscSectionGetOffset(globalCoarse, p, &gOff);

4077:     rowOffsets[0]  = 0;
4078:     offsetsCopy[0] = 0;
4079:     if (numFields) {
4080:       PetscInt f;

4082:       for (f = 0; f < numFields; f++) {
4083:         PetscInt fDof;
4084:         PetscSectionGetFieldDof(localCoarse, p, f, &fDof);
4085:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
4086:       }
4087:       DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices);
4088:     } else {
4089:       DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices);
4090:       rowOffsets[1] = offsetsCopy[0];
4091:     }

4093:     PetscSectionGetDof(multiRootSec, p, &numLeaves);
4094:     PetscSectionGetOffset(multiRootSec, p, &leafStart);
4095:     leafEnd = leafStart + numLeaves;
4096:     for (l = 0; l < dof; l++) parentValues[l] = 0.;
4097:     for (l = leafStart; l < leafEnd; l++) {
4098:       PetscInt           numIndices, childId, offset;
4099:       const PetscScalar *childValues;

4101:       PetscSectionGetDof(rootIndicesSec, l, &numIndices);
4102:       PetscSectionGetOffset(rootIndicesSec, l, &offset);
4103:       childId     = (PetscInt)PetscRealPart(rootValues[offset++]);
4104:       childValues = &rootValues[offset];
4105:       numIndices--;

4107:       if (childId == -2) { /* skip */
4108:         continue;
4109:       } else if (childId == -1) { /* equivalent points: scatter */
4110:         PetscInt m;

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

4117:         contribute = PETSC_TRUE;
4118:         DMPlexGetTreeParent(refTree, childId, &parentId, NULL);

4120:         lim        = PetscMax(1, numFields);
4121:         offsets[0] = 0;
4122:         if (numFields) {
4123:           PetscInt f;

4125:           for (f = 0; f < numFields; f++) {
4126:             PetscInt fDof;
4127:             PetscSectionGetFieldDof(cSecRef, childId, f, &fDof);

4129:             offsets[f + 1] = fDof + offsets[f];
4130:           }
4131:         } else {
4132:           PetscInt cDof;

4134:           PetscSectionGetDof(cSecRef, childId, &cDof);
4135:           offsets[1] = cDof;
4136:         }
4137:         for (f = 0; f < lim; f++) {
4138:           PetscScalar       *childMat = &childrenMats[childId - pRefStart][f][0];
4139:           PetscInt           n        = offsets[f + 1] - offsets[f];
4140:           PetscInt           m        = rowOffsets[f + 1] - rowOffsets[f];
4141:           PetscInt           i, j;
4142:           const PetscScalar *colValues = &childValues[offsets[f]];

4144:           for (i = 0; i < m; i++) {
4145:             PetscScalar val = 0.;
4146:             for (j = 0; j < n; j++) val += childMat[n * i + j] * colValues[j];
4147:             parentValues[rowOffsets[f] + i] += val;
4148:           }
4149:         }
4150:       }
4151:     }
4152:     if (contribute) VecSetValues(vecCoarse, dof, parentIndices, parentValues, INSERT_VALUES);
4153:   }
4154:   PetscSectionDestroy(&multiRootSec);
4155:   PetscSectionDestroy(&rootIndicesSec);
4156:   PetscFree2(parentIndices, parentValues);
4157:   DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats);
4158:   PetscFree(rootValues);
4159:   PetscFree3(offsets, offsetsCopy, rowOffsets);
4160:   return 0;
4161: }

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

4168:   collective

4170:   Input Parameters:
4171: + dmIn        - The DMPlex mesh for the input vector
4172: . vecIn       - The input vector
4173: . sfRefine    - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in
4174:                 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn
4175: . sfCoarsen   - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in
4176:                 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn
4177: . cidsRefine  - The childIds of the points in dmOut.  These childIds relate back to the reference tree: childid[j] = k implies
4178:                 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference
4179:                 tree was refined from its parent.  childid[j] = -1 indicates that the point j in dmOut is exactly
4180:                 equivalent to its root in dmIn, so no interpolation is necessary.  childid[j] = -2 indicates that this
4181:                 point j in dmOut is not a leaf of sfRefine.
4182: . cidsCoarsen - The childIds of the points in dmIn.  These childIds relate back to the reference tree: childid[j] = k implies
4183:                 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference
4184:                 tree coarsens to its parent.  childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen.
4185: . useBCs      - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer.
4186: - time        - Used if boundary values are time dependent.

4188:   Output Parameters:
4189: . vecOut      - Using interpolation and injection operators calculated on the reference tree, the transferred
4190:                 projection of vecIn from dmIn to dmOut.  Note that any field discretized with a PetscFV finite volume
4191:                 method that uses gradient reconstruction will use reconstructed gradients when interpolating from
4192:                 coarse points to fine points.

4194:   Level: developer

4196: .seealso: `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()`
4197: @*/
4198: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4199: {
4200:   VecSet(vecOut, 0.0);
4201:   if (sfRefine) {
4202:     Vec vecInLocal;
4203:     DM  dmGrad   = NULL;
4204:     Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;

4206:     DMGetLocalVector(dmIn, &vecInLocal);
4207:     VecSet(vecInLocal, 0.0);
4208:     {
4209:       PetscInt numFields, i;

4211:       DMGetNumFields(dmIn, &numFields);
4212:       for (i = 0; i < numFields; i++) {
4213:         PetscObject  obj;
4214:         PetscClassId classid;

4216:         DMGetField(dmIn, i, NULL, &obj);
4217:         PetscObjectGetClassId(obj, &classid);
4218:         if (classid == PETSCFV_CLASSID) {
4219:           DMPlexGetDataFVM(dmIn, (PetscFV)obj, &cellGeom, &faceGeom, &dmGrad);
4220:           break;
4221:         }
4222:       }
4223:     }
4224:     if (useBCs) DMPlexInsertBoundaryValues(dmIn, PETSC_TRUE, vecInLocal, time, faceGeom, cellGeom, NULL);
4225:     DMGlobalToLocalBegin(dmIn, vecIn, INSERT_VALUES, vecInLocal);
4226:     DMGlobalToLocalEnd(dmIn, vecIn, INSERT_VALUES, vecInLocal);
4227:     if (dmGrad) {
4228:       DMGetGlobalVector(dmGrad, &grad);
4229:       DMPlexReconstructGradientsFVM(dmIn, vecInLocal, grad);
4230:     }
4231:     DMPlexTransferVecTree_Interpolate(dmIn, vecInLocal, dmOut, vecOut, sfRefine, cidsRefine, grad, cellGeom);
4232:     DMRestoreLocalVector(dmIn, &vecInLocal);
4233:     if (dmGrad) DMRestoreGlobalVector(dmGrad, &grad);
4234:   }
4235:   if (sfCoarsen) DMPlexTransferVecTree_Inject(dmIn, vecIn, dmOut, vecOut, sfCoarsen, cidsCoarsen);
4236:   VecAssemblyBegin(vecOut);
4237:   VecAssemblyEnd(vecOut);
4238:   return 0;
4239: }