Actual source code: plextree.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petsc/private/isimpl.h>
  3:  #include <petsc/private/petscfeimpl.h>
  4:  #include <petscsf.h>
  5:  #include <petscds.h>

  7: /** hierarchy routines */

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

 12:   Not collective

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

 18:   Level: intermediate

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

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

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

 39:   Not collective

 41:   Input Parameters:
 42: . dm - The DMPlex object

 44:   Output Parameters
 45: . ref - The reference tree DMPlex object

 47:   Level: intermediate

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

 58:   *ref = mesh->referenceTree;
 59:   return(0);
 60: }

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

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

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

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

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

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

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

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

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

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

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

193:   Level: developer

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

204:   if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented");
205:   mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);
206:   return(0);
207: }

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

211: PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
212: {

216:   DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);
217:   return(0);
218: }

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

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

258:     DMPlexGetHeightStratum(K, d, &cStart, &cEnd);
259:     for (c = cStart; c < cEnd; c++) {
260:       permvals[offset++] = c;
261:     }

263:     DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);
264:     for (c = cStart; c < cEnd; c++) {
265:       permvals[offset++] = c + (pEnd - pStart);
266:     }
267:   }
268:   ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);
269:   PetscSectionSetPermutation(unionSection,perm);
270:   PetscSectionSetUp(unionSection);
271:   PetscSectionGetStorageSize(unionSection,&numUnionPoints);
272:   PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);
273:   /* count dimension points */
274:   for (d = 0; d <= dim; d++) {
275:     PetscInt cStart, cOff, cOff2;
276:     DMPlexGetHeightStratum(K,d,&cStart,NULL);
277:     PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);
278:     if (d < dim) {
279:       DMPlexGetHeightStratum(K,d+1,&cStart,NULL);
280:       PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);
281:     }
282:     else {
283:       cOff2 = numUnionPoints;
284:     }
285:     numDimPoints[dim - d] = cOff2 - cOff;
286:   }
287:   PetscSectionCreate(comm, &unionConeSection);
288:   PetscSectionSetChart(unionConeSection, 0, numUnionPoints);
289:   /* count the cones in the union */
290:   for (p = pStart; p < pEnd; p++) {
291:     PetscInt dof, uOff;

293:     DMPlexGetConeSize(K, p, &dof);
294:     PetscSectionGetOffset(unionSection, p - pStart,&uOff);
295:     PetscSectionSetDof(unionConeSection, uOff, dof);
296:     coneSizes[uOff] = dof;
297:   }
298:   for (p = pRefStart; p < pRefEnd; p++) {
299:     PetscInt dof, uDof, uOff;

301:     DMPlexGetConeSize(Kref, p, &dof);
302:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
303:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
304:     if (uDof) {
305:       PetscSectionSetDof(unionConeSection, uOff, dof);
306:       coneSizes[uOff] = dof;
307:     }
308:   }
309:   PetscSectionSetUp(unionConeSection);
310:   PetscSectionGetStorageSize(unionConeSection,&numCones);
311:   PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);
312:   /* write the cones in the union */
313:   for (p = pStart; p < pEnd; p++) {
314:     PetscInt dof, uOff, c, cOff;
315:     const PetscInt *cone, *orientation;

317:     DMPlexGetConeSize(K, p, &dof);
318:     DMPlexGetCone(K, p, &cone);
319:     DMPlexGetConeOrientation(K, p, &orientation);
320:     PetscSectionGetOffset(unionSection, p - pStart,&uOff);
321:     PetscSectionGetOffset(unionConeSection,uOff,&cOff);
322:     for (c = 0; c < dof; c++) {
323:       PetscInt e, eOff;
324:       e                           = cone[c];
325:       PetscSectionGetOffset(unionSection, e - pStart, &eOff);
326:       unionCones[cOff + c]        = eOff;
327:       unionOrientations[cOff + c] = orientation[c];
328:     }
329:   }
330:   for (p = pRefStart; p < pRefEnd; p++) {
331:     PetscInt dof, uDof, uOff, c, cOff;
332:     const PetscInt *cone, *orientation;

334:     DMPlexGetConeSize(Kref, p, &dof);
335:     DMPlexGetCone(Kref, p, &cone);
336:     DMPlexGetConeOrientation(Kref, p, &orientation);
337:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
338:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
339:     if (uDof) {
340:       PetscSectionGetOffset(unionConeSection,uOff,&cOff);
341:       for (c = 0; c < dof; c++) {
342:         PetscInt e, eOff, eDof;

344:         e    = cone[c];
345:         PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);
346:         if (eDof) {
347:           PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);
348:         }
349:         else {
350:           DMLabelGetValue(identityRef, e, &e);
351:           PetscSectionGetOffset(unionSection, e - pStart, &eOff);
352:         }
353:         unionCones[cOff + c]        = eOff;
354:         unionOrientations[cOff + c] = orientation[c];
355:       }
356:     }
357:   }
358:   /* get the coordinates */
359:   {
360:     PetscInt     vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff;
361:     PetscSection KcoordsSec, KrefCoordsSec;
362:     Vec          KcoordsVec, KrefCoordsVec;
363:     PetscScalar *Kcoords;

365:     DMGetCoordinateSection(K, &KcoordsSec);
366:     DMGetCoordinatesLocal(K, &KcoordsVec);
367:     DMGetCoordinateSection(Kref, &KrefCoordsSec);
368:     DMGetCoordinatesLocal(Kref, &KrefCoordsVec);

370:     numVerts = numDimPoints[0];
371:     PetscMalloc1(numVerts * dim,&unionCoords);
372:     DMPlexGetDepthStratum(K,0,&vStart,&vEnd);

374:     offset = 0;
375:     for (v = vStart; v < vEnd; v++) {
376:       PetscSectionGetOffset(unionSection,v - pStart,&vOff);
377:       VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);
378:       for (d = 0; d < dim; d++) {
379:         unionCoords[offset * dim + d] = Kcoords[d];
380:       }
381:       offset++;
382:     }
383:     DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);
384:     for (v = vRefStart; v < vRefEnd; v++) {
385:       PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);
386:       PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);
387:       VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);
388:       if (vDof) {
389:         for (d = 0; d < dim; d++) {
390:           unionCoords[offset * dim + d] = Kcoords[d];
391:         }
392:         offset++;
393:       }
394:     }
395:   }
396:   DMCreate(comm,ref);
397:   DMSetType(*ref,DMPLEX);
398:   DMSetDimension(*ref,dim);
399:   DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);
400:   /* set the tree */
401:   PetscSectionCreate(comm,&parentSection);
402:   PetscSectionSetChart(parentSection,0,numUnionPoints);
403:   for (p = pRefStart; p < pRefEnd; p++) {
404:     PetscInt uDof, uOff;

406:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
407:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
408:     if (uDof) {
409:       PetscSectionSetDof(parentSection,uOff,1);
410:     }
411:   }
412:   PetscSectionSetUp(parentSection);
413:   PetscSectionGetStorageSize(parentSection,&parentSize);
414:   PetscMalloc2(parentSize,&parents,parentSize,&childIDs);
415:   for (p = pRefStart; p < pRefEnd; p++) {
416:     PetscInt uDof, uOff;

418:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
419:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
420:     if (uDof) {
421:       PetscInt pOff, parent, parentU;
422:       PetscSectionGetOffset(parentSection,uOff,&pOff);
423:       DMLabelGetValue(identityRef,p,&parent);
424:       PetscSectionGetOffset(unionSection, parent - pStart,&parentU);
425:       parents[pOff] = parentU;
426:       childIDs[pOff] = uOff;
427:     }
428:   }
429:   DMPlexCreateReferenceTree_SetTree(*ref,parentSection,parents,childIDs);
430:   PetscSectionDestroy(&parentSection);
431:   PetscFree2(parents,childIDs);

433:   /* clean up */
434:   PetscSectionDestroy(&unionSection);
435:   PetscSectionDestroy(&unionConeSection);
436:   ISDestroy(&perm);
437:   PetscFree(unionCoords);
438:   PetscFree2(unionCones,unionOrientations);
439:   PetscFree2(coneSizes,numDimPoints);
440:   return(0);
441: }

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

446:   Collective on comm

448:   Input Parameters:
449: + comm    - the MPI communicator
450: . dim     - the spatial dimension
451: - simplex - Flag for simplex, otherwise use a tensor-product cell

453:   Output Parameters:
454: . ref     - the reference tree DMPlex object

456:   Level: intermediate

458: .keywords: reference cell
459: .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree()
460: @*/
461: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
462: {
463:   DM_Plex       *mesh;
464:   DM             K, Kref;
465:   PetscInt       p, pStart, pEnd;
466:   DMLabel        identity;

470: #if 1
471:   comm = PETSC_COMM_SELF;
472: #endif
473:   /* create a reference element */
474:   DMPlexCreateReferenceCell(comm, dim, simplex, &K);
475:   DMCreateLabel(K, "identity");
476:   DMGetLabel(K, "identity", &identity);
477:   DMPlexGetChart(K, &pStart, &pEnd);
478:   for (p = pStart; p < pEnd; p++) {
479:     DMLabelSetValue(identity, p, p);
480:   }
481:   /* refine it */
482:   DMRefine(K,comm,&Kref);

484:   /* the reference tree is the union of these two, without duplicating
485:    * points that appear in both */
486:   DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref);
487:   mesh = (DM_Plex *) (*ref)->data;
488:   mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
489:   DMDestroy(&K);
490:   DMDestroy(&Kref);
491:   return(0);
492: }

494: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
495: {
496:   DM_Plex        *mesh = (DM_Plex *)dm->data;
497:   PetscSection   childSec, pSec;
498:   PetscInt       p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT;
499:   PetscInt       *offsets, *children, pStart, pEnd;

504:   PetscSectionDestroy(&mesh->childSection);
505:   PetscFree(mesh->children);
506:   pSec = mesh->parentSection;
507:   if (!pSec) return(0);
508:   PetscSectionGetStorageSize(pSec,&pSize);
509:   for (p = 0; p < pSize; p++) {
510:     PetscInt par = mesh->parents[p];

512:     parMax = PetscMax(parMax,par+1);
513:     parMin = PetscMin(parMin,par);
514:   }
515:   if (parMin > parMax) {
516:     parMin = -1;
517:     parMax = -1;
518:   }
519:   PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);
520:   PetscSectionSetChart(childSec,parMin,parMax);
521:   for (p = 0; p < pSize; p++) {
522:     PetscInt par = mesh->parents[p];

524:     PetscSectionAddDof(childSec,par,1);
525:   }
526:   PetscSectionSetUp(childSec);
527:   PetscSectionGetStorageSize(childSec,&cSize);
528:   PetscMalloc1(cSize,&children);
529:   PetscCalloc1(parMax-parMin,&offsets);
530:   PetscSectionGetChart(pSec,&pStart,&pEnd);
531:   for (p = pStart; p < pEnd; p++) {
532:     PetscInt dof, off, i;

534:     PetscSectionGetDof(pSec,p,&dof);
535:     PetscSectionGetOffset(pSec,p,&off);
536:     for (i = 0; i < dof; i++) {
537:       PetscInt par = mesh->parents[off + i], cOff;

539:       PetscSectionGetOffset(childSec,par,&cOff);
540:       children[cOff + offsets[par-parMin]++] = p;
541:     }
542:   }
543:   mesh->childSection = childSec;
544:   mesh->children = children;
545:   PetscFree(offsets);
546:   return(0);
547: }

549: static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
550: {
551:   PetscInt       pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
552:   const PetscInt *vals;
553:   PetscSection   secNew;
554:   PetscBool      anyNew, globalAnyNew;
555:   PetscBool      compress;

559:   PetscSectionGetChart(section,&pStart,&pEnd);
560:   ISGetLocalSize(is,&size);
561:   ISGetIndices(is,&vals);
562:   PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);
563:   PetscSectionSetChart(secNew,pStart,pEnd);
564:   for (i = 0; i < size; i++) {
565:     PetscInt dof;

567:     p = vals[i];
568:     if (p < pStart || p >= pEnd) continue;
569:     PetscSectionGetDof(section, p, &dof);
570:     if (dof) break;
571:   }
572:   if (i == size) {
573:     PetscSectionSetUp(secNew);
574:     anyNew   = PETSC_FALSE;
575:     compress = PETSC_FALSE;
576:     sizeNew  = 0;
577:   }
578:   else {
579:     anyNew = PETSC_TRUE;
580:     for (p = pStart; p < pEnd; p++) {
581:       PetscInt dof, off;

583:       PetscSectionGetDof(section, p, &dof);
584:       PetscSectionGetOffset(section, p, &off);
585:       for (i = 0; i < dof; i++) {
586:         PetscInt q = vals[off + i], qDof = 0;

588:         if (q >= pStart && q < pEnd) {
589:           PetscSectionGetDof(section, q, &qDof);
590:         }
591:         if (qDof) {
592:           PetscSectionAddDof(secNew, p, qDof);
593:         }
594:         else {
595:           PetscSectionAddDof(secNew, p, 1);
596:         }
597:       }
598:     }
599:     PetscSectionSetUp(secNew);
600:     PetscSectionGetStorageSize(secNew,&sizeNew);
601:     PetscMalloc1(sizeNew,&valsNew);
602:     compress = PETSC_FALSE;
603:     for (p = pStart; p < pEnd; p++) {
604:       PetscInt dof, off, count, offNew, dofNew;

606:       PetscSectionGetDof(section, p, &dof);
607:       PetscSectionGetOffset(section, p, &off);
608:       PetscSectionGetDof(secNew, p, &dofNew);
609:       PetscSectionGetOffset(secNew, p, &offNew);
610:       count = 0;
611:       for (i = 0; i < dof; i++) {
612:         PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;

614:         if (q >= pStart && q < pEnd) {
615:           PetscSectionGetDof(section, q, &qDof);
616:           PetscSectionGetOffset(section, q, &qOff);
617:         }
618:         if (qDof) {
619:           PetscInt oldCount = count;

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

624:             for (k = 0; k < oldCount; k++) {
625:               if (valsNew[offNew + k] == r) {
626:                 break;
627:               }
628:             }
629:             if (k == oldCount) {
630:               valsNew[offNew + count++] = r;
631:             }
632:           }
633:         }
634:         else {
635:           PetscInt k, oldCount = count;

637:           for (k = 0; k < oldCount; k++) {
638:             if (valsNew[offNew + k] == q) {
639:               break;
640:             }
641:           }
642:           if (k == oldCount) {
643:             valsNew[offNew + count++] = q;
644:           }
645:         }
646:       }
647:       if (count < dofNew) {
648:         PetscSectionSetDof(secNew, p, count);
649:         compress = PETSC_TRUE;
650:       }
651:     }
652:   }
653:   ISRestoreIndices(is,&vals);
654:   MPIU_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
655:   if (!globalAnyNew) {
656:     PetscSectionDestroy(&secNew);
657:     *sectionNew = NULL;
658:     *isNew = NULL;
659:   }
660:   else {
661:     PetscBool globalCompress;

663:     MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
664:     if (compress) {
665:       PetscSection secComp;
666:       PetscInt *valsComp = NULL;

668:       PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);
669:       PetscSectionSetChart(secComp,pStart,pEnd);
670:       for (p = pStart; p < pEnd; p++) {
671:         PetscInt dof;

673:         PetscSectionGetDof(secNew, p, &dof);
674:         PetscSectionSetDof(secComp, p, dof);
675:       }
676:       PetscSectionSetUp(secComp);
677:       PetscSectionGetStorageSize(secComp,&sizeNew);
678:       PetscMalloc1(sizeNew,&valsComp);
679:       for (p = pStart; p < pEnd; p++) {
680:         PetscInt dof, off, offNew, j;

682:         PetscSectionGetDof(secNew, p, &dof);
683:         PetscSectionGetOffset(secNew, p, &off);
684:         PetscSectionGetOffset(secComp, p, &offNew);
685:         for (j = 0; j < dof; j++) {
686:           valsComp[offNew + j] = valsNew[off + j];
687:         }
688:       }
689:       PetscSectionDestroy(&secNew);
690:       secNew  = secComp;
691:       PetscFree(valsNew);
692:       valsNew = valsComp;
693:     }
694:     ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);
695:   }
696:   return(0);
697: }

699: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
700: {
701:   PetscInt       p, pStart, pEnd, *anchors, size;
702:   PetscInt       aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
703:   PetscSection   aSec;
704:   DMLabel        canonLabel;
705:   IS             aIS;

710:   DMPlexGetChart(dm,&pStart,&pEnd);
711:   DMGetLabel(dm,"canonical",&canonLabel);
712:   for (p = pStart; p < pEnd; p++) {
713:     PetscInt parent;

715:     if (canonLabel) {
716:       PetscInt canon;

718:       DMLabelGetValue(canonLabel,p,&canon);
719:       if (p != canon) continue;
720:     }
721:     DMPlexGetTreeParent(dm,p,&parent,NULL);
722:     if (parent != p) {
723:       aMin = PetscMin(aMin,p);
724:       aMax = PetscMax(aMax,p+1);
725:     }
726:   }
727:   if (aMin > aMax) {
728:     aMin = -1;
729:     aMax = -1;
730:   }
731:   PetscSectionCreate(PETSC_COMM_SELF,&aSec);
732:   PetscSectionSetChart(aSec,aMin,aMax);
733:   for (p = aMin; p < aMax; p++) {
734:     PetscInt parent, ancestor = p;

736:     if (canonLabel) {
737:       PetscInt canon;

739:       DMLabelGetValue(canonLabel,p,&canon);
740:       if (p != canon) continue;
741:     }
742:     DMPlexGetTreeParent(dm,p,&parent,NULL);
743:     while (parent != ancestor) {
744:       ancestor = parent;
745:       DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
746:     }
747:     if (ancestor != p) {
748:       PetscInt closureSize, *closure = NULL;

750:       DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
751:       PetscSectionSetDof(aSec,p,closureSize);
752:       DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
753:     }
754:   }
755:   PetscSectionSetUp(aSec);
756:   PetscSectionGetStorageSize(aSec,&size);
757:   PetscMalloc1(size,&anchors);
758:   for (p = aMin; p < aMax; p++) {
759:     PetscInt parent, ancestor = p;

761:     if (canonLabel) {
762:       PetscInt canon;

764:       DMLabelGetValue(canonLabel,p,&canon);
765:       if (p != canon) continue;
766:     }
767:     DMPlexGetTreeParent(dm,p,&parent,NULL);
768:     while (parent != ancestor) {
769:       ancestor = parent;
770:       DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
771:     }
772:     if (ancestor != p) {
773:       PetscInt j, closureSize, *closure = NULL, aOff;

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

777:       DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
778:       for (j = 0; j < closureSize; j++) {
779:         anchors[aOff + j] = closure[2*j];
780:       }
781:       DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
782:     }
783:   }
784:   ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);
785:   {
786:     PetscSection aSecNew = aSec;
787:     IS           aISNew  = aIS;

789:     PetscObjectReference((PetscObject)aSec);
790:     PetscObjectReference((PetscObject)aIS);
791:     while (aSecNew) {
792:       PetscSectionDestroy(&aSec);
793:       ISDestroy(&aIS);
794:       aSec    = aSecNew;
795:       aIS     = aISNew;
796:       aSecNew = NULL;
797:       aISNew  = NULL;
798:       AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);
799:     }
800:   }
801:   DMPlexSetAnchors(dm,aSec,aIS);
802:   PetscSectionDestroy(&aSec);
803:   ISDestroy(&aIS);
804:   return(0);
805: }

