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: [](chapter_unstructured), `DM`, `DMPLEX`,`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: [](chapter_unstructured), `DM`, `DMPLEX`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `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: [](chapter_unstructured), `DM`, `DMPLEX`,`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: [](chapter_unstructured), `DM`, `DMPLEX`, `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: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeParent()`
1075: @*/
1076: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1077: {
1078:   DM_Plex     *mesh = (DM_Plex *)dm->data;
1079:   PetscSection childSec;
1080:   PetscInt     dof = 0;

1083:   childSec = mesh->childSection;
1084:   if (childSec && point >= childSec->pStart && point < childSec->pEnd) PetscSectionGetDof(childSec, point, &dof);
1085:   if (numChildren) *numChildren = dof;
1086:   if (children) {
1087:     if (dof) {
1088:       PetscInt off;

1090:       PetscSectionGetOffset(childSec, point, &off);
1091:       *children = &mesh->children[off];
1092:     } else {
1093:       *children = NULL;
1094:     }
1095:   }
1096:   return 0;
1097: }

1099: 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)
1100: {
1101:   PetscInt f, b, p, c, offset, qPoints;

1103:   PetscSpaceEvaluate(space, nPoints, points, work, NULL, NULL);
1104:   for (f = 0, offset = 0; f < nFunctionals; f++) {
1105:     qPoints = pointsPerFn[f];
1106:     for (b = 0; b < nBasis; b++) {
1107:       PetscScalar val = 0.;

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

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

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

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

1157:     PetscDSGetDiscretization(ds, f, &disc);
1158:     PetscObjectGetClassId(disc, &id);
1159:     if (id == PETSCFE_CLASSID) {
1160:       PetscFE fe = (PetscFE)disc;

1162:       PetscFEGetBasisSpace(fe, &bspace);
1163:       PetscFEGetDualSpace(fe, &dspace);
1164:       PetscDualSpaceGetDimension(dspace, &fSize);
1165:       PetscFEGetNumComponents(fe, &Nc);
1166:     } else if (id == PETSCFV_CLASSID) {
1167:       PetscFV fv = (PetscFV)disc;

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

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

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

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

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

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

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

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

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

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

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

1303:           if (numFields) {
1304:             PetscSectionGetFieldDof(section, a, f, &aSecDof);
1305:             PetscSectionGetFieldOffset(section, a, f, &aSecOff);
1306:           } else {
1307:             PetscSectionGetDof(section, a, &aSecDof);
1308:             PetscSectionGetOffset(section, a, &aSecOff);
1309:           }
1310:           if (!aSecDof) continue;

1312:           for (j = 0; j < closureSizeP; j++) {
1313:             PetscInt q  = closureP[2 * j];
1314:             PetscInt oq = closureP[2 * j + 1];

1316:             if (q == a) {
1317:               PetscInt           r, s, nWorkP;
1318:               const PetscInt    *permP;
1319:               const PetscScalar *flipP;

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

1365:   return 0;
1366: }

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

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

1397:     DMPlexGetTreeParent(refTree, p, &parent, NULL);
1398:     PetscSectionGetDof(refConSec, p, &pDof);
1399:     if (!pDof || parent == p) continue;

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

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

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

1424:         if (numFields) {
1425:           PetscSectionGetFieldDof(refSection, q, f, &aDof);
1426:           PetscSectionGetFieldOffset(refSection, q, f, &aOff);
1427:         } else {
1428:           PetscSectionGetDof(refSection, q, &aDof);
1429:           PetscSectionGetOffset(refSection, q, &aOff);
1430:         }

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

1440:         for (i = 0; i < closureSize; i++) {
1441:           PetscInt           q = closure[2 * i];
1442:           PetscInt           aDof, aOff, j;
1443:           const PetscScalar *flip = flips ? flips[i] : NULL;

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

1477: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1478: {
1479:   PetscDS        ds;
1480:   PetscInt     **refPointFieldN;
1481:   PetscScalar ***refPointFieldMats;
1482:   PetscInt       numFields, maxFields, pRefStart, pRefEnd, p, f;
1483:   PetscSection   refConSec;

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

1497:     DMPlexGetTreeParent(refTree, p, &parent, NULL);
1498:     PetscSectionGetDof(refConSec, p, &pDof);
1499:     if (!pDof || parent == p) continue;

1501:     for (f = 0; f < maxFields; f++) {
1502:       PetscInt cDof;

1504:       if (numFields) {
1505:         PetscSectionGetFieldDof(refConSec, p, f, &cDof);
1506:       } else {
1507:         PetscSectionGetDof(refConSec, p, &cDof);
1508:       }

1510:       PetscFree(refPointFieldMats[p - pRefStart][f]);
1511:     }
1512:     PetscFree(refPointFieldMats[p - pRefStart]);
1513:     PetscFree(refPointFieldN[p - pRefStart]);
1514:   }
1515:   PetscFree(refPointFieldMats);
1516:   PetscFree(refPointFieldN);
1517:   return 0;
1518: }

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

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

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

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

1561:     DMPlexGetTreeParent(dm, point, &parent, NULL);
1562:     if (parent == point) {
1563:       p++;
1564:     } else {
1565:       PetscInt size, closureSize, *closure = NULL, i;

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

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

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

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

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

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

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

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

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

1727:   /* clean up */
1728:   ISRestoreIndices(anIS, &anchors);
1729:   PetscFree2(perm, iperm);
1730:   PetscFree(pointWork);
1731:   DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN);
1732:   return 0;
1733: }

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

1749:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1750:   DMGetDimension(dm, &dim);
1751:   DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);
1752:   DMSetDimension(*ncdm, dim);

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

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

1790:     PetscMalloc1(pNewEnd[dim], &newConeSizes);
1791:     {
1792:       DMPolytopeType pct, qct;
1793:       PetscInt       kStart, kEnd, k, closureSizeK, *closureK = NULL, j;

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

1798:       for (k = kStart; k < kEnd; k++) {
1799:         perm[k - kStart]      = k;
1800:         iperm[k - kStart]     = k - kStart;
1801:         preOrient[k - kStart] = 0;
1802:       }

1804:       DMPlexGetTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK);
1805:       for (j = 1; j < closureSizeK; j++) {
1806:         PetscInt parentOrientA = closureK[2 * j + 1];
1807:         PetscInt parentOrientB = cellClosure[2 * j + 1];
1808:         PetscInt p, q;

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

1823:           DMPlexGetTreeChildren(K, p, &numChildren, &children);
1824:           for (i = 0; i < numChildren; i++) {
1825:             PetscInt kPerm, oPerm;

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

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

1856:       DMPlexGetHeightStratum(K, d, &kStart, &kEnd);
1857:       for (k = kStart; k < kEnd; k++) {
1858:         PetscInt kParent;

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

1871:     PetscSectionSetUp(parentSection);
1872:     PetscSectionGetStorageSize(parentSection, &numPointsWithParents);
1873:     PetscMalloc2(numNewCones, &newCones, numNewCones, &newOrientations);
1874:     PetscMalloc2(numPointsWithParents, &parents, numPointsWithParents, &childIDs);

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

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

1897:       DMPlexGetHeightStratum(K, d, &kStart, &kEnd);
1898:       for (k = kStart; k < kEnd; k++) {
1899:         PetscInt kPerm = perm[k], kParent;
1900:         PetscInt preO  = preOrient[k];

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

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

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

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

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

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

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

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

1992:       tmp                = pNewCount[d];
1993:       pNewCount[d]       = pNewCount[dim - d];
1994:       pNewCount[dim - d] = tmp;
1995:     }

1997:     DMPlexCreateFromDAG(*ncdm, dim, pNewCount, newConeSizes, newCones, newOrientations, newVertexCoords);
1998:     DMPlexSetReferenceTree(*ncdm, K);
1999:     DMPlexSetTree(*ncdm, parentSection, parents, childIDs);

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

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

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

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

2037:   return 0;
2038: }

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

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

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

2091:   DMGetLocalSection(coarse, &localCoarse);
2092:   DMGetGlobalSection(coarse, &globalCoarse);

2094:   DMPlexGetAnchors(coarse, &aSec, &aIS);
2095:   ISGetIndices(aIS, &anchors);
2096:   PetscSectionGetChart(aSec, &aStart, &aEnd);

2098:   DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL);
2099:   PetscSectionGetChart(cSec, &cStart, &cEnd);

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

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

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

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

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

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

2168:         PetscSectionGetOffset(aSec, p, &aOff);
2169:         for (f = 0; f < numFields; f++) {
2170:           PetscInt fDof;

2172:           PetscSectionGetFieldDof(localCoarse, p, f, &fDof);
2173:           offsets[f + 1] = fDof;
2174:         }
2175:         for (a = 0; a < aDof; a++) {
2176:           PetscInt anchor = anchors[a + aOff], aLocalDof;

2178:           PetscSectionGetDof(localCoarse, anchor, &aLocalDof);
2179:           numColIndices += aLocalDof;
2180:           for (f = 0; f < numFields; f++) {
2181:             PetscInt fDof;

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

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

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

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

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

2267:               for (f = 0; f < numFields; f++) {
2268:                 PetscInt fDof;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2828:       PetscSectionGetDof(section, child, &dof);
2829:       numChildDof += dof;
2830:     }
2831:     PetscSectionGetDof(section, p, &numSelfDof);
2832:     if (!numChildDof || !numSelfDof) continue;
2833:     for (f = 0; f < numFields; f++) {
2834:       PetscInt selfOff;

2836:       if (numSecFields) { /* count the dofs for just this field */
2837:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
2838:           PetscInt child = children[i];
2839:           PetscInt dof;

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

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

2872:       PetscSectionGetDof(section, child, &dof);
2873:       numChildDof += dof;
2874:     }
2875:     PetscSectionGetDof(section, p, &numSelfDof);
2876:     if (!numChildDof || !numSelfDof) continue;

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

2890:       if (numSecFields) {
2891:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
2892:           PetscInt child = children[i];
2893:           PetscInt dof;

2895:           PetscSectionGetFieldDof(section, child, f, &dof);
2896:           numChildDof += dof;
2897:         }
2898:         PetscSectionGetFieldDof(section, p, f, &numSelfDof);
2899:         PetscSectionGetFieldOffset(section, p, f, &selfOff);
2900:       } else {
2901:         PetscSectionGetOffset(section, p, &selfOff);
2902:       }

2904:       /* find a cell whose closure contains p */
2905:       if (p >= cStart && p < cEnd) {
2906:         parentCell = p;
2907:       } else {
2908:         PetscInt *star = NULL;
2909:         PetscInt  numStar;

2911:         parentCell = -1;
2912:         DMPlexGetTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star);
2913:         for (i = numStar - 1; i >= 0; i--) {
2914:           PetscInt c = star[2 * i];

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

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

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

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

2962:         for (i = 0; i < numChildren; i++) {
2963:           PetscInt child = children[i];
2964:           PetscInt dof, off, j;

2966:           if (numSecFields) {
2967:             PetscSectionGetFieldDof(cSection, child, f, &dof);
2968:             PetscSectionGetFieldOffset(cSection, child, f, &off);
2969:           } else {
2970:             PetscSectionGetDof(cSection, child, &dof);
2971:             PetscSectionGetOffset(cSection, child, &off);
2972:           }

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

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

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

3014:             for (d = 0; d < dim; d++) pointScalar[d] = points[dim * j + d];
3015:             point = pointScalar;
3016: #else
3017:             point = pointReal;
3018: #endif

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

3022:             for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3023:               PetscInt  child = children[k];
3024:               PetscInt *star  = NULL;
3025:               PetscInt  numStar, s;

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

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

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

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

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

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

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

3092:         DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);
3093:         for (i = 0, childCell = 0; i < numChildren; i++) {
3094:           PetscInt  child = children[i], j;
3095:           PetscReal childVol;

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

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

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

3136:     DMPlexGetTreeParent(refTree, p, &parent, NULL);
3137:     PetscSectionGetDof(refConSec, p, &pDof);
3138:     PetscSectionGetDof(refSection, parent, &parentDof);
3139:     if (!pDof || !parentDof || parent == p) continue;

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

3145:       if (numFields > 1) {
3146:         PetscSectionGetFieldDof(refConSec, p, f, &cDof);
3147:         PetscSectionGetFieldOffset(refConSec, p, f, &cOff);
3148:       } else {
3149:         PetscSectionGetDof(refConSec, p, &cDof);
3150:         PetscSectionGetOffset(refConSec, p, &cOff);
3151:       }

3153:       for (r = 0; r < cDof; r++) rows[r] = cOff + r;
3154:       numCols = 0;
3155:       {
3156:         PetscInt aDof, aOff, j;

3158:         if (numFields > 1) {
3159:           PetscSectionGetFieldDof(refSection, parent, f, &aDof);
3160:           PetscSectionGetFieldOffset(refSection, parent, f, &aOff);
3161:         } else {
3162:           PetscSectionGetDof(refSection, parent, &aDof);
3163:           PetscSectionGetOffset(refSection, parent, &aOff);
3164:         }

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

3179: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3180: {
3181:   PetscDS        ds;
3182:   PetscScalar ***refPointFieldMats;
3183:   PetscInt       numFields, pRefStart, pRefEnd, p, f;
3184:   PetscSection   refConSec, refSection;

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

3196:     DMPlexGetTreeParent(refTree, p, &parent, NULL);
3197:     PetscSectionGetDof(refConSec, p, &pDof);
3198:     PetscSectionGetDof(refSection, parent, &parentDof);
3199:     if (!pDof || !parentDof || parent == p) continue;

3201:     for (f = 0; f < numFields; f++) {
3202:       PetscInt cDof;

3204:       if (numFields > 1) {
3205:         PetscSectionGetFieldDof(refConSec, p, f, &cDof);
3206:       } else {
3207:         PetscSectionGetDof(refConSec, p, &cDof);
3208:       }

3210:       PetscFree(refPointFieldMats[p - pRefStart][f]);
3211:     }
3212:     PetscFree(refPointFieldMats[p - pRefStart]);
3213:   }
3214:   PetscFree(refPointFieldMats);
3215:   return 0;
3216: }

3218: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree, Mat *injRef)
3219: {
3220:   Mat         cMatRef;
3221:   PetscObject injRefObj;

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

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

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

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

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

3283:         pointsWithDofs[offset++] = l;

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

3287:         pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3288:         if (gatheredValues) {
3289:           PetscInt i;

3291:           pVal = &leafVals[off + 1];
3292:           for (i = 0; i < dof; i++) pVal[i] = 0.;
3293:         }
3294:         PetscSectionGetOffset(globalFine, p, &gOff);

3296:         offsets[0] = 0;
3297:         if (numFields) {
3298:           PetscInt f;

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

3316:   DMPlexGetChart(coarse, &pStartC, &pEndC);
3317:   DMGetLocalSection(coarse, &localCoarse);
3318:   DMGetGlobalSection(coarse, &globalCoarse);

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

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

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

3368:       PetscSectionGetDof(leafIndicesSec, p, &dof);
3369:       if (dof) {
3370:         PetscInt off;

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

3395:     coarseToFineEmbedded = sfToParents;

3397:     PetscFree2(parentNodeAndIdCoarse, parentNodeAndIdFine);
3398:     MPI_Type_free(&threeInt);
3399:   }

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

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

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

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

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


3482:   /* get the templates for the fine-to-coarse injection from the reference tree */
3483:   DMPlexGetReferenceTree(coarse, &refTree);
3484:   DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL);
3485:   PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd);
3486:   DMPlexReferenceTreeGetInjector(refTree, &injRef);

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

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

3503:   PetscMalloc1(maxDof, &parentIndices);

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

3515:     PetscSectionGetDof(globalCoarse, p, &dof);
3516:     PetscSectionGetConstraintDof(globalCoarse, p, &cdof);
3517:     if ((dof - cdof) <= 0) continue;
3518:     PetscSectionGetOffset(globalCoarse, p, &gOff);

3520:     rowOffsets[0]  = 0;
3521:     offsetsCopy[0] = 0;
3522:     if (numFields) {
3523:       PetscInt f;

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

3536:     PetscSectionGetDof(multiRootSec, p, &numLeaves);
3537:     PetscSectionGetOffset(multiRootSec, p, &leafStart);
3538:     leafEnd = leafStart + numLeaves;
3539:     for (l = leafStart; l < leafEnd; l++) {
3540:       PetscInt        numIndices, childId, offset;
3541:       const PetscInt *childIndices;

3543:       PetscSectionGetDof(rootIndicesSec, l, &numIndices);
3544:       PetscSectionGetOffset(rootIndicesSec, l, &offset);
3545:       childId      = rootIndices[offset++];
3546:       childIndices = &rootIndices[offset];
3547:       numIndices--;

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

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

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

3568:         lim        = PetscMax(1, numFields);
3569:         offsets[0] = 0;
3570:         if (numFields) {
3571:           PetscInt f;

3573:           for (f = 0; f < numFields; f++) {
3574:             PetscInt fDof;
3575:             PetscSectionGetFieldDof(cSecRef, childId, f, &fDof);

3577:             offsets[f + 1] = fDof + offsets[f];
3578:           }
3579:         } else {
3580:           PetscInt cDof;

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

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

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

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

3619:     PetscSectionGetDof(globalCoarse, p, &dof);
3620:     PetscSectionGetConstraintDof(globalCoarse, p, &cdof);
3621:     if ((dof - cdof) <= 0) continue;
3622:     PetscSectionGetOffset(globalCoarse, p, &gOff);

3624:     rowOffsets[0]  = 0;
3625:     offsetsCopy[0] = 0;
3626:     if (numFields) {
3627:       PetscInt f;

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

3640:     PetscSectionGetDof(multiRootSec, p, &numLeaves);
3641:     PetscSectionGetOffset(multiRootSec, p, &leafStart);
3642:     leafEnd = leafStart + numLeaves;
3643:     for (l = leafStart; l < leafEnd; l++) {
3644:       PetscInt        numIndices, childId, offset;
3645:       const PetscInt *childIndices;

3647:       PetscSectionGetDof(rootIndicesSec, l, &numIndices);
3648:       PetscSectionGetOffset(rootIndicesSec, l, &offset);
3649:       childId      = rootIndices[offset++];
3650:       childIndices = &rootIndices[offset];
3651:       numIndices--;

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

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

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

3662:         lim        = PetscMax(1, numFields);
3663:         offsets[0] = 0;
3664:         if (numFields) {
3665:           PetscInt f;

3667:           for (f = 0; f < numFields; f++) {
3668:             PetscInt fDof;
3669:             PetscSectionGetFieldDof(cSecRef, childId, f, &fDof);

3671:             offsets[f + 1] = fDof + offsets[f];
3672:           }
3673:         } else {
3674:           PetscInt cDof;

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

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

3696:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
3697:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
3698:   return 0;
3699: }

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

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

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

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

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

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

3761:   DMGetLocalSection(coarse, &localCoarse);
3762:   DMGetGlobalSection(coarse, &globalCoarse);

3764:   DMPlexGetAnchors(coarse, &aSec, &aIS);
3765:   ISGetIndices(aIS, &anchors);
3766:   PetscSectionGetChart(aSec, &aStart, &aEnd);

3768:   DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL);
3769:   PetscSectionGetChart(cSec, &cStart, &cEnd);

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

3782:     VecGetDM(cellGeom, &cellDM);
3783:     VecGetArrayRead(cellGeom, &cellGeomArray);
3784:     VecGetDM(grad, &gradDM);
3785:     VecGetArrayRead(grad, &gradArray);
3786:     for (i = 0; i < PetscMax(1, numFields); i++) {
3787:       PetscObject  obj;
3788:       PetscClassId id;

3790:       DMGetField(coarse, i, NULL, &obj);
3791:       PetscObjectGetClassId(obj, &id);
3792:       if (id == PETSCFV_CLASSID) {
3793:         fv = (PetscFV)obj;
3794:         PetscFVGetNumComponents(fv, &numFVcomps);
3795:         fvField = i;
3796:         break;
3797:       }
3798:     }
3799:   }

3801:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
3802:     PetscInt dof;
3803:     PetscInt maxChildId = maxChildIds[p - pStartC];
3804:     PetscInt numValues  = 0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3986:             DMPlexComputeCellGeometryFVM(fine, p, NULL, centroid, NULL);
3987:             for (i = 0; i < dim; i++) diff[i] = centroid[i] - parentCentroid[i];
3988:             for (i = 0; i < numFVcomps; i++) {
3989:               PetscScalar val = 0.;

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

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


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

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

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

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

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

4067:     PetscSectionGetDof(globalCoarse, p, &dof);
4068:     PetscSectionGetConstraintDof(globalCoarse, p, &cdof);
4069:     if ((dof - cdof) <= 0) continue;
4070:     PetscSectionGetDof(localCoarse, p, &dof);
4071:     PetscSectionGetOffset(globalCoarse, p, &gOff);

4073:     rowOffsets[0]  = 0;
4074:     offsetsCopy[0] = 0;
4075:     if (numFields) {
4076:       PetscInt f;

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

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

4097:       PetscSectionGetDof(rootIndicesSec, l, &numIndices);
4098:       PetscSectionGetOffset(rootIndicesSec, l, &offset);
4099:       childId     = (PetscInt)PetscRealPart(rootValues[offset++]);
4100:       childValues = &rootValues[offset];
4101:       numIndices--;

4103:       if (childId == -2) { /* skip */
4104:         continue;
4105:       } else if (childId == -1) { /* equivalent points: scatter */
4106:         PetscInt m;

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

4113:         contribute = PETSC_TRUE;
4114:         DMPlexGetTreeParent(refTree, childId, &parentId, NULL);

4116:         lim        = PetscMax(1, numFields);
4117:         offsets[0] = 0;
4118:         if (numFields) {
4119:           PetscInt f;

4121:           for (f = 0; f < numFields; f++) {
4122:             PetscInt fDof;
4123:             PetscSectionGetFieldDof(cSecRef, childId, f, &fDof);

4125:             offsets[f + 1] = fDof + offsets[f];
4126:           }
4127:         } else {
4128:           PetscInt cDof;

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

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

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

4164:   Collective on dmIn

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

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

4190:   Level: developer

4192: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscSF`, `Vec`, `PetscFV`, `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()`
4193: @*/
4194: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4195: {
4196:   VecSet(vecOut, 0.0);
4197:   if (sfRefine) {
4198:     Vec vecInLocal;
4199:     DM  dmGrad   = NULL;
4200:     Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;

4202:     DMGetLocalVector(dmIn, &vecInLocal);
4203:     VecSet(vecInLocal, 0.0);
4204:     {
4205:       PetscInt numFields, i;

4207:       DMGetNumFields(dmIn, &numFields);
4208:       for (i = 0; i < numFields; i++) {
4209:         PetscObject  obj;
4210:         PetscClassId classid;

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