807: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp)
808: {

812:   if (numTrueSupp[p] == -1) {
813:     PetscInt i, alldof;
814:     const PetscInt *supp;
815:     PetscInt count = 0;

817:     DMPlexGetSupportSize(dm,p,&alldof);
818:     DMPlexGetSupport(dm,p,&supp);
819:     for (i = 0; i < alldof; i++) {
820:       PetscInt q = supp[i], numCones, j;
821:       const PetscInt *cone;

823:       DMPlexGetConeSize(dm,q,&numCones);
824:       DMPlexGetCone(dm,q,&cone);
825:       for (j = 0; j < numCones; j++) {
826:         if (cone[j] == p) break;
827:       }
828:       if (j < numCones) count++;
829:     }
830:     numTrueSupp[p] = count;
831:   }
832:   *dof = numTrueSupp[p];
833:   return(0);
834: }

836: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
837: {
838:   DM_Plex *mesh = (DM_Plex *)dm->data;
839:   PetscSection newSupportSection;
840:   PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth;
841:   PetscInt *numTrueSupp;
842:   PetscInt *offsets;

847:   /* symmetrize the hierarchy */
848:   DMPlexGetDepth(dm,&depth);
849:   PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);
850:   DMPlexGetChart(dm,&pStart,&pEnd);
851:   PetscSectionSetChart(newSupportSection,pStart,pEnd);
852:   PetscCalloc1(pEnd,&offsets);
853:   PetscMalloc1(pEnd,&numTrueSupp);
854:   for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
855:   /* if a point is in the (true) support of q, it should be in the support of
856:    * parent(q) */
857:   for (d = 0; d <= depth; d++) {
858:     DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
859:     for (p = pStart; p < pEnd; ++p) {
860:       PetscInt dof, q, qdof, parent;

862:       DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);
863:       PetscSectionAddDof(newSupportSection, p, dof);
864:       q    = p;
865:       DMPlexGetTreeParent(dm,q,&parent,NULL);
866:       while (parent != q && parent >= pStart && parent < pEnd) {
867:         q = parent;

869:         DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);
870:         PetscSectionAddDof(newSupportSection,p,qdof);
871:         PetscSectionAddDof(newSupportSection,q,dof);
872:         DMPlexGetTreeParent(dm,q,&parent,NULL);
873:       }
874:     }
875:   }
876:   PetscSectionSetUp(newSupportSection);
877:   PetscSectionGetStorageSize(newSupportSection,&newSize);
878:   PetscMalloc1(newSize,&newSupports);
879:   for (d = 0; d <= depth; d++) {
880:     DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
881:     for (p = pStart; p < pEnd; p++) {
882:       PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;

884:       PetscSectionGetDof(mesh->supportSection, p, &dof);
885:       PetscSectionGetOffset(mesh->supportSection, p, &off);
886:       PetscSectionGetDof(newSupportSection, p, &newDof);
887:       PetscSectionGetOffset(newSupportSection, p, &newOff);
888:       for (i = 0; i < dof; i++) {
889:         PetscInt numCones, j;
890:         const PetscInt *cone;
891:         PetscInt q = mesh->supports[off + i];

893:         DMPlexGetConeSize(dm,q,&numCones);
894:         DMPlexGetCone(dm,q,&cone);
895:         for (j = 0; j < numCones; j++) {
896:           if (cone[j] == p) break;
897:         }
898:         if (j < numCones) newSupports[newOff+offsets[p]++] = q;
899:       }
900:       mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof);

902:       q    = p;
903:       DMPlexGetTreeParent(dm,q,&parent,NULL);
904:       while (parent != q && parent >= pStart && parent < pEnd) {
905:         q = parent;
906:         PetscSectionGetDof(mesh->supportSection, q, &qdof);
907:         PetscSectionGetOffset(mesh->supportSection, q, &qoff);
908:         PetscSectionGetOffset(newSupportSection, q, &newqOff);
909:         for (i = 0; i < qdof; i++) {
910:           PetscInt numCones, j;
911:           const PetscInt *cone;
912:           PetscInt r = mesh->supports[qoff + i];

914:           DMPlexGetConeSize(dm,r,&numCones);
915:           DMPlexGetCone(dm,r,&cone);
916:           for (j = 0; j < numCones; j++) {
917:             if (cone[j] == q) break;
918:           }
919:           if (j < numCones) newSupports[newOff+offsets[p]++] = r;
920:         }
921:         for (i = 0; i < dof; i++) {
922:           PetscInt numCones, j;
923:           const PetscInt *cone;
924:           PetscInt r = mesh->supports[off + i];

926:           DMPlexGetConeSize(dm,r,&numCones);
927:           DMPlexGetCone(dm,r,&cone);
928:           for (j = 0; j < numCones; j++) {
929:             if (cone[j] == p) break;
930:           }
931:           if (j < numCones) newSupports[newqOff+offsets[q]++] = r;
932:         }
933:         DMPlexGetTreeParent(dm,q,&parent,NULL);
934:       }
935:     }
936:   }
937:   PetscSectionDestroy(&mesh->supportSection);
938:   mesh->supportSection = newSupportSection;
939:   PetscFree(mesh->supports);
940:   mesh->supports = newSupports;
941:   PetscFree(offsets);
942:   PetscFree(numTrueSupp);

944:   return(0);
945: }

947: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat);
948: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat);

950: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
951: {
952:   DM_Plex       *mesh = (DM_Plex *)dm->data;
953:   DM             refTree;
954:   PetscInt       size;

960:   PetscObjectReference((PetscObject)parentSection);
961:   PetscSectionDestroy(&mesh->parentSection);
962:   mesh->parentSection = parentSection;
963:   PetscSectionGetStorageSize(parentSection,&size);
964:   if (parents != mesh->parents) {
965:     PetscFree(mesh->parents);
966:     PetscMalloc1(size,&mesh->parents);
967:     PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));
968:   }
969:   if (childIDs != mesh->childIDs) {
970:     PetscFree(mesh->childIDs);
971:     PetscMalloc1(size,&mesh->childIDs);
972:     PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));
973:   }
974:   DMPlexGetReferenceTree(dm,&refTree);
975:   if (refTree) {
976:     DMLabel canonLabel;

978:     DMGetLabel(refTree,"canonical",&canonLabel);
979:     if (canonLabel) {
980:       PetscInt i;

982:       for (i = 0; i < size; i++) {
983:         PetscInt canon;
984:         DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);
985:         if (canon >= 0) {
986:           mesh->childIDs[i] = canon;
987:         }
988:       }
989:     }
990:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
991:   }
992:   else {
993:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
994:   }
995:   DMPlexTreeSymmetrize(dm);
996:   if (computeCanonical) {
997:     PetscInt d, dim;

999:     /* add the canonical label */
1000:     DMGetDimension(dm,&dim);
1001:     DMCreateLabel(dm,"canonical");
1002:     for (d = 0; d <= dim; d++) {
1003:       PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
1004:       const PetscInt *cChildren;

1006:       DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
1007:       for (p = dStart; p < dEnd; p++) {
1008:         DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);
1009:         if (cNumChildren) {
1010:           canon = p;
1011:           break;
1012:         }
1013:       }
1014:       if (canon == -1) continue;
1015:       for (p = dStart; p < dEnd; p++) {
1016:         PetscInt numChildren, i;
1017:         const PetscInt *children;

1019:         DMPlexGetTreeChildren(dm,p,&numChildren,&children);
1020:         if (numChildren) {
1021:           if (numChildren != cNumChildren) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"All parent points in a stratum should have the same number of children: %d != %d", numChildren, cNumChildren);
1022:           DMSetLabelValue(dm,"canonical",p,canon);
1023:           for (i = 0; i < numChildren; i++) {
1024:             DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);
1025:           }
1026:         }
1027:       }
1028:     }
1029:   }
1030:   if (exchangeSupports) {
1031:     DMPlexTreeExchangeSupports(dm);
1032:   }
1033:   mesh->createanchors = DMPlexCreateAnchors_Tree;
1034:   /* reset anchors */
1035:   DMPlexSetAnchors(dm,NULL,NULL);
1036:   return(0);
1037: }

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

1044:   Collective on dm

1046:   Input Parameters:
1047: + dm - the DMPlex object
1048: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1049:                   offset indexes the parent and childID list; the reference count of parentSection is incremented
1050: . parents - a list of the point parents; copied, can be destroyed
1051: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1052:              the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed

1054:   Level: intermediate

1056: .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1057: @*/
1058: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
1059: {

1063:   DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);
1064:   return(0);
1065: }

1067: /*@
1068:   DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1069:   Collective on dm

1071:   Input Parameters:
1072: . dm - the DMPlex object

1074:   Output Parameters:
1075: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1076:                   offset indexes the parent and childID list
1077: . parents - a list of the point parents
1078: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1079:              the child corresponds to the point in the reference tree with index childID
1080: . childSection - the inverse of the parent section
1081: - children - a list of the point children

1083:   Level: intermediate

1085: .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1086: @*/
1087: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1088: {
1089:   DM_Plex        *mesh = (DM_Plex *)dm->data;

1093:   if (parentSection) *parentSection = mesh->parentSection;
1094:   if (parents)       *parents       = mesh->parents;
1095:   if (childIDs)      *childIDs      = mesh->childIDs;
1096:   if (childSection)  *childSection  = mesh->childSection;
1097:   if (children)      *children      = mesh->children;
1098:   return(0);
1099: }

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

1104:   Input Parameters:
1105: + dm - the DMPlex object
1106: - point - the query point

1108:   Output Parameters:
1109: + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent
1110: - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point
1111:             does not have a parent

1113:   Level: intermediate

1115: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren()
1116: @*/
1117: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1118: {
1119:   DM_Plex       *mesh = (DM_Plex *)dm->data;
1120:   PetscSection   pSec;

1125:   pSec = mesh->parentSection;
1126:   if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1127:     PetscInt dof;

1129:     PetscSectionGetDof (pSec, point, &dof);
1130:     if (dof) {
1131:       PetscInt off;

1133:       PetscSectionGetOffset (pSec, point, &off);
1134:       if (parent)  *parent = mesh->parents[off];
1135:       if (childID) *childID = mesh->childIDs[off];
1136:       return(0);
1137:     }
1138:   }
1139:   if (parent) {
1140:     *parent = point;
1141:   }
1142:   if (childID) {
1143:     *childID = 0;
1144:   }
1145:   return(0);
1146: }

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

1151:   Input Parameters:
1152: + dm - the DMPlex object
1153: - point - the query point

1155:   Output Parameters:
1156: + numChildren - if not NULL, set to the number of children
1157: - children - if not NULL, set to a list children, or set to NULL if the point has no children

1159:   Level: intermediate

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

1165: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent()
1166: @*/
1167: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1168: {
1169:   DM_Plex       *mesh = (DM_Plex *)dm->data;
1170:   PetscSection   childSec;
1171:   PetscInt       dof = 0;

1176:   childSec = mesh->childSection;
1177:   if (childSec && point >= childSec->pStart && point < childSec->pEnd) {
1178:     PetscSectionGetDof (childSec, point, &dof);
1179:   }
1180:   if (numChildren) *numChildren = dof;
1181:   if (children) {
1182:     if (dof) {
1183:       PetscInt off;

1185:       PetscSectionGetOffset (childSec, point, &off);
1186:       *children = &mesh->children[off];
1187:     }
1188:     else {
1189:       *children = NULL;
1190:     }
1191:   }
1192:   return(0);
1193: }

1195: 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)
1196: {
1197:   PetscInt       f, b, p, c, offset, qPoints;

1201:   PetscSpaceEvaluate(space,nPoints,points,work,NULL,NULL);
1202:   for (f = 0, offset = 0; f < nFunctionals; f++) {
1203:     qPoints = pointsPerFn[f];
1204:     for (b = 0; b < nBasis; b++) {
1205:       PetscScalar val = 0.;

1207:       for (p = 0; p < qPoints; p++) {
1208:         for (c = 0; c < nComps; c++) {
1209:           val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c];
1210:         }
1211:       }
1212:       MatSetValue(basisAtPoints,b,f,val,INSERT_VALUES);
1213:     }
1214:     offset += qPoints;
1215:   }
1216:   MatAssemblyBegin(basisAtPoints,MAT_FINAL_ASSEMBLY);
1217:   MatAssemblyEnd(basisAtPoints,MAT_FINAL_ASSEMBLY);
1218:   return(0);
1219: }

1221: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
1222: {
1223:   PetscDS        ds;
1224:   PetscInt       spdim;
1225:   PetscInt       numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
1226:   const PetscInt *anchors;
1227:   PetscSection   aSec;
1228:   PetscReal      *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
1229:   IS             aIS;

1233:   DMPlexGetChart(dm,&pStart,&pEnd);
1234:   DMGetDS(dm,&ds);
1235:   PetscDSGetNumFields(ds,&numFields);
1236:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1237:   DMPlexGetAnchors(dm,&aSec,&aIS);
1238:   ISGetIndices(aIS,&anchors);
1239:   PetscSectionGetChart(cSec,&conStart,&conEnd);
1240:   DMGetDimension(dm,&spdim);
1241:   PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);

1243:   for (f = 0; f < numFields; f++) {
1244:     PetscObject disc;
1245:     PetscClassId id;
1246:     PetscSpace     bspace;
1247:     PetscDualSpace dspace;
1248:     PetscInt i, j, k, nPoints, Nc, offset;
1249:     PetscInt fSize, maxDof;
1250:     PetscReal   *weights, *pointsRef, *pointsReal, *work;
1251:     PetscScalar *scwork, *X;
1252:     PetscInt  *sizes, *workIndRow, *workIndCol;
1253:     Mat Amat, Bmat, Xmat;
1254:     const PetscInt    *numDof  = NULL;
1255:     const PetscInt    ***perms = NULL;
1256:     const PetscScalar ***flips = NULL;

1258:     PetscDSGetDiscretization(ds,f,&disc);
1259:     PetscObjectGetClassId(disc,&id);
1260:     if (id == PETSCFE_CLASSID) {
1261:       PetscFE fe = (PetscFE) disc;

1263:       PetscFEGetBasisSpace(fe,&bspace);
1264:       PetscFEGetDualSpace(fe,&dspace);
1265:       PetscDualSpaceGetDimension(dspace,&fSize);
1266:       PetscFEGetNumComponents(fe,&Nc);
1267:     }
1268:     else if (id == PETSCFV_CLASSID) {
1269:       PetscFV fv = (PetscFV) disc;

1271:       PetscFVGetNumComponents(fv,&Nc);
1272:       PetscSpaceCreate(PetscObjectComm((PetscObject)fv),&bspace);
1273:       PetscSpaceSetType(bspace,PETSCSPACEPOLYNOMIAL);
1274:       PetscSpaceSetDegree(bspace,0,PETSC_DETERMINE);
1275:       PetscSpaceSetNumComponents(bspace,Nc);
1276:       PetscSpaceSetNumVariables(bspace,spdim);
1277:       PetscSpaceSetUp(bspace);
1278:       PetscFVGetDualSpace(fv,&dspace);
1279:       PetscDualSpaceGetDimension(dspace,&fSize);
1280:     }
1281:     else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1282:     PetscDualSpaceGetNumDof(dspace,&numDof);
1283:     for (i = 0, maxDof = 0; i <= spdim; i++) {maxDof = PetscMax(maxDof,numDof[i]);}
1284:     PetscDualSpaceGetSymmetries(dspace,&perms,&flips);

1286:     MatCreate(PETSC_COMM_SELF,&Amat);
1287:     MatSetSizes(Amat,fSize,fSize,fSize,fSize);
1288:     MatSetType(Amat,MATSEQDENSE);
1289:     MatSetUp(Amat);
1290:     MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);
1291:     MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);
1292:     nPoints = 0;
1293:     for (i = 0; i < fSize; i++) {
1294:       PetscInt        qPoints, thisNc;
1295:       PetscQuadrature quad;

1297:       PetscDualSpaceGetFunctional(dspace,i,&quad);
1298:       PetscQuadratureGetData(quad,NULL,&thisNc,&qPoints,NULL,NULL);
1299:       if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc);
1300:       nPoints += qPoints;
1301:     }
1302:     PetscMalloc7(fSize,&sizes,nPoints*Nc,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal,nPoints*fSize*Nc,&work,maxDof,&workIndRow,maxDof,&workIndCol);
1303:     PetscMalloc1(maxDof * maxDof,&scwork);
1304:     offset = 0;
1305:     for (i = 0; i < fSize; i++) {
1306:       PetscInt        qPoints;
1307:       const PetscReal    *p, *w;
1308:       PetscQuadrature quad;

1310:       PetscDualSpaceGetFunctional(dspace,i,&quad);
1311:       PetscQuadratureGetData(quad,NULL,NULL,&qPoints,&p,&w);
1312:       PetscMemcpy(weights+Nc*offset,w,Nc*qPoints*sizeof(*w));
1313:       PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));
1314:       sizes[i] = qPoints;
1315:       offset  += qPoints;
1316:     }
1317:     EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsRef,weights,work,Amat);
1318:     MatLUFactor(Amat,NULL,NULL,NULL);
1319:     for (c = cStart; c < cEnd; c++) {
1320:       PetscInt        parent;
1321:       PetscInt        closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1322:       PetscInt        *childOffsets, *parentOffsets;

1324:       DMPlexGetTreeParent(dm,c,&parent,NULL);
1325:       if (parent == c) continue;
1326:       DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1327:       for (i = 0; i < closureSize; i++) {
1328:         PetscInt p = closure[2*i];
1329:         PetscInt conDof;

1331:         if (p < conStart || p >= conEnd) continue;
1332:         if (numFields) {
1333:           PetscSectionGetFieldDof(cSec,p,f,&conDof);
1334:         }
1335:         else {
1336:           PetscSectionGetDof(cSec,p,&conDof);
1337:         }
1338:         if (conDof) break;
1339:       }
1340:       if (i == closureSize) {
1341:         DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1342:         continue;
1343:       }

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

1350:         CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i*spdim],vtmp);
1351:         CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);
1352:       }
1353:       EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsReal,weights,work,Bmat);
1354:       MatMatSolve(Amat,Bmat,Xmat);
1355:       MatDenseGetArray(Xmat,&X);
1356:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1357:       PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);
1358:       childOffsets[0] = 0;
1359:       for (i = 0; i < closureSize; i++) {
1360:         PetscInt p = closure[2*i];
1361:         PetscInt dof;

1363:         if (numFields) {
1364:           PetscSectionGetFieldDof(section,p,f,&dof);
1365:         }
1366:         else {
1367:           PetscSectionGetDof(section,p,&dof);
1368:         }
1369:         childOffsets[i+1]=childOffsets[i]+dof;
1370:       }
1371:       parentOffsets[0] = 0;
1372:       for (i = 0; i < closureSizeP; i++) {
1373:         PetscInt p = closureP[2*i];
1374:         PetscInt dof;

1376:         if (numFields) {
1377:           PetscSectionGetFieldDof(section,p,f,&dof);
1378:         }
1379:         else {
1380:           PetscSectionGetDof(section,p,&dof);
1381:         }
1382:         parentOffsets[i+1]=parentOffsets[i]+dof;
1383:       }
1384:       for (i = 0; i < closureSize; i++) {
1385:         PetscInt conDof, conOff, aDof, aOff, nWork;
1386:         PetscInt p = closure[2*i];
1387:         PetscInt o = closure[2*i+1];
1388:         const PetscInt    *perm;
1389:         const PetscScalar *flip;

1391:         if (p < conStart || p >= conEnd) continue;
1392:         if (numFields) {
1393:           PetscSectionGetFieldDof(cSec,p,f,&conDof);
1394:           PetscSectionGetFieldOffset(cSec,p,f,&conOff);
1395:         }
1396:         else {
1397:           PetscSectionGetDof(cSec,p,&conDof);
1398:           PetscSectionGetOffset(cSec,p,&conOff);
1399:         }
1400:         if (!conDof) continue;
1401:         perm  = (perms && perms[i]) ? perms[i][o] : NULL;
1402:         flip  = (flips && flips[i]) ? flips[i][o] : NULL;
1403:         PetscSectionGetDof(aSec,p,&aDof);
1404:         PetscSectionGetOffset(aSec,p,&aOff);
1405:         nWork = childOffsets[i+1]-childOffsets[i];
1406:         for (k = 0; k < aDof; k++) {
1407:           PetscInt a = anchors[aOff + k];
1408:           PetscInt aSecDof, aSecOff;

1410:           if (numFields) {
1411:             PetscSectionGetFieldDof(section,a,f,&aSecDof);
1412:             PetscSectionGetFieldOffset(section,a,f,&aSecOff);
1413:           }
1414:           else {
1415:             PetscSectionGetDof(section,a,&aSecDof);
1416:             PetscSectionGetOffset(section,a,&aSecOff);
1417:           }
1418:           if (!aSecDof) continue;

1420:           for (j = 0; j < closureSizeP; j++) {
1421:             PetscInt q = closureP[2*j];
1422:             PetscInt oq = closureP[2*j+1];

1424:             if (q == a) {
1425:               PetscInt           r, s, nWorkP;
1426:               const PetscInt    *permP;
1427:               const PetscScalar *flipP;

1429:               permP  = (perms && perms[j]) ? perms[j][oq] : NULL;
1430:               flipP  = (flips && flips[j]) ? flips[j][oq] : NULL;
1431:               nWorkP = parentOffsets[j+1]-parentOffsets[j];
1432:               /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the
1433:                * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArray is
1434:                * column-major, so transpose-transpose = do nothing */
1435:               for (r = 0; r < nWork; r++) {
1436:                 for (s = 0; s < nWorkP; s++) {
1437:                   scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])];
1438:                 }
1439:               }
1440:               for (r = 0; r < nWork; r++)  {workIndRow[perm  ? perm[r]  : r] = conOff  + r;}
1441:               for (s = 0; s < nWorkP; s++) {workIndCol[permP ? permP[s] : s] = aSecOff + s;}
1442:               if (flip) {
1443:                 for (r = 0; r < nWork; r++) {
1444:                   for (s = 0; s < nWorkP; s++) {
1445:                     scwork[r * nWorkP + s] *= flip[r];
1446:                   }
1447:                 }
1448:               }
1449:               if (flipP) {
1450:                 for (r = 0; r < nWork; r++) {
1451:                   for (s = 0; s < nWorkP; s++) {
1452:                     scwork[r * nWorkP + s] *= flipP[s];
1453:                   }
1454:                 }
1455:               }
1456:               MatSetValues(cMat,nWork,workIndRow,nWorkP,workIndCol,scwork,INSERT_VALUES);
1457:               break;
1458:             }
1459:           }
1460:         }
1461:       }
1462:       MatDenseRestoreArray(Xmat,&X);
1463:       PetscFree2(childOffsets,parentOffsets);
1464:       DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1465:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1466:     }
1467:     MatDestroy(&Amat);
1468:     MatDestroy(&Bmat);
1469:     MatDestroy(&Xmat);
1470:     PetscFree(scwork);
1471:     PetscFree7(sizes,weights,pointsRef,pointsReal,work,workIndRow,workIndCol);
1472:     if (id == PETSCFV_CLASSID) {
1473:       PetscSpaceDestroy(&bspace);
1474:     }
1475:   }
1476:   MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1477:   MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1478:   PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);
1479:   ISRestoreIndices(aIS,&anchors);

1481:   return(0);
1482: }

1484: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1485: {
1486:   Mat               refCmat;
1487:   PetscDS           ds;
1488:   PetscInt          numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1489:   PetscScalar       ***refPointFieldMats;
1490:   PetscSection      refConSec, refAnSec, refSection;
1491:   IS                refAnIS;
1492:   const PetscInt    *refAnchors;
1493:   const PetscInt    **perms;
1494:   const PetscScalar **flips;
1495:   PetscErrorCode    ierr;

1498:   DMGetDS(refTree,&ds);
1499:   PetscDSGetNumFields(ds,&numFields);
1500:   maxFields = PetscMax(1,numFields);
1501:   DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1502:   DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1503:   ISGetIndices(refAnIS,&refAnchors);
1504:   DMGetSection(refTree,&refSection);
1505:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1506:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
1507:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);
1508:   PetscSectionGetMaxDof(refConSec,&maxDof);
1509:   PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1510:   PetscMalloc1(maxDof,&rows);
1511:   PetscMalloc1(maxDof*maxAnDof,&cols);
1512:   for (p = pRefStart; p < pRefEnd; p++) {
1513:     PetscInt parent, closureSize, *closure = NULL, pDof;

1515:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
1516:     PetscSectionGetDof(refConSec,p,&pDof);
1517:     if (!pDof || parent == p) continue;

1519:     PetscMalloc1(maxFields,&refPointFieldMats[p-pRefStart]);
1520:     PetscCalloc1(maxFields,&refPointFieldN[p-pRefStart]);
1521:     DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1522:     for (f = 0; f < maxFields; f++) {
1523:       PetscInt cDof, cOff, numCols, r, i;

1525:       if (f < numFields) {
1526:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1527:         PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
1528:         PetscSectionGetFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);
1529:       } else {
1530:         PetscSectionGetDof(refConSec,p,&cDof);
1531:         PetscSectionGetOffset(refConSec,p,&cOff);
1532:         PetscSectionGetPointSyms(refSection,closureSize,closure,&perms,&flips);
1533:       }

1535:       for (r = 0; r < cDof; r++) {
1536:         rows[r] = cOff + r;
1537:       }
1538:       numCols = 0;
1539:       for (i = 0; i < closureSize; i++) {
1540:         PetscInt          q = closure[2*i];
1541:         PetscInt          aDof, aOff, j;
1542:         const PetscInt    *perm = perms ? perms[i] : NULL;

1544:         if (numFields) {
1545:           PetscSectionGetFieldDof(refSection,q,f,&aDof);
1546:           PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1547:         }
1548:         else {
1549:           PetscSectionGetDof(refSection,q,&aDof);
1550:           PetscSectionGetOffset(refSection,q,&aOff);
1551:         }

1553:         for (j = 0; j < aDof; j++) {
1554:           cols[numCols++] = aOff + (perm ? perm[j] : j);
1555:         }
1556:       }
1557:       refPointFieldN[p-pRefStart][f] = numCols;
1558:       PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
1559:       MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);
1560:       if (flips) {
1561:         PetscInt colOff = 0;

1563:         for (i = 0; i < closureSize; i++) {
1564:           PetscInt          q = closure[2*i];
1565:           PetscInt          aDof, aOff, j;
1566:           const PetscScalar *flip = flips ? flips[i] : NULL;

1568:           if (numFields) {
1569:             PetscSectionGetFieldDof(refSection,q,f,&aDof);
1570:             PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1571:           }
1572:           else {
1573:             PetscSectionGetDof(refSection,q,&aDof);
1574:             PetscSectionGetOffset(refSection,q,&aOff);
1575:           }
1576:           if (flip) {
1577:             PetscInt k;
1578:             for (k = 0; k < cDof; k++) {
1579:               for (j = 0; j < aDof; j++) {
1580:                 refPointFieldMats[p-pRefStart][f][k * numCols + colOff + j] *= flip[j];
1581:               }
1582:             }
1583:           }
1584:           colOff += aDof;
1585:         }
1586:       }
1587:       if (numFields) {
1588:         PetscSectionRestoreFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);
1589:       } else {
1590:         PetscSectionRestorePointSyms(refSection,closureSize,closure,&perms,&flips);
1591:       }
1592:     }
1593:     DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1594:   }
1595:   *childrenMats = refPointFieldMats;
1596:   *childrenN = refPointFieldN;
1597:   ISRestoreIndices(refAnIS,&refAnchors);
1598:   PetscFree(rows);
1599:   PetscFree(cols);
1600:   return(0);
1601: }

1603: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1604: {
1605:   PetscDS        ds;
1606:   PetscInt       **refPointFieldN;
1607:   PetscScalar    ***refPointFieldMats;
1608:   PetscInt       numFields, maxFields, pRefStart, pRefEnd, p, f;
1609:   PetscSection   refConSec;

1613:   refPointFieldN = *childrenN;
1614:   *childrenN = NULL;
1615:   refPointFieldMats = *childrenMats;
1616:   *childrenMats = NULL;
1617:   DMGetDS(refTree,&ds);
1618:   PetscDSGetNumFields(ds,&numFields);
1619:   maxFields = PetscMax(1,numFields);
1620:   DMGetDefaultConstraints(refTree,&refConSec,NULL);
1621:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1622:   for (p = pRefStart; p < pRefEnd; p++) {
1623:     PetscInt parent, pDof;

1625:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
1626:     PetscSectionGetDof(refConSec,p,&pDof);
1627:     if (!pDof || parent == p) continue;

1629:     for (f = 0; f < maxFields; f++) {
1630:       PetscInt cDof;

1632:       if (numFields) {
1633:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1634:       }
1635:       else {
1636:         PetscSectionGetDof(refConSec,p,&cDof);
1637:       }

1639:       PetscFree(refPointFieldMats[p - pRefStart][f]);
1640:     }
1641:     PetscFree(refPointFieldMats[p - pRefStart]);
1642:     PetscFree(refPointFieldN[p - pRefStart]);
1643:   }
1644:   PetscFree(refPointFieldMats);
1645:   PetscFree(refPointFieldN);
1646:   return(0);
1647: }

1649: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1650: {
1651:   DM             refTree;
1652:   PetscDS        ds;
1653:   Mat            refCmat;
1654:   PetscInt       numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1655:   PetscScalar ***refPointFieldMats, *pointWork;
1656:   PetscSection   refConSec, refAnSec, anSec;
1657:   IS             refAnIS, anIS;
1658:   const PetscInt *anchors;

1663:   DMGetDS(dm,&ds);
1664:   PetscDSGetNumFields(ds,&numFields);
1665:   maxFields = PetscMax(1,numFields);
1666:   DMPlexGetReferenceTree(dm,&refTree);
1667:   DMCopyDisc(dm,refTree);
1668:   DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1669:   DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1670:   DMPlexGetAnchors(dm,&anSec,&anIS);
1671:   ISGetIndices(anIS,&anchors);
1672:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1673:   PetscSectionGetChart(conSec,&conStart,&conEnd);
1674:   PetscSectionGetMaxDof(refConSec,&maxDof);
1675:   PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1676:   PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);

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

1681:   /* step 2: compute the preorder */
1682:   DMPlexGetChart(dm,&pStart,&pEnd);
1683:   PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);
1684:   for (p = pStart; p < pEnd; p++) {
1685:     perm[p - pStart] = p;
1686:     iperm[p - pStart] = p-pStart;
1687:   }
1688:   for (p = 0; p < pEnd - pStart;) {
1689:     PetscInt point = perm[p];
1690:     PetscInt parent;

1692:     DMPlexGetTreeParent(dm,point,&parent,NULL);
1693:     if (parent == point) {
1694:       p++;
1695:     }
1696:     else {
1697:       PetscInt size, closureSize, *closure = NULL, i;

1699:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1700:       for (i = 0; i < closureSize; i++) {
1701:         PetscInt q = closure[2*i];
1702:         if (iperm[q-pStart] > iperm[point-pStart]) {
1703:           /* swap */
1704:           perm[p]               = q;
1705:           perm[iperm[q-pStart]] = point;
1706:           iperm[point-pStart]   = iperm[q-pStart];
1707:           iperm[q-pStart]       = p;
1708:           break;
1709:         }
1710:       }
1711:       size = closureSize;
1712:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1713:       if (i == size) {
1714:         p++;
1715:       }
1716:     }
1717:   }

1719:   /* step 3: fill the constraint matrix */
1720:   /* we are going to use a preorder progressive fill strategy.  Mat doesn't
1721:    * allow progressive fill without assembly, so we are going to set up the
1722:    * values outside of the Mat first.
1723:    */
1724:   {
1725:     PetscInt nRows, row, nnz;
1726:     PetscBool done;
1727:     const PetscInt *ia, *ja;
1728:     PetscScalar *vals;

1730:     MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1731:     if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix");
1732:     nnz  = ia[nRows];
1733:     /* malloc and then zero rows right before we fill them: this way valgrind
1734:      * can tell if we are doing progressive fill in the wrong order */
1735:     PetscMalloc1(nnz,&vals);
1736:     for (p = 0; p < pEnd - pStart; p++) {
1737:       PetscInt        parent, childid, closureSize, *closure = NULL;
1738:       PetscInt        point = perm[p], pointDof;

1740:       DMPlexGetTreeParent(dm,point,&parent,&childid);
1741:       if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1742:       PetscSectionGetDof(conSec,point,&pointDof);
1743:       if (!pointDof) continue;
1744:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1745:       for (f = 0; f < maxFields; f++) {
1746:         PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
1747:         PetscScalar *pointMat;
1748:         const PetscInt    **perms;
1749:         const PetscScalar **flips;

1751:         if (numFields) {
1752:           PetscSectionGetFieldDof(conSec,point,f,&cDof);
1753:           PetscSectionGetFieldOffset(conSec,point,f,&cOff);
1754:         }
1755:         else {
1756:           PetscSectionGetDof(conSec,point,&cDof);
1757:           PetscSectionGetOffset(conSec,point,&cOff);
1758:         }
1759:         if (!cDof) continue;
1760:         if (numFields) {PetscSectionGetFieldPointSyms(section,f,closureSize,closure,&perms,&flips);}
1761:         else           {PetscSectionGetPointSyms(section,closureSize,closure,&perms,&flips);}

1763:         /* make sure that every row for this point is the same size */
1764: #if defined(PETSC_USE_DEBUG)
1765:         for (r = 0; r < cDof; r++) {
1766:           if (cDof > 1 && r) {
1767:             if ((ia[cOff+r+1]-ia[cOff+r]) != (ia[cOff+r]-ia[cOff+r-1])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Two point rows have different nnz: %D vs. %D", (ia[cOff+r+1]-ia[cOff+r]), (ia[cOff+r]-ia[cOff+r-1]));
1768:           }
1769:         }
1770: #endif
1771:         /* zero rows */
1772:         for (i = ia[cOff] ; i< ia[cOff+cDof];i++) {
1773:           vals[i] = 0.;
1774:         }
1775:         matOffset = ia[cOff];
1776:         numFillCols = ia[cOff+1] - matOffset;
1777:         pointMat = refPointFieldMats[childid-pRefStart][f];
1778:         numCols = refPointFieldN[childid-pRefStart][f];
1779:         offset = 0;
1780:         for (i = 0; i < closureSize; i++) {
1781:           PetscInt q = closure[2*i];
1782:           PetscInt aDof, aOff, j, k, qConDof, qConOff;
1783:           const PetscInt    *perm = perms ? perms[i] : NULL;
1784:           const PetscScalar *flip = flips ? flips[i] : NULL;

1786:           qConDof = qConOff = 0;
1787:           if (numFields) {
1788:             PetscSectionGetFieldDof(section,q,f,&aDof);
1789:             PetscSectionGetFieldOffset(section,q,f,&aOff);
1790:             if (q >= conStart && q < conEnd) {
1791:               PetscSectionGetFieldDof(conSec,q,f,&qConDof);
1792:               PetscSectionGetFieldOffset(conSec,q,f,&qConOff);
1793:             }
1794:           }
1795:           else {
1796:             PetscSectionGetDof(section,q,&aDof);
1797:             PetscSectionGetOffset(section,q,&aOff);
1798:             if (q >= conStart && q < conEnd) {
1799:               PetscSectionGetDof(conSec,q,&qConDof);
1800:               PetscSectionGetOffset(conSec,q,&qConOff);
1801:             }
1802:           }
1803:           if (!aDof) continue;
1804:           if (qConDof) {
1805:             /* this point has anchors: its rows of the matrix should already
1806:              * be filled, thanks to preordering */
1807:             /* first multiply into pointWork, then set in matrix */
1808:             PetscInt aMatOffset = ia[qConOff];
1809:             PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1810:             for (r = 0; r < cDof; r++) {
1811:               for (j = 0; j < aNumFillCols; j++) {
1812:                 PetscScalar inVal = 0;
1813:                 for (k = 0; k < aDof; k++) {
1814:                   PetscInt col = perm ? perm[k] : k;

1816:                   inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.);
1817:                 }
1818:                 pointWork[r * aNumFillCols + j] = inVal;
1819:               }
1820:             }
1821:             /* assume that the columns are sorted, spend less time searching */
1822:             for (j = 0, k = 0; j < aNumFillCols; j++) {
1823:               PetscInt col = ja[aMatOffset + j];
1824:               for (;k < numFillCols; k++) {
1825:                 if (ja[matOffset + k] == col) {
1826:                   break;
1827:                 }
1828:               }
1829:               if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col);
1830:               for (r = 0; r < cDof; r++) {
1831:                 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1832:               }
1833:             }
1834:           }
1835:           else {
1836:             /* find where to put this portion of pointMat into the matrix */
1837:             for (k = 0; k < numFillCols; k++) {
1838:               if (ja[matOffset + k] == aOff) {
1839:                 break;
1840:               }
1841:             }
1842:             if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff);
1843:             for (r = 0; r < cDof; r++) {
1844:               for (j = 0; j < aDof; j++) {
1845:                 PetscInt col = perm ? perm[j] : j;

1847:                 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
1848:               }
1849:             }
1850:           }
1851:           offset += aDof;
1852:         }
1853:         if (numFields) {
1854:           PetscSectionRestoreFieldPointSyms(section,f,closureSize,closure,&perms,&flips);
1855:         } else {
1856:           PetscSectionRestorePointSyms(section,closureSize,closure,&perms,&flips);
1857:         }
1858:       }
1859:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1860:     }
1861:     for (row = 0; row < nRows; row++) {
1862:       MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);
1863:     }
1864:     MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1865:     if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix");
1866:     MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1867:     MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1868:     PetscFree(vals);
1869:   }

1871:   /* clean up */
1872:   ISRestoreIndices(anIS,&anchors);
1873:   PetscFree2(perm,iperm);
1874:   PetscFree(pointWork);
1875:   DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
1876:   return(0);
1877: }

1879: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1880:  * a non-conforming mesh.  Local refinement comes later */
1881: PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm)
1882: {
1883:   DM K;
1884:   PetscMPIInt rank;
1885:   PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1886:   PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1887:   PetscInt *Kembedding;
1888:   PetscInt *cellClosure=NULL, nc;
1889:   PetscScalar *newVertexCoords;
1890:   PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1891:   PetscSection parentSection;

1895:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1896:   DMGetDimension(dm,&dim);
1897:   DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);
1898:   DMSetDimension(*ncdm,dim);

1900:   DMPlexGetChart(dm, &pStart, &pEnd);
1901:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);
1902:   DMPlexGetReferenceTree(dm,&K);
1903:   if (!rank) {
1904:     /* compute the new charts */
1905:     PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);
1906:     offset = 0;
1907:     for (d = 0; d <= dim; d++) {
1908:       PetscInt pOldCount, kStart, kEnd, k;

1910:       pNewStart[d] = offset;
1911:       DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);
1912:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1913:       pOldCount = pOldEnd[d] - pOldStart[d];
1914:       /* adding the new points */
1915:       pNewCount[d] = pOldCount + kEnd - kStart;
1916:       if (!d) {
1917:         /* removing the cell */
1918:         pNewCount[d]--;
1919:       }
1920:       for (k = kStart; k < kEnd; k++) {
1921:         PetscInt parent;
1922:         DMPlexGetTreeParent(K,k,&parent,NULL);
1923:         if (parent == k) {
1924:           /* avoid double counting points that won't actually be new */
1925:           pNewCount[d]--;
1926:         }
1927:       }
1928:       pNewEnd[d] = pNewStart[d] + pNewCount[d];
1929:       offset = pNewEnd[d];

1931:     }
1932:     if (cell < pOldStart[0] || cell >= pOldEnd[0]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"%d not in cell range [%d, %d)", cell, pOldStart[0], pOldEnd[0]);
1933:     /* get the current closure of the cell that we are removing */
1934:     DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);

1936:     PetscMalloc1(pNewEnd[dim],&newConeSizes);
1937:     {
1938:       PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;

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

1943:       for (k = kStart; k < kEnd; k++) {
1944:         perm[k - kStart] = k;
1945:         iperm [k - kStart] = k - kStart;
1946:         preOrient[k - kStart] = 0;
1947:       }

1949:       DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1950:       for (j = 1; j < closureSizeK; j++) {
1951:         PetscInt parentOrientA = closureK[2*j+1];
1952:         PetscInt parentOrientB = cellClosure[2*j+1];
1953:         PetscInt p, q;

1955:         p = closureK[2*j];
1956:         q = cellClosure[2*j];
1957:         for (d = 0; d <= dim; d++) {
1958:           if (q >= pOldStart[d] && q < pOldEnd[d]) {
1959:             Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1960:           }
1961:         }
1962:         if (parentOrientA != parentOrientB) {
1963:           PetscInt numChildren, i;
1964:           const PetscInt *children;

1966:           DMPlexGetTreeChildren(K,p,&numChildren,&children);
1967:           for (i = 0; i < numChildren; i++) {
1968:             PetscInt kPerm, oPerm;

1970:             k    = children[i];
1971:             DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);
1972:             /* perm = what refTree position I'm in */
1973:             perm[kPerm-kStart]      = k;
1974:             /* iperm = who is at this position */
1975:             iperm[k-kStart]         = kPerm-kStart;
1976:             preOrient[kPerm-kStart] = oPerm;
1977:           }
1978:         }
1979:       }
1980:       DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1981:     }
1982:     PetscSectionSetChart(parentSection,0,pNewEnd[dim]);
1983:     offset = 0;
1984:     numNewCones = 0;
1985:     for (d = 0; d <= dim; d++) {
1986:       PetscInt kStart, kEnd, k;
1987:       PetscInt p;
1988:       PetscInt size;

1990:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1991:         /* skip cell 0 */
1992:         if (p == cell) continue;
1993:         /* old cones to new cones */
1994:         DMPlexGetConeSize(dm,p,&size);
1995:         newConeSizes[offset++] = size;
1996:         numNewCones += size;
1997:       }

1999:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
2000:       for (k = kStart; k < kEnd; k++) {
2001:         PetscInt kParent;

2003:         DMPlexGetTreeParent(K,k,&kParent,NULL);
2004:         if (kParent != k) {
2005:           Kembedding[k] = offset;
2006:           DMPlexGetConeSize(K,k,&size);
2007:           newConeSizes[offset++] = size;
2008:           numNewCones += size;
2009:           if (kParent != 0) {
2010:             PetscSectionSetDof(parentSection,Kembedding[k],1);
2011:           }
2012:         }
2013:       }
2014:     }

2016:     PetscSectionSetUp(parentSection);
2017:     PetscSectionGetStorageSize(parentSection,&numPointsWithParents);
2018:     PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);
2019:     PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);

2021:     /* fill new cones */
2022:     offset = 0;
2023:     for (d = 0; d <= dim; d++) {
2024:       PetscInt kStart, kEnd, k, l;
2025:       PetscInt p;
2026:       PetscInt size;
2027:       const PetscInt *cone, *orientation;

2029:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
2030:         /* skip cell 0 */
2031:         if (p == cell) continue;
2032:         /* old cones to new cones */
2033:         DMPlexGetConeSize(dm,p,&size);
2034:         DMPlexGetCone(dm,p,&cone);
2035:         DMPlexGetConeOrientation(dm,p,&orientation);
2036:         for (l = 0; l < size; l++) {
2037:           newCones[offset]          = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
2038:           newOrientations[offset++] = orientation[l];
2039:         }
2040:       }

2042:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
2043:       for (k = kStart; k < kEnd; k++) {
2044:         PetscInt kPerm = perm[k], kParent;
2045:         PetscInt preO  = preOrient[k];

2047:         DMPlexGetTreeParent(K,k,&kParent,NULL);
2048:         if (kParent != k) {
2049:           /* embed new cones */
2050:           DMPlexGetConeSize(K,k,&size);
2051:           DMPlexGetCone(K,kPerm,&cone);
2052:           DMPlexGetConeOrientation(K,kPerm,&orientation);
2053:           for (l = 0; l < size; l++) {
2054:             PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size);
2055:             PetscInt newO, lSize, oTrue;

2057:             q                         = iperm[cone[m]];
2058:             newCones[offset]          = Kembedding[q];
2059:             DMPlexGetConeSize(K,q,&lSize);
2060:             oTrue                     = orientation[m];
2061:             oTrue                     = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
2062:             newO                      = DihedralCompose(lSize,oTrue,preOrient[q]);
2063:             newOrientations[offset++] = newO;
2064:           }
2065:           if (kParent != 0) {
2066:             PetscInt newPoint = Kembedding[kParent];
2067:             PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);
2068:             parents[pOffset]  = newPoint;
2069:             childIDs[pOffset] = k;
2070:           }
2071:         }
2072:       }
2073:     }

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

2077:     /* fill coordinates */
2078:     offset = 0;
2079:     {
2080:       PetscInt kStart, kEnd, l;
2081:       PetscSection vSection;
2082:       PetscInt v;
2083:       Vec coords;
2084:       PetscScalar *coordvals;
2085:       PetscInt dof, off;
2086:       PetscReal v0[3], J[9], detJ;

2088: #if defined(PETSC_USE_DEBUG)
2089:       {
2090:         PetscInt k;
2091:         DMPlexGetHeightStratum(K,0,&kStart,&kEnd);
2092:         for (k = kStart; k < kEnd; k++) {
2093:           DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);
2094:           if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k);
2095:         }
2096:       }
2097: #endif
2098:       DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);
2099:       DMGetCoordinateSection(dm,&vSection);
2100:       DMGetCoordinatesLocal(dm,&coords);
2101:       VecGetArray(coords,&coordvals);
2102:       for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {

2104:         PetscSectionGetDof(vSection,v,&dof);
2105:         PetscSectionGetOffset(vSection,v,&off);
2106:         for (l = 0; l < dof; l++) {
2107:           newVertexCoords[offset++] = coordvals[off + l];
2108:         }
2109:       }
2110:       VecRestoreArray(coords,&coordvals);

2112:       DMGetCoordinateSection(K,&vSection);
2113:       DMGetCoordinatesLocal(K,&coords);
2114:       VecGetArray(coords,&coordvals);
2115:       DMPlexGetDepthStratum(K,0,&kStart,&kEnd);
2116:       for (v = kStart; v < kEnd; v++) {
2117:         PetscReal coord[3], newCoord[3];
2118:         PetscInt  vPerm = perm[v];
2119:         PetscInt  kParent;
2120:         const PetscReal xi0[3] = {-1.,-1.,-1.};

2122:         DMPlexGetTreeParent(K,v,&kParent,NULL);
2123:         if (kParent != v) {
2124:           /* this is a new vertex */
2125:           PetscSectionGetOffset(vSection,vPerm,&off);
2126:           for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]);
2127:           CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
2128:           for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l];
2129:           offset += dim;
2130:         }
2131:       }
2132:       VecRestoreArray(coords,&coordvals);
2133:     }

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

2139:       tmp = pNewCount[d];
2140:       pNewCount[d] = pNewCount[dim - d];
2141:       pNewCount[dim - d] = tmp;
2142:     }

2144:     DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);
2145:     DMPlexSetReferenceTree(*ncdm,K);
2146:     DMPlexSetTree(*ncdm,parentSection,parents,childIDs);

2148:     /* clean up */
2149:     DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
2150:     PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);
2151:     PetscFree(newConeSizes);
2152:     PetscFree2(newCones,newOrientations);
2153:     PetscFree(newVertexCoords);
2154:     PetscFree2(parents,childIDs);
2155:     PetscFree4(Kembedding,perm,iperm,preOrient);
2156:   }
2157:   else {
2158:     PetscInt    p, counts[4];
2159:     PetscInt    *coneSizes, *cones, *orientations;
2160:     Vec         coordVec;
2161:     PetscScalar *coords;

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

2166:       DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
2167:       counts[d] = dEnd - dStart;
2168:     }
2169:     PetscMalloc1(pEnd-pStart,&coneSizes);
2170:     for (p = pStart; p < pEnd; p++) {
2171:       DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);
2172:     }
2173:     DMPlexGetCones(dm, &cones);
2174:     DMPlexGetConeOrientations(dm, &orientations);
2175:     DMGetCoordinatesLocal(dm,&coordVec);
2176:     VecGetArray(coordVec,&coords);

2178:     PetscSectionSetChart(parentSection,pStart,pEnd);
2179:     PetscSectionSetUp(parentSection);
2180:     DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);
2181:     DMPlexSetReferenceTree(*ncdm,K);
2182:     DMPlexSetTree(*ncdm,parentSection,NULL,NULL);
2183:     VecRestoreArray(coordVec,&coords);
2184:   }
2185:   PetscSectionDestroy(&parentSection);

2187:   return(0);
2188: }

2190: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2191: {
2192:   PetscSF           coarseToFineEmbedded;
2193:   PetscSection      globalCoarse, globalFine;
2194:   PetscSection      localCoarse, localFine;
2195:   PetscSection      aSec, cSec;
2196:   PetscSection      rootIndicesSec, rootMatricesSec;
2197:   PetscSection      leafIndicesSec, leafMatricesSec;
2198:   PetscInt          *rootIndices, *leafIndices;
2199:   PetscScalar       *rootMatrices, *leafMatrices;
2200:   IS                aIS;
2201:   const PetscInt    *anchors;
2202:   Mat               cMat;
2203:   PetscInt          numFields, maxFields;
2204:   PetscInt          pStartC, pEndC, pStartF, pEndF, p;
2205:   PetscInt          aStart, aEnd, cStart, cEnd;
2206:   PetscInt          *maxChildIds;
2207:   PetscInt          *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2208:   const PetscInt    ***perms;
2209:   const PetscScalar ***flips;
2210:   PetscErrorCode    ierr;

2213:   DMPlexGetChart(coarse,&pStartC,&pEndC);
2214:   DMPlexGetChart(fine,&pStartF,&pEndF);
2215:   DMGetGlobalSection(fine,&globalFine);
2216:   { /* winnow fine points that don't have global dofs out of the sf */
2217:     PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2218:     const PetscInt *leaves;

2220:     PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);
2221:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2222:       p = leaves ? leaves[l] : l;
2223:       PetscSectionGetDof(globalFine,p,&dof);
2224:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
2225:       if ((dof - cdof) > 0) {
2226:         numPointsWithDofs++;
2227:       }
2228:     }
2229:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
2230:     for (l = 0, offset = 0; l < nleaves; l++) {
2231:       p = leaves ? leaves[l] : l;
2232:       PetscSectionGetDof(globalFine,p,&dof);
2233:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
2234:       if ((dof - cdof) > 0) {
2235:         pointsWithDofs[offset++] = l;
2236:       }
2237:     }
2238:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
2239:     PetscFree(pointsWithDofs);
2240:   }
2241:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2242:   PetscMalloc1(pEndC-pStartC,&maxChildIds);
2243:   for (p = pStartC; p < pEndC; p++) {
2244:     maxChildIds[p - pStartC] = -2;
2245:   }
2246:   PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);
2247:   PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);

2249:   DMGetSection(coarse,&localCoarse);
2250:   DMGetGlobalSection(coarse,&globalCoarse);

2252:   DMPlexGetAnchors(coarse,&aSec,&aIS);
2253:   ISGetIndices(aIS,&anchors);
2254:   PetscSectionGetChart(aSec,&aStart,&aEnd);

2256:   DMGetDefaultConstraints(coarse,&cSec,&cMat);
2257:   PetscSectionGetChart(cSec,&cStart,&cEnd);

2259:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2260:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
2261:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);
2262:   PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);
2263:   PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);
2264:   PetscSectionGetNumFields(localCoarse,&numFields);
2265:   maxFields = PetscMax(1,numFields);
2266:   PetscMalloc7(maxFields+1,&offsets,maxFields+1,&offsetsCopy,maxFields+1,&newOffsets,maxFields+1,&newOffsetsCopy,maxFields+1,&rowOffsets,maxFields+1,&numD,maxFields+1,&numO);
2267:   PetscMalloc2(maxFields+1,(PetscInt****)&perms,maxFields+1,(PetscScalar****)&flips);
2268:   PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **));
2269:   PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **));

2271:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2272:     PetscInt dof, matSize   = 0;
2273:     PetscInt aDof           = 0;
2274:     PetscInt cDof           = 0;
2275:     PetscInt maxChildId     = maxChildIds[p - pStartC];
2276:     PetscInt numRowIndices  = 0;
2277:     PetscInt numColIndices  = 0;
2278:     PetscInt f;

2280:     PetscSectionGetDof(globalCoarse,p,&dof);
2281:     if (dof < 0) {
2282:       dof = -(dof + 1);
2283:     }
2284:     if (p >= aStart && p < aEnd) {
2285:       PetscSectionGetDof(aSec,p,&aDof);
2286:     }
2287:     if (p >= cStart && p < cEnd) {
2288:       PetscSectionGetDof(cSec,p,&cDof);
2289:     }
2290:     for (f = 0; f <= numFields; f++) offsets[f] = 0;
2291:     for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2292:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2293:       PetscInt *closure = NULL, closureSize, cl;

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

2299:         PetscSectionGetDof(localCoarse,c,&clDof);
2300:         numRowIndices += clDof;
2301:         for (f = 0; f < numFields; f++) {
2302:           PetscSectionGetFieldDof(localCoarse,c,f,&clDof);
2303:           offsets[f + 1] += clDof;
2304:         }
2305:       }
2306:       for (f = 0; f < numFields; f++) {
2307:         offsets[f + 1]   += offsets[f];
2308:         newOffsets[f + 1] = offsets[f + 1];
2309:       }
2310:       /* get the number of indices needed and their field offsets */
2311:       DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);
2312:       DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
2313:       if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2314:         numColIndices = numRowIndices;
2315:         matSize = 0;
2316:       }
2317:       else if (numFields) { /* we send one submat for each field: sum their sizes */
2318:         matSize = 0;
2319:         for (f = 0; f < numFields; f++) {
2320:           PetscInt numRow, numCol;

2322:           numRow = offsets[f + 1] - offsets[f];
2323:           numCol = newOffsets[f + 1] - newOffsets[f];
2324:           matSize += numRow * numCol;
2325:         }
2326:       }
2327:       else {
2328:         matSize = numRowIndices * numColIndices;
2329:       }
2330:     } else if (maxChildId == -1) {
2331:       if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2332:         PetscInt aOff, a;

2334:         PetscSectionGetOffset(aSec,p,&aOff);
2335:         for (f = 0; f < numFields; f++) {
2336:           PetscInt fDof;

2338:           PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2339:           offsets[f+1] = fDof;
2340:         }
2341:         for (a = 0; a < aDof; a++) {
2342:           PetscInt anchor = anchors[a + aOff], aLocalDof;

2344:           PetscSectionGetDof(localCoarse,anchor,&aLocalDof);
2345:           numColIndices += aLocalDof;
2346:           for (f = 0; f < numFields; f++) {
2347:             PetscInt fDof;

2349:             PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2350:             newOffsets[f+1] += fDof;
2351:           }
2352:         }
2353:         if (numFields) {
2354:           matSize = 0;
2355:           for (f = 0; f < numFields; f++) {
2356:             matSize += offsets[f+1] * newOffsets[f+1];
2357:           }
2358:         }
2359:         else {
2360:           matSize = numColIndices * dof;
2361:         }
2362:       }
2363:       else { /* no children, and no constraints on dofs: just get the global indices */
2364:         numColIndices = dof;
2365:         matSize       = 0;
2366:       }
2367:     }
2368:     /* we will pack the column indices with the field offsets */
2369:     PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);
2370:     PetscSectionSetDof(rootMatricesSec,p,matSize);
2371:   }
2372:   PetscSectionSetUp(rootIndicesSec);
2373:   PetscSectionSetUp(rootMatricesSec);
2374:   {
2375:     PetscInt numRootIndices, numRootMatrices;

2377:     PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
2378:     PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);
2379:     PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);
2380:     for (p = pStartC; p < pEndC; p++) {
2381:       PetscInt    numRowIndices, numColIndices, matSize, dof;
2382:       PetscInt    pIndOff, pMatOff, f;
2383:       PetscInt    *pInd;
2384:       PetscInt    maxChildId = maxChildIds[p - pStartC];
2385:       PetscScalar *pMat = NULL;

2387:       PetscSectionGetDof(rootIndicesSec,p,&numColIndices);
2388:       if (!numColIndices) {
2389:         continue;
2390:       }
2391:       for (f = 0; f <= numFields; f++) {
2392:         offsets[f]        = 0;
2393:         newOffsets[f]     = 0;
2394:         offsetsCopy[f]    = 0;
2395:         newOffsetsCopy[f] = 0;
2396:       }
2397:       numColIndices -= 2 * numFields;
2398:       PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);
2399:       pInd = &(rootIndices[pIndOff]);
2400:       PetscSectionGetDof(rootMatricesSec,p,&matSize);
2401:       if (matSize) {
2402:         PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);
2403:         pMat = &rootMatrices[pMatOff];
2404:       }
2405:       PetscSectionGetDof(globalCoarse,p,&dof);
2406:       if (dof < 0) {
2407:         dof = -(dof + 1);
2408:       }
2409:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2410:         PetscInt i, j;
2411:         PetscInt numRowIndices = matSize / numColIndices;

2413:         if (!numRowIndices) { /* don't need to calculate the mat, just the indices */
2414:           PetscInt numIndices, *indices;
2415:           DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2416:           if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations");
2417:           for (i = 0; i < numColIndices; i++) {
2418:             pInd[i] = indices[i];
2419:           }
2420:           for (i = 0; i < numFields; i++) {
2421:             pInd[numColIndices + i]             = offsets[i+1];
2422:             pInd[numColIndices + numFields + i] = offsets[i+1];
2423:           }
2424:           DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2425:         }
2426:         else {
2427:           PetscInt closureSize, *closure = NULL, cl;
2428:           PetscScalar *pMatIn, *pMatModified;
2429:           PetscInt numPoints,*points;

2431:           DMGetWorkArray(coarse,numRowIndices * numRowIndices,MPIU_SCALAR,&pMatIn);
2432:           for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2433:             for (j = 0; j < numRowIndices; j++) {
2434:               pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2435:             }
2436:           }
2437:           DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2438:           for (f = 0; f < maxFields; f++) {
2439:             if (numFields) {PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);}
2440:             else           {PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);}
2441:           }
2442:           if (numFields) {
2443:             for (cl = 0; cl < closureSize; cl++) {
2444:               PetscInt c = closure[2 * cl];

2446:               for (f = 0; f < numFields; f++) {
2447:                 PetscInt fDof;

2449:                 PetscSectionGetFieldDof(localCoarse,c,f,&fDof);
2450:                 offsets[f + 1] += fDof;
2451:               }
2452:             }
2453:             for (f = 0; f < numFields; f++) {
2454:               offsets[f + 1]   += offsets[f];
2455:               newOffsets[f + 1] = offsets[f + 1];
2456:             }
2457:           }
2458:           /* TODO : flips here ? */
2459:           /* apply hanging node constraints on the right, get the new points and the new offsets */
2460:           DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);
2461:           for (f = 0; f < maxFields; f++) {
2462:             if (numFields) {PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);}
2463:             else           {PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);}
2464:           }
2465:           for (f = 0; f < maxFields; f++) {
2466:             if (numFields) {PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);}
2467:             else           {PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);}
2468:           }
2469:           if (!numFields) {
2470:             for (i = 0; i < numRowIndices * numColIndices; i++) {
2471:               pMat[i] = pMatModified[i];
2472:             }
2473:           }
2474:           else {
2475:             PetscInt i, j, count;
2476:             for (f = 0, count = 0; f < numFields; f++) {
2477:               for (i = offsets[f]; i < offsets[f+1]; i++) {
2478:                 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) {
2479:                   pMat[count] = pMatModified[i * numColIndices + j];
2480:                 }
2481:               }
2482:             }
2483:           }
2484:           DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatModified);
2485:           DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2486:           DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatIn);
2487:           if (numFields) {
2488:             for (f = 0; f < numFields; f++) {
2489:               pInd[numColIndices + f]             = offsets[f+1];
2490:               pInd[numColIndices + numFields + f] = newOffsets[f+1];
2491:             }
2492:             for (cl = 0; cl < numPoints; cl++) {
2493:               PetscInt globalOff, c = points[2*cl];
2494:               PetscSectionGetOffset(globalCoarse, c, &globalOff);
2495:               DMPlexGetIndicesPointFields_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, pInd);
2496:             }
2497:           } else {
2498:             for (cl = 0; cl < numPoints; cl++) {
2499:               PetscInt c = points[2*cl], globalOff;
2500:               const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;

2502:               PetscSectionGetOffset(globalCoarse, c, &globalOff);
2503:               DMPlexGetIndicesPoint_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, pInd);
2504:             }
2505:           }
2506:           for (f = 0; f < maxFields; f++) {
2507:             if (numFields) {PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);}
2508:             else           {PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);}
2509:           }
2510:           DMRestoreWorkArray(coarse,numPoints,MPIU_SCALAR,&points);
2511:         }
2512:       }
2513:       else if (matSize) {
2514:         PetscInt cOff;
2515:         PetscInt *rowIndices, *colIndices, a, aDof, aOff;

2517:         numRowIndices = matSize / numColIndices;
2518:         if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs");
2519:         DMGetWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);
2520:         DMGetWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);
2521:         PetscSectionGetOffset(cSec,p,&cOff);
2522:         PetscSectionGetDof(aSec,p,&aDof);
2523:         PetscSectionGetOffset(aSec,p,&aOff);
2524:         if (numFields) {
2525:           for (f = 0; f < numFields; f++) {
2526:             PetscInt fDof;

2528:             PetscSectionGetFieldDof(cSec,p,f,&fDof);
2529:             offsets[f + 1] = fDof;
2530:             for (a = 0; a < aDof; a++) {
2531:               PetscInt anchor = anchors[a + aOff];
2532:               PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2533:               newOffsets[f + 1] += fDof;
2534:             }
2535:           }
2536:           for (f = 0; f < numFields; f++) {
2537:             offsets[f + 1]       += offsets[f];
2538:             offsetsCopy[f + 1]    = offsets[f + 1];
2539:             newOffsets[f + 1]    += newOffsets[f];
2540:             newOffsetsCopy[f + 1] = newOffsets[f + 1];
2541:           }
2542:           DMPlexGetIndicesPointFields_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1,rowIndices);
2543:           for (a = 0; a < aDof; a++) {
2544:             PetscInt anchor = anchors[a + aOff], lOff;
2545:             PetscSectionGetOffset(localCoarse,anchor,&lOff);
2546:             DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1,colIndices);
2547:           }
2548:         }
2549:         else {
2550:           DMPlexGetIndicesPoint_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,rowIndices);
2551:           for (a = 0; a < aDof; a++) {
2552:             PetscInt anchor = anchors[a + aOff], lOff;
2553:             PetscSectionGetOffset(localCoarse,anchor,&lOff);
2554:             DMPlexGetIndicesPoint_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,colIndices);
2555:           }
2556:         }
2557:         if (numFields) {
2558:           PetscInt count, a;

2560:           for (f = 0, count = 0; f < numFields; f++) {
2561:             PetscInt iSize = offsets[f + 1] - offsets[f];
2562:             PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2563:             MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);
2564:             count += iSize * jSize;
2565:             pInd[numColIndices + f]             = offsets[f+1];
2566:             pInd[numColIndices + numFields + f] = newOffsets[f+1];
2567:           }
2568:           for (a = 0; a < aDof; a++) {
2569:             PetscInt anchor = anchors[a + aOff];
2570:             PetscInt gOff;
2571:             PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2572:             DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1,pInd);
2573:           }
2574:         }
2575:         else {
2576:           PetscInt a;
2577:           MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);
2578:           for (a = 0; a < aDof; a++) {
2579:             PetscInt anchor = anchors[a + aOff];
2580:             PetscInt gOff;
2581:             PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2582:             DMPlexGetIndicesPoint_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,pInd);
2583:           }
2584:         }
2585:         DMRestoreWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);
2586:         DMRestoreWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);
2587:       }
2588:       else {
2589:         PetscInt gOff;

2591:         PetscSectionGetOffset(globalCoarse,p,&gOff);
2592:         if (numFields) {
2593:           for (f = 0; f < numFields; f++) {
2594:             PetscInt fDof;
2595:             PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2596:             offsets[f + 1] = fDof + offsets[f];
2597:           }
2598:           for (f = 0; f < numFields; f++) {
2599:             pInd[numColIndices + f]             = offsets[f+1];
2600:             pInd[numColIndices + numFields + f] = offsets[f+1];
2601:           }
2602:           DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);
2603:         } else {
2604:           DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);
2605:         }
2606:       }
2607:     }
2608:     PetscFree(maxChildIds);
2609:   }
2610:   {
2611:     PetscSF  indicesSF, matricesSF;
2612:     PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;

2614:     PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
2615:     PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);
2616:     PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);
2617:     PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);
2618:     PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);
2619:     PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);
2620:     PetscSFDestroy(&coarseToFineEmbedded);
2621:     PetscFree(remoteOffsetsIndices);
2622:     PetscFree(remoteOffsetsMatrices);
2623:     PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);
2624:     PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);
2625:     PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);
2626:     PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices);
2627:     PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);
2628:     PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices);
2629:     PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);
2630:     PetscSFDestroy(&matricesSF);
2631:     PetscSFDestroy(&indicesSF);
2632:     PetscFree2(rootIndices,rootMatrices);
2633:     PetscSectionDestroy(&rootIndicesSec);
2634:     PetscSectionDestroy(&rootMatricesSec);
2635:   }
2636:   /* count to preallocate */
2637:   DMGetSection(fine,&localFine);
2638:   {
2639:     PetscInt    nGlobal;
2640:     PetscInt    *dnnz, *onnz;
2641:     PetscLayout rowMap, colMap;
2642:     PetscInt    rowStart, rowEnd, colStart, colEnd;
2643:     PetscInt    maxDof;
2644:     PetscInt    *rowIndices;
2645:     DM           refTree;
2646:     PetscInt     **refPointFieldN;
2647:     PetscScalar  ***refPointFieldMats;
2648:     PetscSection refConSec, refAnSec;
2649:     PetscInt     pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd;
2650:     PetscScalar  *pointWork;

2652:     PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);
2653:     PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);
2654:     MatGetLayouts(mat,&rowMap,&colMap);
2655:     PetscLayoutSetUp(rowMap);
2656:     PetscLayoutSetUp(colMap);
2657:     PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
2658:     PetscLayoutGetRange(colMap,&colStart,&colEnd);
2659:     PetscSectionGetMaxDof(globalFine,&maxDof);
2660:     PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd);
2661:     DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
2662:     for (p = leafStart; p < leafEnd; p++) {
2663:       PetscInt    gDof, gcDof, gOff;
2664:       PetscInt    numColIndices, pIndOff, *pInd;
2665:       PetscInt    matSize;
2666:       PetscInt    i;

2668:       PetscSectionGetDof(globalFine,p,&gDof);
2669:       PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2670:       if ((gDof - gcDof) <= 0) {
2671:         continue;
2672:       }
2673:       PetscSectionGetOffset(globalFine,p,&gOff);
2674:       if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset");
2675:       if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs");
2676:       PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2677:       PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2678:       numColIndices -= 2 * numFields;
2679:       if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from");
2680:       pInd = &leafIndices[pIndOff];
2681:       offsets[0]        = 0;
2682:       offsetsCopy[0]    = 0;
2683:       newOffsets[0]     = 0;
2684:       newOffsetsCopy[0] = 0;
2685:       if (numFields) {
2686:         PetscInt f;
2687:         for (f = 0; f < numFields; f++) {
2688:           PetscInt rowDof;

2690:           PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2691:           offsets[f + 1]        = offsets[f] + rowDof;
2692:           offsetsCopy[f + 1]    = offsets[f + 1];
2693:           newOffsets[f + 1]     = pInd[numColIndices + numFields + f];
2694:           numD[f] = 0;
2695:           numO[f] = 0;
2696:         }
2697:         DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);
2698:         for (f = 0; f < numFields; f++) {
2699:           PetscInt colOffset    = newOffsets[f];
2700:           PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];

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

2705:             if (gInd >= colStart && gInd < colEnd) {
2706:               numD[f]++;
2707:             }
2708:             else if (gInd >= 0) { /* negative means non-entry */
2709:               numO[f]++;
2710:             }
2711:           }
2712:         }
2713:       }
2714:       else {
2715:         DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);
2716:         numD[0] = 0;
2717:         numO[0] = 0;
2718:         for (i = 0; i < numColIndices; i++) {
2719:           PetscInt gInd = pInd[i];

2721:           if (gInd >= colStart && gInd < colEnd) {
2722:             numD[0]++;
2723:           }
2724:           else if (gInd >= 0) { /* negative means non-entry */
2725:             numO[0]++;
2726:           }
2727:         }
2728:       }
2729:       PetscSectionGetDof(leafMatricesSec,p,&matSize);
2730:       if (!matSize) { /* incoming matrix is identity */
2731:         PetscInt childId;

2733:         childId = childIds[p-pStartF];
2734:         if (childId < 0) { /* no child interpolation: one nnz per */
2735:           if (numFields) {
2736:             PetscInt f;
2737:             for (f = 0; f < numFields; f++) {
2738:               PetscInt numRows = offsets[f+1] - offsets[f], row;
2739:               for (row = 0; row < numRows; row++) {
2740:                 PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2741:                 PetscInt gIndFine   = rowIndices[offsets[f] + row];
2742:                 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2743:                   if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2744:                   dnnz[gIndFine - rowStart] = 1;
2745:                 }
2746:                 else if (gIndCoarse >= 0) { /* remote */
2747:                   if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2748:                   onnz[gIndFine - rowStart] = 1;
2749:                 }
2750:                 else { /* constrained */
2751:                   if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2752:                 }
2753:               }
2754:             }
2755:           }
2756:           else {
2757:             PetscInt i;
2758:             for (i = 0; i < gDof; i++) {
2759:               PetscInt gIndCoarse = pInd[i];
2760:               PetscInt gIndFine   = rowIndices[i];
2761:               if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2762:                 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2763:                 dnnz[gIndFine - rowStart] = 1;
2764:               }
2765:               else if (gIndCoarse >= 0) { /* remote */
2766:                 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2767:                 onnz[gIndFine - rowStart] = 1;
2768:               }
2769:               else { /* constrained */
2770:                 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2771:               }
2772:             }
2773:           }
2774:         }
2775:         else { /* interpolate from all */
2776:           if (numFields) {
2777:             PetscInt f;
2778:             for (f = 0; f < numFields; f++) {
2779:               PetscInt numRows = offsets[f+1] - offsets[f], row;
2780:               for (row = 0; row < numRows; row++) {
2781:                 PetscInt gIndFine = rowIndices[offsets[f] + row];
2782:                 if (gIndFine >= 0) {
2783:                   if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2784:                   dnnz[gIndFine - rowStart] = numD[f];
2785:                   onnz[gIndFine - rowStart] = numO[f];
2786:                 }
2787:               }
2788:             }
2789:           }
2790:           else {
2791:             PetscInt i;
2792:             for (i = 0; i < gDof; i++) {
2793:               PetscInt gIndFine = rowIndices[i];
2794:               if (gIndFine >= 0) {
2795:                 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2796:                 dnnz[gIndFine - rowStart] = numD[0];
2797:                 onnz[gIndFine - rowStart] = numO[0];
2798:               }
2799:             }
2800:           }
2801:         }
2802:       }
2803:       else { /* interpolate from all */
2804:         if (numFields) {
2805:           PetscInt f;
2806:           for (f = 0; f < numFields; f++) {
2807:             PetscInt numRows = offsets[f+1] - offsets[f], row;
2808:             for (row = 0; row < numRows; row++) {
2809:               PetscInt gIndFine = rowIndices[offsets[f] + row];
2810:               if (gIndFine >= 0) {
2811:                 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2812:                 dnnz[gIndFine - rowStart] = numD[f];
2813:                 onnz[gIndFine - rowStart] = numO[f];
2814:               }
2815:             }
2816:           }
2817:         }
2818:         else { /* every dof get a full row */
2819:           PetscInt i;
2820:           for (i = 0; i < gDof; i++) {
2821:             PetscInt gIndFine = rowIndices[i];
2822:             if (gIndFine >= 0) {
2823:               if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2824:               dnnz[gIndFine - rowStart] = numD[0];
2825:               onnz[gIndFine - rowStart] = numO[0];
2826:             }
2827:           }
2828:         }
2829:       }
2830:     }
2831:     MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);
2832:     PetscFree2(dnnz,onnz);

2834:     DMPlexGetReferenceTree(fine,&refTree);
2835:     DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
2836:     DMGetDefaultConstraints(refTree,&refConSec,NULL);
2837:     DMPlexGetAnchors(refTree,&refAnSec,NULL);
2838:     PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
2839:     PetscSectionGetMaxDof(refConSec,&maxConDof);
2840:     PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);
2841:     PetscMalloc1(maxConDof*maxColumns,&pointWork);
2842:     for (p = leafStart; p < leafEnd; p++) {
2843:       PetscInt gDof, gcDof, gOff;
2844:       PetscInt numColIndices, pIndOff, *pInd;
2845:       PetscInt matSize;
2846:       PetscInt childId;

2848:       PetscSectionGetDof(globalFine,p,&gDof);
2849:       PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2850:       if ((gDof - gcDof) <= 0) {
2851:         continue;
2852:       }
2853:       childId = childIds[p-pStartF];
2854:       PetscSectionGetOffset(globalFine,p,&gOff);
2855:       PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2856:       PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2857:       numColIndices -= 2 * numFields;
2858:       pInd = &leafIndices[pIndOff];
2859:       offsets[0]        = 0;
2860:       offsetsCopy[0]    = 0;
2861:       newOffsets[0]     = 0;
2862:       newOffsetsCopy[0] = 0;
2863:       rowOffsets[0]     = 0;
2864:       if (numFields) {
2865:         PetscInt f;
2866:         for (f = 0; f < numFields; f++) {
2867:           PetscInt rowDof;

2869:           PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2870:           offsets[f + 1]     = offsets[f] + rowDof;
2871:           offsetsCopy[f + 1] = offsets[f + 1];
2872:           rowOffsets[f + 1]  = pInd[numColIndices + f];
2873:           newOffsets[f + 1]  = pInd[numColIndices + numFields + f];
2874:         }
2875:         DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);
2876:       }
2877:       else {
2878:         DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);
2879:       }
2880:       PetscSectionGetDof(leafMatricesSec,p,&matSize);
2881:       if (!matSize) { /* incoming matrix is identity */
2882:         if (childId < 0) { /* no child interpolation: scatter */
2883:           if (numFields) {
2884:             PetscInt f;
2885:             for (f = 0; f < numFields; f++) {
2886:               PetscInt numRows = offsets[f+1] - offsets[f], row;
2887:               for (row = 0; row < numRows; row++) {
2888:                 MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);
2889:               }
2890:             }
2891:           }
2892:           else {
2893:             PetscInt numRows = gDof, row;
2894:             for (row = 0; row < numRows; row++) {
2895:               MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);
2896:             }
2897:           }
2898:         }
2899:         else { /* interpolate from all */
2900:           if (numFields) {
2901:             PetscInt f;
2902:             for (f = 0; f < numFields; f++) {
2903:               PetscInt numRows = offsets[f+1] - offsets[f];
2904:               PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2905:               MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);
2906:             }
2907:           }
2908:           else {
2909:             MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);
2910:           }
2911:         }
2912:       }
2913:       else { /* interpolate from all */
2914:         PetscInt    pMatOff;
2915:         PetscScalar *pMat;

2917:         PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);
2918:         pMat = &leafMatrices[pMatOff];
2919:         if (childId < 0) { /* copy the incoming matrix */
2920:           if (numFields) {
2921:             PetscInt f, count;
2922:             for (f = 0, count = 0; f < numFields; f++) {
2923:               PetscInt numRows = offsets[f+1]-offsets[f];
2924:               PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2925:               PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2926:               PetscScalar *inMat = &pMat[count];

2928:               MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);
2929:               count += numCols * numInRows;
2930:             }
2931:           }
2932:           else {
2933:             MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);
2934:           }
2935:         }
2936:         else { /* multiply the incoming matrix by the child interpolation */
2937:           if (numFields) {
2938:             PetscInt f, count;
2939:             for (f = 0, count = 0; f < numFields; f++) {
2940:               PetscInt numRows = offsets[f+1]-offsets[f];
2941:               PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2942:               PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2943:               PetscScalar *inMat = &pMat[count];
2944:               PetscInt i, j, k;
2945:               if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch");
2946:               for (i = 0; i < numRows; i++) {
2947:                 for (j = 0; j < numCols; j++) {
2948:                   PetscScalar val = 0.;
2949:                   for (k = 0; k < numInRows; k++) {
2950:                     val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2951:                   }
2952:                   pointWork[i * numCols + j] = val;
2953:                 }
2954:               }
2955:               MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);
2956:               count += numCols * numInRows;
2957:             }
2958:           }
2959:           else { /* every dof gets a full row */
2960:             PetscInt numRows   = gDof;
2961:             PetscInt numCols   = numColIndices;
2962:             PetscInt numInRows = matSize / numColIndices;
2963:             PetscInt i, j, k;
2964:             if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch");
2965:             for (i = 0; i < numRows; i++) {
2966:               for (j = 0; j < numCols; j++) {
2967:                 PetscScalar val = 0.;
2968:                 for (k = 0; k < numInRows; k++) {
2969:                   val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2970:                 }
2971:                 pointWork[i * numCols + j] = val;
2972:               }
2973:             }
2974:             MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);
2975:           }
2976:         }
2977:       }
2978:     }
2979:     DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
2980:     DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
2981:     PetscFree(pointWork);
2982:   }
2983:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
2984:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
2985:   PetscSectionDestroy(&leafIndicesSec);
2986:   PetscSectionDestroy(&leafMatricesSec);
2987:   PetscFree2(leafIndices,leafMatrices);
2988:   PetscFree2(*(PetscInt****)&perms,*(PetscScalar****)&flips);
2989:   PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);
2990:   ISRestoreIndices(aIS,&anchors);
2991:   return(0);
2992: }

2994: /*
2995:  * Assuming a nodal basis (w.r.t. the dual basis) basis:
2996:  *
2997:  * for each coarse dof \phi^c_i:
2998:  *   for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
2999:  *     for each fine dof \phi^f_j;
3000:  *       a_{i,j} = 0;
3001:  *       for each fine dof \phi^f_k:
3002:  *         a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
3003:  *                    [^^^ this is = \phi^c_i ^^^]
3004:  */
3005: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
3006: {
3007:   PetscDS        ds;
3008:   PetscSection   section, cSection;
3009:   DMLabel        canonical, depth;
3010:   Mat            cMat, mat;
3011:   PetscInt       *nnz;
3012:   PetscInt       f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
3013:   PetscInt       m, n;
3014:   PetscScalar    *pointScalar;
3015:   PetscReal      *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;

3019:   DMGetSection(refTree,&section);
3020:   DMGetDimension(refTree, &dim);
3021:   PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);
3022:   PetscMalloc2(dim,&pointScalar,dim,&pointRef);
3023:   DMGetDS(refTree,&ds);
3024:   PetscDSGetNumFields(ds,&numFields);
3025:   PetscSectionGetNumFields(section,&numSecFields);
3026:   DMGetLabel(refTree,"canonical",&canonical);
3027:   DMGetLabel(refTree,"depth",&depth);
3028:   DMGetDefaultConstraints(refTree,&cSection,&cMat);
3029:   DMPlexGetChart(refTree, &pStart, &pEnd);
3030:   DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);
3031:   MatGetSize(cMat,&n,&m); /* the injector has transpose sizes from the constraint matrix */
3032:   /* Step 1: compute non-zero pattern.  A proper subset of constraint matrix non-zero */
3033:   PetscCalloc1(m,&nnz);
3034:   for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
3035:     const PetscInt *children;
3036:     PetscInt numChildren;
3037:     PetscInt i, numChildDof, numSelfDof;

3039:     if (canonical) {
3040:       PetscInt pCanonical;
3041:       DMLabelGetValue(canonical,p,&pCanonical);
3042:       if (p != pCanonical) continue;
3043:     }
3044:     DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3045:     if (!numChildren) continue;
3046:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
3047:       PetscInt child = children[i];
3048:       PetscInt dof;

3050:       PetscSectionGetDof(section,child,&dof);
3051:       numChildDof += dof;
3052:     }
3053:     PetscSectionGetDof(section,p,&numSelfDof);
3054:     if (!numChildDof || !numSelfDof) continue;
3055:     for (f = 0; f < numFields; f++) {
3056:       PetscInt selfOff;

3058:       if (numSecFields) { /* count the dofs for just this field */
3059:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
3060:           PetscInt child = children[i];
3061:           PetscInt dof;

3063:           PetscSectionGetFieldDof(section,child,f,&dof);
3064:           numChildDof += dof;
3065:         }
3066:         PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3067:         PetscSectionGetFieldOffset(section,p,f,&selfOff);
3068:       }
3069:       else {
3070:         PetscSectionGetOffset(section,p,&selfOff);
3071:       }
3072:       for (i = 0; i < numSelfDof; i++) {
3073:         nnz[selfOff + i] = numChildDof;
3074:       }
3075:     }
3076:   }
3077:   MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);
3078:   PetscFree(nnz);
3079:   /* Setp 2: compute entries */
3080:   for (p = pStart; p < pEnd; p++) {
3081:     const PetscInt *children;
3082:     PetscInt numChildren;
3083:     PetscInt i, numChildDof, numSelfDof;

3085:     /* same conditions about when entries occur */
3086:     if (canonical) {
3087:       PetscInt pCanonical;
3088:       DMLabelGetValue(canonical,p,&pCanonical);
3089:       if (p != pCanonical) continue;
3090:     }
3091:     DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3092:     if (!numChildren) continue;
3093:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
3094:       PetscInt child = children[i];
3095:       PetscInt dof;

3097:       PetscSectionGetDof(section,child,&dof);
3098:       numChildDof += dof;
3099:     }
3100:     PetscSectionGetDof(section,p,&numSelfDof);
3101:     if (!numChildDof || !numSelfDof) continue;

3103:     for (f = 0; f < numFields; f++) {
3104:       PetscInt       pI = -1, cI = -1;
3105:       PetscInt       selfOff, Nc, parentCell;
3106:       PetscInt       cellShapeOff;
3107:       PetscObject    disc;
3108:       PetscDualSpace dsp;
3109:       PetscClassId   classId;
3110:       PetscScalar    *pointMat;
3111:       PetscInt       *matRows, *matCols;
3112:       PetscInt       pO = PETSC_MIN_INT;
3113:       const PetscInt *depthNumDof;

3115:       if (numSecFields) {
3116:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
3117:           PetscInt child = children[i];
3118:           PetscInt dof;

3120:           PetscSectionGetFieldDof(section,child,f,&dof);
3121:           numChildDof += dof;
3122:         }
3123:         PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3124:         PetscSectionGetFieldOffset(section,p,f,&selfOff);
3125:       }
3126:       else {
3127:         PetscSectionGetOffset(section,p,&selfOff);
3128:       }

3130:       /* find a cell whose closure contains p */
3131:       if (p >= cStart && p < cEnd) {
3132:         parentCell = p;
3133:       }
3134:       else {
3135:         PetscInt *star = NULL;
3136:         PetscInt numStar;

3138:         parentCell = -1;
3139:         DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3140:         for (i = numStar - 1; i >= 0; i--) {
3141:           PetscInt c = star[2 * i];

3143:           if (c >= cStart && c < cEnd) {
3144:             parentCell = c;
3145:             break;
3146:           }
3147:         }
3148:         DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3149:       }
3150:       /* determine the offset of p's shape functions withing parentCell's shape functions */
3151:       PetscDSGetDiscretization(ds,f,&disc);
3152:       PetscObjectGetClassId(disc,&classId);
3153:       if (classId == PETSCFE_CLASSID) {
3154:         PetscFEGetDualSpace((PetscFE)disc,&dsp);
3155:       }
3156:       else if (classId == PETSCFV_CLASSID) {
3157:         PetscFVGetDualSpace((PetscFV)disc,&dsp);
3158:       }
3159:       else {
3160:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object");
3161:       }
3162:       PetscDualSpaceGetNumDof(dsp,&depthNumDof);
3163:       PetscDualSpaceGetNumComponents(dsp,&Nc);
3164:       {
3165:         PetscInt *closure = NULL;
3166:         PetscInt numClosure;

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

3172:           pO = closure[2 * i + 1];
3173:           if (point == p) {
3174:             pI = i;
3175:             break;
3176:           }
3177:           DMLabelGetValue(depth,point,&pointDepth);
3178:           cellShapeOff += depthNumDof[pointDepth];
3179:         }
3180:         DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);
3181:       }

3183:       DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);
3184:       DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);
3185:       matCols = matRows + numSelfDof;
3186:       for (i = 0; i < numSelfDof; i++) {
3187:         matRows[i] = selfOff + i;
3188:       }
3189:       for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
3190:       {
3191:         PetscInt colOff = 0;

3193:         for (i = 0; i < numChildren; i++) {
3194:           PetscInt child = children[i];
3195:           PetscInt dof, off, j;

3197:           if (numSecFields) {
3198:             PetscSectionGetFieldDof(cSection,child,f,&dof);
3199:             PetscSectionGetFieldOffset(cSection,child,f,&off);
3200:           }
3201:           else {
3202:             PetscSectionGetDof(cSection,child,&dof);
3203:             PetscSectionGetOffset(cSection,child,&off);
3204:           }

3206:           for (j = 0; j < dof; j++) {
3207:             matCols[colOff++] = off + j;
3208:           }
3209:         }
3210:       }
3211:       if (classId == PETSCFE_CLASSID) {
3212:         PetscFE        fe = (PetscFE) disc;
3213:         PetscInt       fSize;
3214:         const PetscInt ***perms;
3215:         const PetscScalar ***flips;
3216:         const PetscInt *pperms;


3219:         PetscFEGetDualSpace(fe,&dsp);
3220:         PetscDualSpaceGetDimension(dsp,&fSize);
3221:         PetscDualSpaceGetSymmetries(dsp, &perms, &flips);
3222:         pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL;
3223:         for (i = 0; i < numSelfDof; i++) { /* for every shape function */
3224:           PetscQuadrature q;
3225:           PetscInt        dim, thisNc, numPoints, j, k;
3226:           const PetscReal *points;
3227:           const PetscReal *weights;
3228:           PetscInt        *closure = NULL;
3229:           PetscInt        numClosure;
3230:           PetscInt        iCell = pperms ? pperms[i] : i;
3231:           PetscInt        parentCellShapeDof = cellShapeOff + iCell;
3232:           PetscReal       *Bparent;

3234:           PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);
3235:           PetscQuadratureGetData(q,&dim,&thisNc,&numPoints,&points,&weights);
3236:           if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc);
3237:           PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3238:           for (j = 0; j < numPoints; j++) {
3239:             PetscInt          childCell = -1;
3240:             PetscReal         *parentValAtPoint;
3241:             const PetscReal   xi0[3] = {-1.,-1.,-1.};
3242:             const PetscReal   *pointReal = &points[dim * j];
3243:             const PetscScalar *point;
3244:             PetscReal         *Bchild;
3245:             PetscInt          childCellShapeOff, pointMatOff;
3246: #if defined(PETSC_USE_COMPLEX)
3247:             PetscInt          d;

3249:             for (d = 0; d < dim; d++) {
3250:               pointScalar[d] = points[dim * j + d];
3251:             }
3252:             point = pointScalar;
3253: #else
3254:             point = pointReal;
3255: #endif

3257:             parentValAtPoint = &Bparent[(fSize * j + parentCellShapeDof) * Nc];

3259:             for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3260:               PetscInt child = children[k];
3261:               PetscInt *star = NULL;
3262:               PetscInt numStar, s;

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

3268:                 if (c < cStart || c >= cEnd) continue;
3269:                 DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);
3270:                 if (childCell >= 0) break;
3271:               }
3272:               DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);
3273:               if (childCell >= 0) break;
3274:             }
3275:             if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point");
3276:             DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);
3277:             DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);
3278:             CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
3279:             CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);

3281:             PetscFEGetTabulation(fe,1,pointRef,&Bchild,NULL,NULL);
3282:             DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3283:             for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3284:               PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT;
3285:               PetscInt l;
3286:               const PetscInt *cperms;

3288:               DMLabelGetValue(depth,child,&childDepth);
3289:               childDof = depthNumDof[childDepth];
3290:               for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) {
3291:                 PetscInt point = closure[2 * l];
3292:                 PetscInt pointDepth;

3294:                 childO = closure[2 * l + 1];
3295:                 if (point == child) {
3296:                   cI = l;
3297:                   break;
3298:                 }
3299:                 DMLabelGetValue(depth,point,&pointDepth);
3300:                 childCellShapeOff += depthNumDof[pointDepth];
3301:               }
3302:               if (l == numClosure) {
3303:                 pointMatOff += childDof;
3304:                 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3305:               }
3306:               cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL;
3307:               for (l = 0; l < childDof; l++) {
3308:                 PetscInt    lCell = cperms ? cperms[l] : l;
3309:                 PetscInt    childCellDof = childCellShapeOff + lCell;
3310:                 PetscReal   *childValAtPoint;
3311:                 PetscReal   val = 0.;

3313:                 childValAtPoint = &Bchild[childCellDof * Nc];
3314:                 for (m = 0; m < Nc; m++) {
3315:                   val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];
3316:                 }

3318:                 pointMat[i * numChildDof + pointMatOff + l] += val;
3319:               }
3320:               pointMatOff += childDof;
3321:             }
3322:             DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3323:             PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);
3324:           }
3325:           PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);
3326:         }
3327:       }
3328:       else { /* just the volume-weighted averages of the children */
3329:         PetscReal parentVol;
3330:         PetscInt  childCell;

3332:         DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);
3333:         for (i = 0, childCell = 0; i < numChildren; i++) {
3334:           PetscInt  child = children[i], j;
3335:           PetscReal childVol;

3337:           if (child < cStart || child >= cEnd) continue;
3338:           DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);
3339:           for (j = 0; j < Nc; j++) {
3340:             pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3341:           }
3342:           childCell++;
3343:         }
3344:       }
3345:       /* Insert pointMat into mat */
3346:       MatSetValues(mat,numSelfDof,matRows,numChildDof,matCols,pointMat,INSERT_VALUES);
3347:       DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);
3348:       DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);
3349:     }
3350:   }
3351:   PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);
3352:   PetscFree2(pointScalar,pointRef);
3353:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3354:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3355:   *inj = mat;
3356:   return(0);
3357: }

3359: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3360: {
3361:   PetscDS        ds;
3362:   PetscInt       numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3363:   PetscScalar    ***refPointFieldMats;
3364:   PetscSection   refConSec, refSection;

3368:   DMGetDS(refTree,&ds);
3369:   PetscDSGetNumFields(ds,&numFields);
3370:   DMGetDefaultConstraints(refTree,&refConSec,NULL);
3371:   DMGetSection(refTree,&refSection);
3372:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3373:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
3374:   PetscSectionGetMaxDof(refConSec,&maxDof);
3375:   PetscMalloc1(maxDof,&rows);
3376:   PetscMalloc1(maxDof*maxDof,&cols);
3377:   for (p = pRefStart; p < pRefEnd; p++) {
3378:     PetscInt parent, pDof, parentDof;

3380:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
3381:     PetscSectionGetDof(refConSec,p,&pDof);
3382:     PetscSectionGetDof(refSection,parent,&parentDof);
3383:     if (!pDof || !parentDof || parent == p) continue;

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

3389:       if (numFields > 1) {
3390:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3391:         PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
3392:       }
3393:       else {
3394:         PetscSectionGetDof(refConSec,p,&cDof);
3395:         PetscSectionGetOffset(refConSec,p,&cOff);
3396:       }

3398:       for (r = 0; r < cDof; r++) {
3399:         rows[r] = cOff + r;
3400:       }
3401:       numCols = 0;
3402:       {
3403:         PetscInt aDof, aOff, j;

3405:         if (numFields > 1) {
3406:           PetscSectionGetFieldDof(refSection,parent,f,&aDof);
3407:           PetscSectionGetFieldOffset(refSection,parent,f,&aOff);
3408:         }
3409:         else {
3410:           PetscSectionGetDof(refSection,parent,&aDof);
3411:           PetscSectionGetOffset(refSection,parent,&aOff);
3412:         }

3414:         for (j = 0; j < aDof; j++) {
3415:           cols[numCols++] = aOff + j;
3416:         }
3417:       }
3418:       PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
3419:       /* transpose of constraint matrix */
3420:       MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);
3421:     }
3422:   }
3423:   *childrenMats = refPointFieldMats;
3424:   PetscFree(rows);
3425:   PetscFree(cols);
3426:   return(0);
3427: }

3429: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3430: {
3431:   PetscDS        ds;
3432:   PetscScalar    ***refPointFieldMats;
3433:   PetscInt       numFields, pRefStart, pRefEnd, p, f;
3434:   PetscSection   refConSec, refSection;

3438:   refPointFieldMats = *childrenMats;
3439:   *childrenMats = NULL;
3440:   DMGetDS(refTree,&ds);
3441:   DMGetSection(refTree,&refSection);
3442:   PetscDSGetNumFields(ds,&numFields);
3443:   DMGetDefaultConstraints(refTree,&refConSec,NULL);
3444:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3445:   for (p = pRefStart; p < pRefEnd; p++) {
3446:     PetscInt parent, pDof, parentDof;

3448:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
3449:     PetscSectionGetDof(refConSec,p,&pDof);
3450:     PetscSectionGetDof(refSection,parent,&parentDof);
3451:     if (!pDof || !parentDof || parent == p) continue;

3453:     for (f = 0; f < numFields; f++) {
3454:       PetscInt cDof;

3456:       if (numFields > 1) {
3457:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3458:       }
3459:       else {
3460:         PetscSectionGetDof(refConSec,p,&cDof);
3461:       }

3463:       PetscFree(refPointFieldMats[p - pRefStart][f]);
3464:     }
3465:     PetscFree(refPointFieldMats[p - pRefStart]);
3466:   }
3467:   PetscFree(refPointFieldMats);
3468:   return(0);
3469: }

3471: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef)
3472: {
3473:   Mat            cMatRef;
3474:   PetscObject    injRefObj;

3478:   DMGetDefaultConstraints(refTree,NULL,&cMatRef);
3479:   PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);
3480:   *injRef = (Mat) injRefObj;
3481:   if (!*injRef) {
3482:     DMPlexComputeInjectorReferenceTree(refTree,injRef);
3483:     PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);
3484:     /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3485:     PetscObjectDereference((PetscObject)*injRef);
3486:   }
3487:   return(0);
3488: }

3490: 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)
3491: {
3492:   PetscInt       pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3493:   PetscSection   globalCoarse, globalFine;
3494:   PetscSection   localCoarse, localFine, leafIndicesSec;
3495:   PetscSection   multiRootSec, rootIndicesSec;
3496:   PetscInt       *leafInds, *rootInds = NULL;
3497:   const PetscInt *rootDegrees;
3498:   PetscScalar    *leafVals = NULL, *rootVals = NULL;
3499:   PetscSF        coarseToFineEmbedded;

3503:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3504:   DMPlexGetChart(fine,&pStartF,&pEndF);
3505:   DMGetSection(fine,&localFine);
3506:   DMGetGlobalSection(fine,&globalFine);
3507:   PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
3508:   PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);
3509:   PetscSectionGetMaxDof(localFine,&maxDof);
3510:   { /* winnow fine points that don't have global dofs out of the sf */
3511:     PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3512:     const PetscInt *leaves;

3514:     PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);
3515:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3516:       p    = leaves ? leaves[l] : l;
3517:       PetscSectionGetDof(globalFine,p,&dof);
3518:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3519:       if ((dof - cdof) > 0) {
3520:         numPointsWithDofs++;

3522:         PetscSectionGetDof(localFine,p,&dof);
3523:         PetscSectionSetDof(leafIndicesSec,p,dof + 1);
3524:       }
3525:     }
3526:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3527:     PetscSectionSetUp(leafIndicesSec);
3528:     PetscSectionGetStorageSize(leafIndicesSec,&numIndices);
3529:     PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);
3530:     if (gatheredValues)  {PetscMalloc1(numIndices,&leafVals);}
3531:     for (l = 0, offset = 0; l < nleaves; l++) {
3532:       p    = leaves ? leaves[l] : l;
3533:       PetscSectionGetDof(globalFine,p,&dof);
3534:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3535:       if ((dof - cdof) > 0) {
3536:         PetscInt    off, gOff;
3537:         PetscInt    *pInd;
3538:         PetscScalar *pVal = NULL;

3540:         pointsWithDofs[offset++] = l;

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

3544:         pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3545:         if (gatheredValues) {
3546:           PetscInt i;

3548:           pVal = &leafVals[off + 1];
3549:           for (i = 0; i < dof; i++) pVal[i] = 0.;
3550:         }
3551:         PetscSectionGetOffset(globalFine,p,&gOff);

3553:         offsets[0] = 0;
3554:         if (numFields) {
3555:           PetscInt f;

3557:           for (f = 0; f < numFields; f++) {
3558:             PetscInt fDof;
3559:             PetscSectionGetFieldDof(localFine,p,f,&fDof);
3560:             offsets[f + 1] = fDof + offsets[f];
3561:           }
3562:           DMPlexGetIndicesPointFields_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);
3563:         } else {
3564:           DMPlexGetIndicesPoint_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);
3565:         }
3566:         if (gatheredValues) {VecGetValues(fineVec,dof,pInd,pVal);}
3567:       }
3568:     }
3569:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
3570:     PetscFree(pointsWithDofs);
3571:   }

3573:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3574:   DMGetSection(coarse,&localCoarse);
3575:   DMGetGlobalSection(coarse,&globalCoarse);

3577:   { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3578:     MPI_Datatype threeInt;
3579:     PetscMPIInt  rank;
3580:     PetscInt     (*parentNodeAndIdCoarse)[3];
3581:     PetscInt     (*parentNodeAndIdFine)[3];
3582:     PetscInt     p, nleaves, nleavesToParents;
3583:     PetscSF      pointSF, sfToParents;
3584:     const PetscInt *ilocal;
3585:     const PetscSFNode *iremote;
3586:     PetscSFNode  *iremoteToParents;
3587:     PetscInt     *ilocalToParents;

3589:     MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
3590:     MPI_Type_contiguous(3,MPIU_INT,&threeInt);
3591:     MPI_Type_commit(&threeInt);
3592:     PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);
3593:     DMGetPointSF(coarse,&pointSF);
3594:     PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);
3595:     for (p = pStartC; p < pEndC; p++) {
3596:       PetscInt parent, childId;
3597:       DMPlexGetTreeParent(coarse,p,&parent,&childId);
3598:       parentNodeAndIdCoarse[p - pStartC][0] = rank;
3599:       parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3600:       parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3601:       if (nleaves > 0) {
3602:         PetscInt leaf = -1;

3604:         if (ilocal) {
3605:           PetscFindInt(parent,nleaves,ilocal,&leaf);
3606:         }
3607:         else {
3608:           leaf = p - pStartC;
3609:         }
3610:         if (leaf >= 0) {
3611:           parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3612:           parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3613:         }
3614:       }
3615:     }
3616:     for (p = pStartF; p < pEndF; p++) {
3617:       parentNodeAndIdFine[p - pStartF][0] = -1;
3618:       parentNodeAndIdFine[p - pStartF][1] = -1;
3619:       parentNodeAndIdFine[p - pStartF][2] = -1;
3620:     }
3621:     PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3622:     PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3623:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3624:       PetscInt dof;

3626:       PetscSectionGetDof(leafIndicesSec,p,&dof);
3627:       if (dof) {
3628:         PetscInt off;

3630:         PetscSectionGetOffset(leafIndicesSec,p,&off);
3631:         if (gatheredIndices) {
3632:           leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3633:         } else if (gatheredValues) {
3634:           leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3635:         }
3636:       }
3637:       if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3638:         nleavesToParents++;
3639:       }
3640:     }
3641:     PetscMalloc1(nleavesToParents,&ilocalToParents);
3642:     PetscMalloc1(nleavesToParents,&iremoteToParents);
3643:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3644:       if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3645:         ilocalToParents[nleavesToParents] = p - pStartF;
3646:         iremoteToParents[nleavesToParents].rank  = parentNodeAndIdFine[p-pStartF][0];
3647:         iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1];
3648:         nleavesToParents++;
3649:       }
3650:     }
3651:     PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);
3652:     PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);
3653:     PetscSFDestroy(&coarseToFineEmbedded);

3655:     coarseToFineEmbedded = sfToParents;

3657:     PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);
3658:     MPI_Type_free(&threeInt);
3659:   }

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

3665:     for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3666:       PetscSectionGetDof(globalCoarse,p,&dof);
3667:       PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3668:       if ((dof - cdof) > 0) {
3669:         numPointsWithDofs++;
3670:       }
3671:     }
3672:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3673:     for (p = pStartC, offset = 0; p < pEndC; p++) {
3674:       PetscSectionGetDof(globalCoarse,p,&dof);
3675:       PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3676:       if ((dof - cdof) > 0) {
3677:         pointsWithDofs[offset++] = p - pStartC;
3678:       }
3679:     }
3680:     PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);
3681:     PetscSFDestroy(&coarseToFineEmbedded);
3682:     PetscFree(pointsWithDofs);
3683:     coarseToFineEmbedded = sfDofsOnly;
3684:   }

3686:   /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3687:   PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);
3688:   PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);
3689:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);
3690:   PetscSectionSetChart(multiRootSec,pStartC,pEndC);
3691:   for (p = pStartC; p < pEndC; p++) {
3692:     PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);
3693:   }
3694:   PetscSectionSetUp(multiRootSec);
3695:   PetscSectionGetStorageSize(multiRootSec,&numMulti);
3696:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
3697:   { /* distribute the leaf section */
3698:     PetscSF multi, multiInv, indicesSF;
3699:     PetscInt *remoteOffsets, numRootIndices;

3701:     PetscSFGetMultiSF(coarseToFineEmbedded,&multi);
3702:     PetscSFCreateInverseSF(multi,&multiInv);
3703:     PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);
3704:     PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);
3705:     PetscFree(remoteOffsets);
3706:     PetscSFDestroy(&multiInv);
3707:     PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
3708:     if (gatheredIndices) {
3709:       PetscMalloc1(numRootIndices,&rootInds);
3710:       PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds);
3711:       PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds);
3712:     }
3713:     if (gatheredValues) {
3714:       PetscMalloc1(numRootIndices,&rootVals);
3715:       PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals);
3716:       PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals);
3717:     }
3718:     PetscSFDestroy(&indicesSF);
3719:   }
3720:   PetscSectionDestroy(&leafIndicesSec);
3721:   PetscFree(leafInds);
3722:   PetscFree(leafVals);
3723:   PetscSFDestroy(&coarseToFineEmbedded);
3724:   *rootMultiSec = multiRootSec;
3725:   *multiLeafSec = rootIndicesSec;
3726:   if (gatheredIndices) *gatheredIndices = rootInds;
3727:   if (gatheredValues)  *gatheredValues  = rootVals;
3728:   return(0);
3729: }

3731: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3732: {
3733:   DM             refTree;
3734:   PetscSection   multiRootSec, rootIndicesSec;
3735:   PetscSection   globalCoarse, globalFine;
3736:   PetscSection   localCoarse, localFine;
3737:   PetscSection   cSecRef;
3738:   PetscInt       *rootIndices, *parentIndices, pRefStart, pRefEnd;
3739:   Mat            injRef;
3740:   PetscInt       numFields, maxDof;
3741:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
3742:   PetscInt       *offsets, *offsetsCopy, *rowOffsets;
3743:   PetscLayout    rowMap, colMap;
3744:   PetscInt       rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3745:   PetscScalar    ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */


3750:   /* get the templates for the fine-to-coarse injection from the reference tree */
3751:   DMPlexGetReferenceTree(coarse,&refTree);
3752:   DMGetDefaultConstraints(refTree,&cSecRef,NULL);
3753:   PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
3754:   DMPlexReferenceTreeGetInjector(refTree,&injRef);

3756:   DMPlexGetChart(fine,&pStartF,&pEndF);
3757:   DMGetSection(fine,&localFine);
3758:   DMGetGlobalSection(fine,&globalFine);
3759:   PetscSectionGetNumFields(localFine,&numFields);
3760:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3761:   DMGetSection(coarse,&localCoarse);
3762:   DMGetGlobalSection(coarse,&globalCoarse);
3763:   PetscSectionGetMaxDof(localCoarse,&maxDof);
3764:   {
3765:     PetscInt maxFields = PetscMax(1,numFields) + 1;
3766:     PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
3767:   }

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

3771:   PetscMalloc1(maxDof,&parentIndices);

3773:   /* count indices */
3774:   MatGetLayouts(mat,&rowMap,&colMap);
3775:   PetscLayoutSetUp(rowMap);
3776:   PetscLayoutSetUp(colMap);
3777:   PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
3778:   PetscLayoutGetRange(colMap,&colStart,&colEnd);
3779:   PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);
3780:   for (p = pStartC; p < pEndC; p++) {
3781:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3783:     PetscSectionGetDof(globalCoarse,p,&dof);
3784:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3785:     if ((dof - cdof) <= 0) continue;
3786:     PetscSectionGetOffset(globalCoarse,p,&gOff);

3788:     rowOffsets[0] = 0;
3789:     offsetsCopy[0] = 0;
3790:     if (numFields) {
3791:       PetscInt f;

3793:       for (f = 0; f < numFields; f++) {
3794:         PetscInt fDof;
3795:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3796:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3797:       }
3798:       DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);
3799:     } else {
3800:       DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);
3801:       rowOffsets[1] = offsetsCopy[0];
3802:     }

3804:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
3805:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
3806:     leafEnd = leafStart + numLeaves;
3807:     for (l = leafStart; l < leafEnd; l++) {
3808:       PetscInt numIndices, childId, offset;
3809:       const PetscInt *childIndices;

3811:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3812:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
3813:       childId = rootIndices[offset++];
3814:       childIndices = &rootIndices[offset];
3815:       numIndices--;

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

3820:         for (i = 0; i < numIndices; i++) {
3821:           PetscInt colIndex = childIndices[i];
3822:           PetscInt rowIndex = parentIndices[i];
3823:           if (rowIndex < 0) continue;
3824:           if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse");
3825:           if (colIndex >= colStart && colIndex < colEnd) {
3826:             nnzD[rowIndex - rowStart] = 1;
3827:           }
3828:           else {
3829:             nnzO[rowIndex - rowStart] = 1;
3830:           }
3831:         }
3832:       }
3833:       else {
3834:         PetscInt parentId, f, lim;

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

3838:         lim = PetscMax(1,numFields);
3839:         offsets[0] = 0;
3840:         if (numFields) {
3841:           PetscInt f;

3843:           for (f = 0; f < numFields; f++) {
3844:             PetscInt fDof;
3845:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

3847:             offsets[f + 1] = fDof + offsets[f];
3848:           }
3849:         }
3850:         else {
3851:           PetscInt cDof;

3853:           PetscSectionGetDof(cSecRef,childId,&cDof);
3854:           offsets[1] = cDof;
3855:         }
3856:         for (f = 0; f < lim; f++) {
3857:           PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3858:           PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3859:           PetscInt i, numD = 0, numO = 0;

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

3864:             if (colIndex < 0) continue;
3865:             if (colIndex >= colStart && colIndex < colEnd) {
3866:               numD++;
3867:             }
3868:             else {
3869:               numO++;
3870:             }
3871:           }
3872:           for (i = parentStart; i < parentEnd; i++) {
3873:             PetscInt rowIndex = parentIndices[i];

3875:             if (rowIndex < 0) continue;
3876:             nnzD[rowIndex - rowStart] += numD;
3877:             nnzO[rowIndex - rowStart] += numO;
3878:           }
3879:         }
3880:       }
3881:     }
3882:   }
3883:   /* preallocate */
3884:   MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);
3885:   PetscFree2(nnzD,nnzO);
3886:   /* insert values */
3887:   DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3888:   for (p = pStartC; p < pEndC; p++) {
3889:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3891:     PetscSectionGetDof(globalCoarse,p,&dof);
3892:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3893:     if ((dof - cdof) <= 0) continue;
3894:     PetscSectionGetOffset(globalCoarse,p,&gOff);

3896:     rowOffsets[0] = 0;
3897:     offsetsCopy[0] = 0;
3898:     if (numFields) {
3899:       PetscInt f;

3901:       for (f = 0; f < numFields; f++) {
3902:         PetscInt fDof;
3903:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3904:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3905:       }
3906:       DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);
3907:     } else {
3908:       DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);
3909:       rowOffsets[1] = offsetsCopy[0];
3910:     }

3912:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
3913:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
3914:     leafEnd = leafStart + numLeaves;
3915:     for (l = leafStart; l < leafEnd; l++) {
3916:       PetscInt numIndices, childId, offset;
3917:       const PetscInt *childIndices;

3919:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3920:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
3921:       childId = rootIndices[offset++];
3922:       childIndices = &rootIndices[offset];
3923:       numIndices--;

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

3928:         for (i = 0; i < numIndices; i++) {
3929:           MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);
3930:         }
3931:       }
3932:       else {
3933:         PetscInt parentId, f, lim;

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

3937:         lim = PetscMax(1,numFields);
3938:         offsets[0] = 0;
3939:         if (numFields) {
3940:           PetscInt f;

3942:           for (f = 0; f < numFields; f++) {
3943:             PetscInt fDof;
3944:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

3946:             offsets[f + 1] = fDof + offsets[f];
3947:           }
3948:         }
3949:         else {
3950:           PetscInt cDof;

3952:           PetscSectionGetDof(cSecRef,childId,&cDof);
3953:           offsets[1] = cDof;
3954:         }
3955:         for (f = 0; f < lim; f++) {
3956:           PetscScalar    *childMat   = &childrenMats[childId - pRefStart][f][0];
3957:           PetscInt       *rowIndices = &parentIndices[rowOffsets[f]];
3958:           const PetscInt *colIndices = &childIndices[offsets[f]];

3960:           MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);
3961:         }
3962:       }
3963:     }
3964:   }
3965:   PetscSectionDestroy(&multiRootSec);
3966:   PetscSectionDestroy(&rootIndicesSec);
3967:   PetscFree(parentIndices);
3968:   DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3969:   PetscFree(rootIndices);
3970:   PetscFree3(offsets,offsetsCopy,rowOffsets);

3972:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3973:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3974:   return(0);
3975: }

3977: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3978: {
3979:   PetscSF           coarseToFineEmbedded;
3980:   PetscSection      globalCoarse, globalFine;
3981:   PetscSection      localCoarse, localFine;
3982:   PetscSection      aSec, cSec;
3983:   PetscSection      rootValuesSec;
3984:   PetscSection      leafValuesSec;
3985:   PetscScalar       *rootValues, *leafValues;
3986:   IS                aIS;
3987:   const PetscInt    *anchors;
3988:   Mat               cMat;
3989:   PetscInt          numFields;
3990:   PetscInt          pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd, cellEndInterior;
3991:   PetscInt          aStart, aEnd, cStart, cEnd;
3992:   PetscInt          *maxChildIds;
3993:   PetscInt          *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3994:   PetscFV           fv = NULL;
3995:   PetscInt          dim, numFVcomps = -1, fvField = -1;
3996:   DM                cellDM = NULL, gradDM = NULL;
3997:   const PetscScalar *cellGeomArray = NULL;
3998:   const PetscScalar *gradArray = NULL;
3999:   PetscErrorCode    ierr;

4002:   VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
4003:   DMPlexGetChart(coarse,&pStartC,&pEndC);
4004:   DMPlexGetHeightStratum(coarse,0,&cellStart,&cellEnd);
4005:   DMPlexGetHybridBounds(coarse,&cellEndInterior,NULL,NULL,NULL);
4006:   cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior;
4007:   DMPlexGetChart(fine,&pStartF,&pEndF);
4008:   DMGetGlobalSection(fine,&globalFine);
4009:   DMGetCoordinateDim(coarse,&dim);
4010:   { /* winnow fine points that don't have global dofs out of the sf */
4011:     PetscInt       nleaves, l;
4012:     const PetscInt *leaves;
4013:     PetscInt       dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;

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

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

4020:       PetscSectionGetDof(globalFine,p,&dof);
4021:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
4022:       if ((dof - cdof) > 0) {
4023:         numPointsWithDofs++;
4024:       }
4025:     }
4026:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
4027:     for (l = 0, offset = 0; l < nleaves; l++) {
4028:       PetscInt p = leaves ? leaves[l] : l;

4030:       PetscSectionGetDof(globalFine,p,&dof);
4031:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
4032:       if ((dof - cdof) > 0) {
4033:         pointsWithDofs[offset++] = l;
4034:       }
4035:     }
4036:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
4037:     PetscFree(pointsWithDofs);
4038:   }
4039:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
4040:   PetscMalloc1(pEndC-pStartC,&maxChildIds);
4041:   for (p = pStartC; p < pEndC; p++) {
4042:     maxChildIds[p - pStartC] = -2;
4043:   }
4044:   PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);
4045:   PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);

4047:   DMGetSection(coarse,&localCoarse);
4048:   DMGetGlobalSection(coarse,&globalCoarse);

4050:   DMPlexGetAnchors(coarse,&aSec,&aIS);
4051:   ISGetIndices(aIS,&anchors);
4052:   PetscSectionGetChart(aSec,&aStart,&aEnd);

4054:   DMGetDefaultConstraints(coarse,&cSec,&cMat);
4055:   PetscSectionGetChart(cSec,&cStart,&cEnd);

4057:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
4058:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);
4059:   PetscSectionSetChart(rootValuesSec,pStartC,pEndC);
4060:   PetscSectionGetNumFields(localCoarse,&numFields);
4061:   {
4062:     PetscInt maxFields = PetscMax(1,numFields) + 1;
4063:     PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);
4064:   }
4065:   if (grad) {
4066:     PetscInt i;

4068:     VecGetDM(cellGeom,&cellDM);
4069:     VecGetArrayRead(cellGeom,&cellGeomArray);
4070:     VecGetDM(grad,&gradDM);
4071:     VecGetArrayRead(grad,&gradArray);
4072:     for (i = 0; i < PetscMax(1,numFields); i++) {
4073:       PetscObject  obj;
4074:       PetscClassId id;

4076:       DMGetField(coarse, i, NULL, &obj);
4077:       PetscObjectGetClassId(obj,&id);
4078:       if (id == PETSCFV_CLASSID) {
4079:         fv      = (PetscFV) obj;
4080:         PetscFVGetNumComponents(fv,&numFVcomps);
4081:         fvField = i;
4082:         break;
4083:       }
4084:     }
4085:   }

4087:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
4088:     PetscInt dof;
4089:     PetscInt maxChildId     = maxChildIds[p - pStartC];
4090:     PetscInt numValues      = 0;

4092:     PetscSectionGetDof(globalCoarse,p,&dof);
4093:     if (dof < 0) {
4094:       dof = -(dof + 1);
4095:     }
4096:     offsets[0]    = 0;
4097:     newOffsets[0] = 0;
4098:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
4099:       PetscInt *closure = NULL, closureSize, cl;

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

4105:         PetscSectionGetDof(localCoarse,c,&clDof);
4106:         numValues += clDof;
4107:       }
4108:       DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
4109:     }
4110:     else if (maxChildId == -1) {
4111:       PetscSectionGetDof(localCoarse,p,&numValues);
4112:     }
4113:     /* we will pack the column indices with the field offsets */
4114:     if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
4115:       /* also send the centroid, and the gradient */
4116:       numValues += dim * (1 + numFVcomps);
4117:     }
4118:     PetscSectionSetDof(rootValuesSec,p,numValues);
4119:   }
4120:   PetscSectionSetUp(rootValuesSec);
4121:   {
4122:     PetscInt          numRootValues;
4123:     const PetscScalar *coarseArray;

4125:     PetscSectionGetStorageSize(rootValuesSec,&numRootValues);
4126:     PetscMalloc1(numRootValues,&rootValues);
4127:     VecGetArrayRead(vecCoarseLocal,&coarseArray);
4128:     for (p = pStartC; p < pEndC; p++) {
4129:       PetscInt    numValues;
4130:       PetscInt    pValOff;
4131:       PetscScalar *pVal;
4132:       PetscInt    maxChildId = maxChildIds[p - pStartC];

4134:       PetscSectionGetDof(rootValuesSec,p,&numValues);
4135:       if (!numValues) {
4136:         continue;
4137:       }
4138:       PetscSectionGetOffset(rootValuesSec,p,&pValOff);
4139:       pVal = &(rootValues[pValOff]);
4140:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
4141:         PetscInt closureSize = numValues;
4142:         DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);
4143:         if (grad && p >= cellStart && p < cellEnd) {
4144:           PetscFVCellGeom *cg;
4145:           PetscScalar     *gradVals = NULL;
4146:           PetscInt        i;

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

4150:           DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);
4151:           for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
4152:           pVal += dim;
4153:           DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);
4154:           for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
4155:         }
4156:       }
4157:       else if (maxChildId == -1) {
4158:         PetscInt lDof, lOff, i;

4160:         PetscSectionGetDof(localCoarse,p,&lDof);
4161:         PetscSectionGetOffset(localCoarse,p,&lOff);
4162:         for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
4163:       }
4164:     }
4165:     VecRestoreArrayRead(vecCoarseLocal,&coarseArray);
4166:     PetscFree(maxChildIds);
4167:   }
4168:   {
4169:     PetscSF  valuesSF;
4170:     PetscInt *remoteOffsetsValues, numLeafValues;

4172:     PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);
4173:     PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);
4174:     PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);
4175:     PetscSFDestroy(&coarseToFineEmbedded);
4176:     PetscFree(remoteOffsetsValues);
4177:     PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);
4178:     PetscMalloc1(numLeafValues,&leafValues);
4179:     PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues);
4180:     PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues);
4181:     PetscSFDestroy(&valuesSF);
4182:     PetscFree(rootValues);
4183:     PetscSectionDestroy(&rootValuesSec);
4184:   }
4185:   DMGetSection(fine,&localFine);
4186:   {
4187:     PetscInt    maxDof;
4188:     PetscInt    *rowIndices;
4189:     DM           refTree;
4190:     PetscInt     **refPointFieldN;
4191:     PetscScalar  ***refPointFieldMats;
4192:     PetscSection refConSec, refAnSec;
4193:     PetscInt     pRefStart,pRefEnd,leafStart,leafEnd;
4194:     PetscScalar  *pointWork;

4196:     PetscSectionGetMaxDof(localFine,&maxDof);
4197:     DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
4198:     DMGetWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);
4199:     DMPlexGetReferenceTree(fine,&refTree);
4200:     DMCopyDisc(fine,refTree);
4201:     DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
4202:     DMGetDefaultConstraints(refTree,&refConSec,NULL);
4203:     DMPlexGetAnchors(refTree,&refAnSec,NULL);
4204:     PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
4205:     PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);
4206:     DMPlexGetHeightStratum(fine,0,&cellStart,&cellEnd);
4207:     DMPlexGetHybridBounds(fine,&cellEndInterior,NULL,NULL,NULL);
4208:     cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior;
4209:     for (p = leafStart; p < leafEnd; p++) {
4210:       PetscInt          gDof, gcDof, gOff, lDof;
4211:       PetscInt          numValues, pValOff;
4212:       PetscInt          childId;
4213:       const PetscScalar *pVal;
4214:       const PetscScalar *fvGradData = NULL;

4216:       PetscSectionGetDof(globalFine,p,&gDof);
4217:       PetscSectionGetDof(localFine,p,&lDof);
4218:       PetscSectionGetConstraintDof(globalFine,p,&gcDof);
4219:       if ((gDof - gcDof) <= 0) {
4220:         continue;
4221:       }
4222:       PetscSectionGetOffset(globalFine,p,&gOff);
4223:       PetscSectionGetDof(leafValuesSec,p,&numValues);
4224:       if (!numValues) continue;
4225:       PetscSectionGetOffset(leafValuesSec,p,&pValOff);
4226:       pVal = &leafValues[pValOff];
4227:       offsets[0]        = 0;
4228:       offsetsCopy[0]    = 0;
4229:       newOffsets[0]     = 0;
4230:       newOffsetsCopy[0] = 0;
4231:       childId           = cids[p - pStartF];
4232:       if (numFields) {
4233:         PetscInt f;
4234:         for (f = 0; f < numFields; f++) {
4235:           PetscInt rowDof;

4237:           PetscSectionGetFieldDof(localFine,p,f,&rowDof);
4238:           offsets[f + 1]        = offsets[f] + rowDof;
4239:           offsetsCopy[f + 1]    = offsets[f + 1];
4240:           /* TODO: closure indices */
4241:           newOffsets[f + 1]     = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
4242:         }
4243:         DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);
4244:       }
4245:       else {
4246:         offsets[0]    = 0;
4247:         offsets[1]    = lDof;
4248:         newOffsets[0] = 0;
4249:         newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
4250:         DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);
4251:       }
4252:       if (childId == -1) { /* no child interpolation: one nnz per */
4253:         VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);
4254:       } else {
4255:         PetscInt f;

4257:         if (grad && p >= cellStart && p < cellEnd) {
4258:           numValues -= (dim * (1 + numFVcomps));
4259:           fvGradData = &pVal[numValues];
4260:         }
4261:         for (f = 0; f < PetscMax(1,numFields); f++) {
4262:           const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
4263:           PetscInt numRows = offsets[f+1] - offsets[f];
4264:           PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
4265:           const PetscScalar *cVal = &pVal[newOffsets[f]];
4266:           PetscScalar *rVal = &pointWork[offsets[f]];
4267:           PetscInt i, j;

4269: #if 0
4270:           PetscInfo5(coarse,"childId %D, numRows %D, numCols %D, refPointFieldN %D maxDof %D\n",childId,numRows,numCols,refPointFieldN[childId - pRefStart][f], maxDof);
4271: #endif
4272:           for (i = 0; i < numRows; i++) {
4273:             PetscScalar val = 0.;
4274:             for (j = 0; j < numCols; j++) {
4275:               val += childMat[i * numCols + j] * cVal[j];
4276:             }
4277:             rVal[i] = val;
4278:           }
4279:           if (f == fvField && p >= cellStart && p < cellEnd) {
4280:             PetscReal   centroid[3];
4281:             PetscScalar diff[3];
4282:             const PetscScalar *parentCentroid = &fvGradData[0];
4283:             const PetscScalar *gradient       = &fvGradData[dim];

4285:             DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);
4286:             for (i = 0; i < dim; i++) {
4287:               diff[i] = centroid[i] - parentCentroid[i];
4288:             }
4289:             for (i = 0; i < numFVcomps; i++) {
4290:               PetscScalar val = 0.;

4292:               for (j = 0; j < dim; j++) {
4293:                 val += gradient[dim * i + j] * diff[j];
4294:               }
4295:               rVal[i] += val;
4296:             }
4297:           }
4298:           VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);
4299:         }
4300:       }
4301:     }
4302:     DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
4303:     DMRestoreWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);
4304:     DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
4305:   }
4306:   PetscFree(leafValues);
4307:   PetscSectionDestroy(&leafValuesSec);
4308:   PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);
4309:   ISRestoreIndices(aIS,&anchors);
4310:   return(0);
4311: }

4313: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4314: {
4315:   DM             refTree;
4316:   PetscSection   multiRootSec, rootIndicesSec;
4317:   PetscSection   globalCoarse, globalFine;
4318:   PetscSection   localCoarse, localFine;
4319:   PetscSection   cSecRef;
4320:   PetscInt       *parentIndices, pRefStart, pRefEnd;
4321:   PetscScalar    *rootValues, *parentValues;
4322:   Mat            injRef;
4323:   PetscInt       numFields, maxDof;
4324:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
4325:   PetscInt       *offsets, *offsetsCopy, *rowOffsets;
4326:   PetscLayout    rowMap, colMap;
4327:   PetscInt       rowStart, rowEnd, colStart, colEnd;
4328:   PetscScalar    ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */


4333:   /* get the templates for the fine-to-coarse injection from the reference tree */
4334:   VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
4335:   VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
4336:   DMPlexGetReferenceTree(coarse,&refTree);
4337:   DMCopyDisc(coarse,refTree);
4338:   DMGetDefaultConstraints(refTree,&cSecRef,NULL);
4339:   PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
4340:   DMPlexReferenceTreeGetInjector(refTree,&injRef);

4342:   DMPlexGetChart(fine,&pStartF,&pEndF);
4343:   DMGetSection(fine,&localFine);
4344:   DMGetGlobalSection(fine,&globalFine);
4345:   PetscSectionGetNumFields(localFine,&numFields);
4346:   DMPlexGetChart(coarse,&pStartC,&pEndC);
4347:   DMGetSection(coarse,&localCoarse);
4348:   DMGetGlobalSection(coarse,&globalCoarse);
4349:   PetscSectionGetMaxDof(localCoarse,&maxDof);
4350:   {
4351:     PetscInt maxFields = PetscMax(1,numFields) + 1;
4352:     PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
4353:   }

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

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

4359:   /* count indices */
4360:   VecGetLayout(vecFine,&colMap);
4361:   VecGetLayout(vecCoarse,&rowMap);
4362:   PetscLayoutSetUp(rowMap);
4363:   PetscLayoutSetUp(colMap);
4364:   PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
4365:   PetscLayoutGetRange(colMap,&colStart,&colEnd);
4366:   /* insert values */
4367:   DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
4368:   for (p = pStartC; p < pEndC; p++) {
4369:     PetscInt  numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4370:     PetscBool contribute = PETSC_FALSE;

4372:     PetscSectionGetDof(globalCoarse,p,&dof);
4373:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
4374:     if ((dof - cdof) <= 0) continue;
4375:     PetscSectionGetDof(localCoarse,p,&dof);
4376:     PetscSectionGetOffset(globalCoarse,p,&gOff);

4378:     rowOffsets[0] = 0;
4379:     offsetsCopy[0] = 0;
4380:     if (numFields) {
4381:       PetscInt f;

4383:       for (f = 0; f < numFields; f++) {
4384:         PetscInt fDof;
4385:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
4386:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
4387:       }
4388:       DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);
4389:     } else {
4390:       DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);
4391:       rowOffsets[1] = offsetsCopy[0];
4392:     }

4394:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
4395:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
4396:     leafEnd = leafStart + numLeaves;
4397:     for (l = 0; l < dof; l++) parentValues[l] = 0.;
4398:     for (l = leafStart; l < leafEnd; l++) {
4399:       PetscInt numIndices, childId, offset;
4400:       const PetscScalar *childValues;

4402:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
4403:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
4404:       childId = (PetscInt) PetscRealPart(rootValues[offset++]);
4405:       childValues = &rootValues[offset];
4406:       numIndices--;

4408:       if (childId == -2) { /* skip */
4409:         continue;
4410:       } else if (childId == -1) { /* equivalent points: scatter */
4411:         PetscInt m;

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

4418:         contribute = PETSC_TRUE;
4419:         DMPlexGetTreeParent(refTree,childId,&parentId,NULL);

4421:         lim = PetscMax(1,numFields);
4422:         offsets[0] = 0;
4423:         if (numFields) {
4424:           PetscInt f;

4426:           for (f = 0; f < numFields; f++) {
4427:             PetscInt fDof;
4428:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

4430:             offsets[f + 1] = fDof + offsets[f];
4431:           }
4432:         }
4433:         else {
4434:           PetscInt cDof;

4436:           PetscSectionGetDof(cSecRef,childId,&cDof);
4437:           offsets[1] = cDof;
4438:         }
4439:         for (f = 0; f < lim; f++) {
4440:           PetscScalar       *childMat   = &childrenMats[childId - pRefStart][f][0];
4441:           PetscInt          n           = offsets[f+1]-offsets[f];
4442:           PetscInt          m           = rowOffsets[f+1]-rowOffsets[f];
4443:           PetscInt          i, j;
4444:           const PetscScalar *colValues  = &childValues[offsets[f]];

4446:           for (i = 0; i < m; i++) {
4447:             PetscScalar val = 0.;
4448:             for (j = 0; j < n; j++) {
4449:               val += childMat[n * i + j] * colValues[j];
4450:             }
4451:             parentValues[rowOffsets[f] + i] += val;
4452:           }
4453:         }
4454:       }
4455:     }
4456:     if (contribute) {VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);}
4457:   }
4458:   PetscSectionDestroy(&multiRootSec);
4459:   PetscSectionDestroy(&rootIndicesSec);
4460:   PetscFree2(parentIndices,parentValues);
4461:   DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
4462:   PetscFree(rootValues);
4463:   PetscFree3(offsets,offsetsCopy,rowOffsets);
4464:   return(0);
4465: }

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

4472:   collective

4474:   Input Parameters:
4475: + dmIn        - The DMPlex mesh for the input vector
4476: . vecIn       - The input vector
4477: . sfRefine    - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in
4478:                 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn
4479: . sfCoarsen   - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in
4480:                 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn
4481: . cidsRefine  - The childIds of the points in dmOut.  These childIds relate back to the reference tree: childid[j] = k implies
4482:                 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference
4483:                 tree was refined from its parent.  childid[j] = -1 indicates that the point j in dmOut is exactly
4484:                 equivalent to its root in dmIn, so no interpolation is necessary.  childid[j] = -2 indicates that this
4485:                 point j in dmOut is not a leaf of sfRefine.
4486: . cidsCoarsen - The childIds of the points in dmIn.  These childIds relate back to the reference tree: childid[j] = k implies
4487:                 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference
4488:                 tree coarsens to its parent.  childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen.
4489: . useBCs      - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer.
4490: - time        - Used if boundary values are time dependent.

4492:   Output Parameters:
4493: . vecOut      - Using interpolation and injection operators calculated on the reference tree, the transfered
4494:                 projection of vecIn from dmIn to dmOut.  Note that any field discretized with a PetscFV finite volume
4495:                 method that uses gradient reconstruction will use reconstructed gradients when interpolating from
4496:                 coarse points to fine points.

4498:   Level: developer

4500: .seealso(): DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients()
4501: @*/
4502: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4503: {

4507:   VecSet(vecOut,0.0);
4508:   if (sfRefine) {
4509:     Vec vecInLocal;
4510:     DM  dmGrad = NULL;
4511:     Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;

4513:     DMGetLocalVector(dmIn,&vecInLocal);
4514:     VecSet(vecInLocal,0.0);
4515:     {
4516:       PetscInt  numFields, i;

4518:       DMGetNumFields(dmIn, &numFields);
4519:       for (i = 0; i < numFields; i++) {
4520:         PetscObject  obj;
4521:         PetscClassId classid;

4523:         DMGetField(dmIn, i, NULL, &obj);
4524:         PetscObjectGetClassId(obj, &classid);
4525:         if (classid == PETSCFV_CLASSID) {
4526:           DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);
4527:           break;
4528:         }
4529:       }
4530:     }
4531:     if (useBCs) {
4532:       DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);
4533:     }
4534:     DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);
4535:     DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);
4536:     if (dmGrad) {
4537:       DMGetGlobalVector(dmGrad,&grad);
4538:       DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);
4539:     }
4540:     DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);
4541:     DMRestoreLocalVector(dmIn,&vecInLocal);
4542:     if (dmGrad) {
4543:       DMRestoreGlobalVector(dmGrad,&grad);
4544:     }
4545:   }
4546:   if (sfCoarsen) {
4547:     DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);
4548:   }
4549:   VecAssemblyBegin(vecOut);
4550:   VecAssemblyEnd(vecOut);
4551:   return(0);
4552: }