Actual source code: plextree.c

petsc-3.9.4 2018-09-11
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_Union(DM K, DM Kref, const char *labelName, DM *ref)
212: {
213:   MPI_Comm       comm;
214:   PetscInt       dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs;
215:   PetscInt      *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts;
216:   DMLabel        identity, identityRef;
217:   PetscSection   unionSection, unionConeSection, parentSection;
218:   PetscScalar   *unionCoords;
219:   IS             perm;

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

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

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

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

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

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

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

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

356:     DMGetCoordinateSection(K, &KcoordsSec);
357:     DMGetCoordinatesLocal(K, &KcoordsVec);
358:     DMGetCoordinateSection(Kref, &KrefCoordsSec);
359:     DMGetCoordinatesLocal(Kref, &KrefCoordsVec);

361:     numVerts = numDimPoints[0];
362:     PetscMalloc1(numVerts * dim,&unionCoords);
363:     DMPlexGetDepthStratum(K,0,&vStart,&vEnd);

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

397:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
398:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
399:     if (uDof) {
400:       PetscSectionSetDof(parentSection,uOff,1);
401:     }
402:   }
403:   PetscSectionSetUp(parentSection);
404:   PetscSectionGetStorageSize(parentSection,&parentSize);
405:   PetscMalloc2(parentSize,&parents,parentSize,&childIDs);
406:   for (p = pRefStart; p < pRefEnd; p++) {
407:     PetscInt uDof, uOff;

409:     PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
410:     PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
411:     if (uDof) {
412:       PetscInt pOff, parent, parentU;
413:       PetscSectionGetOffset(parentSection,uOff,&pOff);
414:       DMLabelGetValue(identityRef,p,&parent);
415:       PetscSectionGetOffset(unionSection, parent - pStart,&parentU);
416:       parents[pOff] = parentU;
417:       childIDs[pOff] = uOff;
418:     }
419:   }
420:   DMPlexSetTree_Internal(*ref,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);
421:   PetscSectionDestroy(&parentSection);
422:   PetscFree2(parents,childIDs);

424:   /* clean up */
425:   PetscSectionDestroy(&unionSection);
426:   PetscSectionDestroy(&unionConeSection);
427:   ISDestroy(&perm);
428:   PetscFree(unionCoords);
429:   PetscFree2(unionCones,unionOrientations);
430:   PetscFree2(coneSizes,numDimPoints);
431:   return(0);
432: }

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

437:   Collective on comm

439:   Input Parameters:
440: + comm    - the MPI communicator
441: . dim     - the spatial dimension
442: - simplex - Flag for simplex, otherwise use a tensor-product cell

444:   Output Parameters:
445: . ref     - the reference tree DMPlex object

447:   Level: intermediate

449: .keywords: reference cell
450: .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree()
451: @*/
452: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
453: {
454:   DM_Plex       *mesh;
455:   DM             K, Kref;
456:   PetscInt       p, pStart, pEnd;
457:   DMLabel        identity;

461: #if 1
462:   comm = PETSC_COMM_SELF;
463: #endif
464:   /* create a reference element */
465:   DMPlexCreateReferenceCell(comm, dim, simplex, &K);
466:   DMCreateLabel(K, "identity");
467:   DMGetLabel(K, "identity", &identity);
468:   DMPlexGetChart(K, &pStart, &pEnd);
469:   for (p = pStart; p < pEnd; p++) {
470:     DMLabelSetValue(identity, p, p);
471:   }
472:   /* refine it */
473:   DMRefine(K,comm,&Kref);

475:   /* the reference tree is the union of these two, without duplicating
476:    * points that appear in both */
477:   DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref);
478:   mesh = (DM_Plex *) (*ref)->data;
479:   mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
480:   DMDestroy(&K);
481:   DMDestroy(&Kref);
482:   return(0);
483: }

485: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
486: {
487:   DM_Plex        *mesh = (DM_Plex *)dm->data;
488:   PetscSection   childSec, pSec;
489:   PetscInt       p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT;
490:   PetscInt       *offsets, *children, pStart, pEnd;

495:   PetscSectionDestroy(&mesh->childSection);
496:   PetscFree(mesh->children);
497:   pSec = mesh->parentSection;
498:   if (!pSec) return(0);
499:   PetscSectionGetStorageSize(pSec,&pSize);
500:   for (p = 0; p < pSize; p++) {
501:     PetscInt par = mesh->parents[p];

503:     parMax = PetscMax(parMax,par+1);
504:     parMin = PetscMin(parMin,par);
505:   }
506:   if (parMin > parMax) {
507:     parMin = -1;
508:     parMax = -1;
509:   }
510:   PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);
511:   PetscSectionSetChart(childSec,parMin,parMax);
512:   for (p = 0; p < pSize; p++) {
513:     PetscInt par = mesh->parents[p];

515:     PetscSectionAddDof(childSec,par,1);
516:   }
517:   PetscSectionSetUp(childSec);
518:   PetscSectionGetStorageSize(childSec,&cSize);
519:   PetscMalloc1(cSize,&children);
520:   PetscCalloc1(parMax-parMin,&offsets);
521:   PetscSectionGetChart(pSec,&pStart,&pEnd);
522:   for (p = pStart; p < pEnd; p++) {
523:     PetscInt dof, off, i;

525:     PetscSectionGetDof(pSec,p,&dof);
526:     PetscSectionGetOffset(pSec,p,&off);
527:     for (i = 0; i < dof; i++) {
528:       PetscInt par = mesh->parents[off + i], cOff;

530:       PetscSectionGetOffset(childSec,par,&cOff);
531:       children[cOff + offsets[par-parMin]++] = p;
532:     }
533:   }
534:   mesh->childSection = childSec;
535:   mesh->children = children;
536:   PetscFree(offsets);
537:   return(0);
538: }

540: static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
541: {
542:   PetscInt       pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
543:   const PetscInt *vals;
544:   PetscSection   secNew;
545:   PetscBool      anyNew, globalAnyNew;
546:   PetscBool      compress;

550:   PetscSectionGetChart(section,&pStart,&pEnd);
551:   ISGetLocalSize(is,&size);
552:   ISGetIndices(is,&vals);
553:   PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);
554:   PetscSectionSetChart(secNew,pStart,pEnd);
555:   for (i = 0; i < size; i++) {
556:     PetscInt dof;

558:     p = vals[i];
559:     if (p < pStart || p >= pEnd) continue;
560:     PetscSectionGetDof(section, p, &dof);
561:     if (dof) break;
562:   }
563:   if (i == size) {
564:     PetscSectionSetUp(secNew);
565:     anyNew   = PETSC_FALSE;
566:     compress = PETSC_FALSE;
567:     sizeNew  = 0;
568:   }
569:   else {
570:     anyNew = PETSC_TRUE;
571:     for (p = pStart; p < pEnd; p++) {
572:       PetscInt dof, off;

574:       PetscSectionGetDof(section, p, &dof);
575:       PetscSectionGetOffset(section, p, &off);
576:       for (i = 0; i < dof; i++) {
577:         PetscInt q = vals[off + i], qDof = 0;

579:         if (q >= pStart && q < pEnd) {
580:           PetscSectionGetDof(section, q, &qDof);
581:         }
582:         if (qDof) {
583:           PetscSectionAddDof(secNew, p, qDof);
584:         }
585:         else {
586:           PetscSectionAddDof(secNew, p, 1);
587:         }
588:       }
589:     }
590:     PetscSectionSetUp(secNew);
591:     PetscSectionGetStorageSize(secNew,&sizeNew);
592:     PetscMalloc1(sizeNew,&valsNew);
593:     compress = PETSC_FALSE;
594:     for (p = pStart; p < pEnd; p++) {
595:       PetscInt dof, off, count, offNew, dofNew;

597:       PetscSectionGetDof(section, p, &dof);
598:       PetscSectionGetOffset(section, p, &off);
599:       PetscSectionGetDof(secNew, p, &dofNew);
600:       PetscSectionGetOffset(secNew, p, &offNew);
601:       count = 0;
602:       for (i = 0; i < dof; i++) {
603:         PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;

605:         if (q >= pStart && q < pEnd) {
606:           PetscSectionGetDof(section, q, &qDof);
607:           PetscSectionGetOffset(section, q, &qOff);
608:         }
609:         if (qDof) {
610:           PetscInt oldCount = count;

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

615:             for (k = 0; k < oldCount; k++) {
616:               if (valsNew[offNew + k] == r) {
617:                 break;
618:               }
619:             }
620:             if (k == oldCount) {
621:               valsNew[offNew + count++] = r;
622:             }
623:           }
624:         }
625:         else {
626:           PetscInt k, oldCount = count;

628:           for (k = 0; k < oldCount; k++) {
629:             if (valsNew[offNew + k] == q) {
630:               break;
631:             }
632:           }
633:           if (k == oldCount) {
634:             valsNew[offNew + count++] = q;
635:           }
636:         }
637:       }
638:       if (count < dofNew) {
639:         PetscSectionSetDof(secNew, p, count);
640:         compress = PETSC_TRUE;
641:       }
642:     }
643:   }
644:   ISRestoreIndices(is,&vals);
645:   MPIU_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
646:   if (!globalAnyNew) {
647:     PetscSectionDestroy(&secNew);
648:     *sectionNew = NULL;
649:     *isNew = NULL;
650:   }
651:   else {
652:     PetscBool globalCompress;

654:     MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
655:     if (compress) {
656:       PetscSection secComp;
657:       PetscInt *valsComp = NULL;

659:       PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);
660:       PetscSectionSetChart(secComp,pStart,pEnd);
661:       for (p = pStart; p < pEnd; p++) {
662:         PetscInt dof;

664:         PetscSectionGetDof(secNew, p, &dof);
665:         PetscSectionSetDof(secComp, p, dof);
666:       }
667:       PetscSectionSetUp(secComp);
668:       PetscSectionGetStorageSize(secComp,&sizeNew);
669:       PetscMalloc1(sizeNew,&valsComp);
670:       for (p = pStart; p < pEnd; p++) {
671:         PetscInt dof, off, offNew, j;

673:         PetscSectionGetDof(secNew, p, &dof);
674:         PetscSectionGetOffset(secNew, p, &off);
675:         PetscSectionGetOffset(secComp, p, &offNew);
676:         for (j = 0; j < dof; j++) {
677:           valsComp[offNew + j] = valsNew[off + j];
678:         }
679:       }
680:       PetscSectionDestroy(&secNew);
681:       secNew  = secComp;
682:       PetscFree(valsNew);
683:       valsNew = valsComp;
684:     }
685:     ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);
686:   }
687:   return(0);
688: }

690: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
691: {
692:   PetscInt       p, pStart, pEnd, *anchors, size;
693:   PetscInt       aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
694:   PetscSection   aSec;
695:   DMLabel        canonLabel;
696:   IS             aIS;

701:   DMPlexGetChart(dm,&pStart,&pEnd);
702:   DMGetLabel(dm,"canonical",&canonLabel);
703:   for (p = pStart; p < pEnd; p++) {
704:     PetscInt parent;

706:     if (canonLabel) {
707:       PetscInt canon;

709:       DMLabelGetValue(canonLabel,p,&canon);
710:       if (p != canon) continue;
711:     }
712:     DMPlexGetTreeParent(dm,p,&parent,NULL);
713:     if (parent != p) {
714:       aMin = PetscMin(aMin,p);
715:       aMax = PetscMax(aMax,p+1);
716:     }
717:   }
718:   if (aMin > aMax) {
719:     aMin = -1;
720:     aMax = -1;
721:   }
722:   PetscSectionCreate(PETSC_COMM_SELF,&aSec);
723:   PetscSectionSetChart(aSec,aMin,aMax);
724:   for (p = aMin; p < aMax; p++) {
725:     PetscInt parent, ancestor = p;

727:     if (canonLabel) {
728:       PetscInt canon;

730:       DMLabelGetValue(canonLabel,p,&canon);
731:       if (p != canon) continue;
732:     }
733:     DMPlexGetTreeParent(dm,p,&parent,NULL);
734:     while (parent != ancestor) {
735:       ancestor = parent;
736:       DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
737:     }
738:     if (ancestor != p) {
739:       PetscInt closureSize, *closure = NULL;

741:       DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
742:       PetscSectionSetDof(aSec,p,closureSize);
743:       DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
744:     }
745:   }
746:   PetscSectionSetUp(aSec);
747:   PetscSectionGetStorageSize(aSec,&size);
748:   PetscMalloc1(size,&anchors);
749:   for (p = aMin; p < aMax; p++) {
750:     PetscInt parent, ancestor = p;

752:     if (canonLabel) {
753:       PetscInt canon;

755:       DMLabelGetValue(canonLabel,p,&canon);
756:       if (p != canon) continue;
757:     }
758:     DMPlexGetTreeParent(dm,p,&parent,NULL);
759:     while (parent != ancestor) {
760:       ancestor = parent;
761:       DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
762:     }
763:     if (ancestor != p) {
764:       PetscInt j, closureSize, *closure = NULL, aOff;

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

768:       DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
769:       for (j = 0; j < closureSize; j++) {
770:         anchors[aOff + j] = closure[2*j];
771:       }
772:       DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
773:     }
774:   }
775:   ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);
776:   {
777:     PetscSection aSecNew = aSec;
778:     IS           aISNew  = aIS;

780:     PetscObjectReference((PetscObject)aSec);
781:     PetscObjectReference((PetscObject)aIS);
782:     while (aSecNew) {
783:       PetscSectionDestroy(&aSec);
784:       ISDestroy(&aIS);
785:       aSec    = aSecNew;
786:       aIS     = aISNew;
787:       aSecNew = NULL;
788:       aISNew  = NULL;
789:       AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);
790:     }
791:   }
792:   DMPlexSetAnchors(dm,aSec,aIS);
793:   PetscSectionDestroy(&aSec);
794:   ISDestroy(&aIS);
795:   return(0);
796: }

798: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp)
799: {

803:   if (numTrueSupp[p] == -1) {
804:     PetscInt i, alldof;
805:     const PetscInt *supp;
806:     PetscInt count = 0;

808:     DMPlexGetSupportSize(dm,p,&alldof);
809:     DMPlexGetSupport(dm,p,&supp);
810:     for (i = 0; i < alldof; i++) {
811:       PetscInt q = supp[i], numCones, j;
812:       const PetscInt *cone;

814:       DMPlexGetConeSize(dm,q,&numCones);
815:       DMPlexGetCone(dm,q,&cone);
816:       for (j = 0; j < numCones; j++) {
817:         if (cone[j] == p) break;
818:       }
819:       if (j < numCones) count++;
820:     }
821:     numTrueSupp[p] = count;
822:   }
823:   *dof = numTrueSupp[p];
824:   return(0);
825: }

827: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
828: {
829:   DM_Plex *mesh = (DM_Plex *)dm->data;
830:   PetscSection newSupportSection;
831:   PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth;
832:   PetscInt *numTrueSupp;
833:   PetscInt *offsets;

838:   /* symmetrize the hierarchy */
839:   DMPlexGetDepth(dm,&depth);
840:   PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);
841:   DMPlexGetChart(dm,&pStart,&pEnd);
842:   PetscSectionSetChart(newSupportSection,pStart,pEnd);
843:   PetscCalloc1(pEnd,&offsets);
844:   PetscMalloc1(pEnd,&numTrueSupp);
845:   for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
846:   /* if a point is in the (true) support of q, it should be in the support of
847:    * parent(q) */
848:   for (d = 0; d <= depth; d++) {
849:     DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
850:     for (p = pStart; p < pEnd; ++p) {
851:       PetscInt dof, q, qdof, parent;

853:       DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);
854:       PetscSectionAddDof(newSupportSection, p, dof);
855:       q    = p;
856:       DMPlexGetTreeParent(dm,q,&parent,NULL);
857:       while (parent != q && parent >= pStart && parent < pEnd) {
858:         q = parent;

860:         DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);
861:         PetscSectionAddDof(newSupportSection,p,qdof);
862:         PetscSectionAddDof(newSupportSection,q,dof);
863:         DMPlexGetTreeParent(dm,q,&parent,NULL);
864:       }
865:     }
866:   }
867:   PetscSectionSetUp(newSupportSection);
868:   PetscSectionGetStorageSize(newSupportSection,&newSize);
869:   PetscMalloc1(newSize,&newSupports);
870:   for (d = 0; d <= depth; d++) {
871:     DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
872:     for (p = pStart; p < pEnd; p++) {
873:       PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;

875:       PetscSectionGetDof(mesh->supportSection, p, &dof);
876:       PetscSectionGetOffset(mesh->supportSection, p, &off);
877:       PetscSectionGetDof(newSupportSection, p, &newDof);
878:       PetscSectionGetOffset(newSupportSection, p, &newOff);
879:       for (i = 0; i < dof; i++) {
880:         PetscInt numCones, j;
881:         const PetscInt *cone;
882:         PetscInt q = mesh->supports[off + i];

884:         DMPlexGetConeSize(dm,q,&numCones);
885:         DMPlexGetCone(dm,q,&cone);
886:         for (j = 0; j < numCones; j++) {
887:           if (cone[j] == p) break;
888:         }
889:         if (j < numCones) newSupports[newOff+offsets[p]++] = q;
890:       }
891:       mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof);

893:       q    = p;
894:       DMPlexGetTreeParent(dm,q,&parent,NULL);
895:       while (parent != q && parent >= pStart && parent < pEnd) {
896:         q = parent;
897:         PetscSectionGetDof(mesh->supportSection, q, &qdof);
898:         PetscSectionGetOffset(mesh->supportSection, q, &qoff);
899:         PetscSectionGetOffset(newSupportSection, q, &newqOff);
900:         for (i = 0; i < qdof; i++) {
901:           PetscInt numCones, j;
902:           const PetscInt *cone;
903:           PetscInt r = mesh->supports[qoff + i];

905:           DMPlexGetConeSize(dm,r,&numCones);
906:           DMPlexGetCone(dm,r,&cone);
907:           for (j = 0; j < numCones; j++) {
908:             if (cone[j] == q) break;
909:           }
910:           if (j < numCones) newSupports[newOff+offsets[p]++] = r;
911:         }
912:         for (i = 0; i < dof; i++) {
913:           PetscInt numCones, j;
914:           const PetscInt *cone;
915:           PetscInt r = mesh->supports[off + i];

917:           DMPlexGetConeSize(dm,r,&numCones);
918:           DMPlexGetCone(dm,r,&cone);
919:           for (j = 0; j < numCones; j++) {
920:             if (cone[j] == p) break;
921:           }
922:           if (j < numCones) newSupports[newqOff+offsets[q]++] = r;
923:         }
924:         DMPlexGetTreeParent(dm,q,&parent,NULL);
925:       }
926:     }
927:   }
928:   PetscSectionDestroy(&mesh->supportSection);
929:   mesh->supportSection = newSupportSection;
930:   PetscFree(mesh->supports);
931:   mesh->supports = newSupports;
932:   PetscFree(offsets);
933:   PetscFree(numTrueSupp);

935:   return(0);
936: }

938: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat);
939: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat);

941: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
942: {
943:   DM_Plex       *mesh = (DM_Plex *)dm->data;
944:   DM             refTree;
945:   PetscInt       size;

951:   PetscObjectReference((PetscObject)parentSection);
952:   PetscSectionDestroy(&mesh->parentSection);
953:   mesh->parentSection = parentSection;
954:   PetscSectionGetStorageSize(parentSection,&size);
955:   if (parents != mesh->parents) {
956:     PetscFree(mesh->parents);
957:     PetscMalloc1(size,&mesh->parents);
958:     PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));
959:   }
960:   if (childIDs != mesh->childIDs) {
961:     PetscFree(mesh->childIDs);
962:     PetscMalloc1(size,&mesh->childIDs);
963:     PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));
964:   }
965:   DMPlexGetReferenceTree(dm,&refTree);
966:   if (refTree) {
967:     DMLabel canonLabel;

969:     DMGetLabel(refTree,"canonical",&canonLabel);
970:     if (canonLabel) {
971:       PetscInt i;

973:       for (i = 0; i < size; i++) {
974:         PetscInt canon;
975:         DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);
976:         if (canon >= 0) {
977:           mesh->childIDs[i] = canon;
978:         }
979:       }
980:     }
981:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
982:   }
983:   else {
984:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
985:   }
986:   DMPlexTreeSymmetrize(dm);
987:   if (computeCanonical) {
988:     PetscInt d, dim;

990:     /* add the canonical label */
991:     DMGetDimension(dm,&dim);
992:     DMCreateLabel(dm,"canonical");
993:     for (d = 0; d <= dim; d++) {
994:       PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
995:       const PetscInt *cChildren;

997:       DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
998:       for (p = dStart; p < dEnd; p++) {
999:         DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);
1000:         if (cNumChildren) {
1001:           canon = p;
1002:           break;
1003:         }
1004:       }
1005:       if (canon == -1) continue;
1006:       for (p = dStart; p < dEnd; p++) {
1007:         PetscInt numChildren, i;
1008:         const PetscInt *children;

1010:         DMPlexGetTreeChildren(dm,p,&numChildren,&children);
1011:         if (numChildren) {
1012:           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);
1013:           DMSetLabelValue(dm,"canonical",p,canon);
1014:           for (i = 0; i < numChildren; i++) {
1015:             DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);
1016:           }
1017:         }
1018:       }
1019:     }
1020:   }
1021:   if (exchangeSupports) {
1022:     DMPlexTreeExchangeSupports(dm);
1023:   }
1024:   mesh->createanchors = DMPlexCreateAnchors_Tree;
1025:   /* reset anchors */
1026:   DMPlexSetAnchors(dm,NULL,NULL);
1027:   return(0);
1028: }

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

1035:   Collective on dm

1037:   Input Parameters:
1038: + dm - the DMPlex object
1039: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1040:                   offset indexes the parent and childID list; the reference count of parentSection is incremented
1041: . parents - a list of the point parents; copied, can be destroyed
1042: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1043:              the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed

1045:   Level: intermediate

1047: .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1048: @*/
1049: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
1050: {

1054:   DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);
1055:   return(0);
1056: }

1058: /*@
1059:   DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1060:   Collective on dm

1062:   Input Parameters:
1063: . dm - the DMPlex object

1065:   Output Parameters:
1066: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1067:                   offset indexes the parent and childID list
1068: . parents - a list of the point parents
1069: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1070:              the child corresponds to the point in the reference tree with index childID
1071: . childSection - the inverse of the parent section
1072: - children - a list of the point children

1074:   Level: intermediate

1076: .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1077: @*/
1078: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1079: {
1080:   DM_Plex        *mesh = (DM_Plex *)dm->data;

1084:   if (parentSection) *parentSection = mesh->parentSection;
1085:   if (parents)       *parents       = mesh->parents;
1086:   if (childIDs)      *childIDs      = mesh->childIDs;
1087:   if (childSection)  *childSection  = mesh->childSection;
1088:   if (children)      *children      = mesh->children;
1089:   return(0);
1090: }

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

1095:   Input Parameters:
1096: + dm - the DMPlex object
1097: - point - the query point

1099:   Output Parameters:
1100: + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent
1101: - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point
1102:             does not have a parent

1104:   Level: intermediate

1106: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren()
1107: @*/
1108: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1109: {
1110:   DM_Plex       *mesh = (DM_Plex *)dm->data;
1111:   PetscSection   pSec;

1116:   pSec = mesh->parentSection;
1117:   if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1118:     PetscInt dof;

1120:     PetscSectionGetDof (pSec, point, &dof);
1121:     if (dof) {
1122:       PetscInt off;

1124:       PetscSectionGetOffset (pSec, point, &off);
1125:       if (parent)  *parent = mesh->parents[off];
1126:       if (childID) *childID = mesh->childIDs[off];
1127:       return(0);
1128:     }
1129:   }
1130:   if (parent) {
1131:     *parent = point;
1132:   }
1133:   if (childID) {
1134:     *childID = 0;
1135:   }
1136:   return(0);
1137: }

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

1142:   Input Parameters:
1143: + dm - the DMPlex object
1144: - point - the query point

1146:   Output Parameters:
1147: + numChildren - if not NULL, set to the number of children
1148: - children - if not NULL, set to a list children, or set to NULL if the point has no children

1150:   Level: intermediate

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

1156: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent()
1157: @*/
1158: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1159: {
1160:   DM_Plex       *mesh = (DM_Plex *)dm->data;
1161:   PetscSection   childSec;
1162:   PetscInt       dof = 0;

1167:   childSec = mesh->childSection;
1168:   if (childSec && point >= childSec->pStart && point < childSec->pEnd) {
1169:     PetscSectionGetDof (childSec, point, &dof);
1170:   }
1171:   if (numChildren) *numChildren = dof;
1172:   if (children) {
1173:     if (dof) {
1174:       PetscInt off;

1176:       PetscSectionGetOffset (childSec, point, &off);
1177:       *children = &mesh->children[off];
1178:     }
1179:     else {
1180:       *children = NULL;
1181:     }
1182:   }
1183:   return(0);
1184: }

1186: 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)
1187: {
1188:   PetscInt       f, b, p, c, offset, qPoints;

1192:   PetscSpaceEvaluate(space,nPoints,points,work,NULL,NULL);
1193:   for (f = 0, offset = 0; f < nFunctionals; f++) {
1194:     qPoints = pointsPerFn[f];
1195:     for (b = 0; b < nBasis; b++) {
1196:       PetscScalar val = 0.;

1198:       for (p = 0; p < qPoints; p++) {
1199:         for (c = 0; c < nComps; c++) {
1200:           val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c];
1201:         }
1202:       }
1203:       MatSetValue(basisAtPoints,b,f,val,INSERT_VALUES);
1204:     }
1205:     offset += qPoints;
1206:   }
1207:   MatAssemblyBegin(basisAtPoints,MAT_FINAL_ASSEMBLY);
1208:   MatAssemblyEnd(basisAtPoints,MAT_FINAL_ASSEMBLY);
1209:   return(0);
1210: }

1212: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
1213: {
1214:   PetscDS        ds;
1215:   PetscInt       spdim;
1216:   PetscInt       numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
1217:   const PetscInt *anchors;
1218:   PetscSection   aSec;
1219:   PetscReal      *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
1220:   IS             aIS;

1224:   DMPlexGetChart(dm,&pStart,&pEnd);
1225:   DMGetDS(dm,&ds);
1226:   PetscDSGetNumFields(ds,&numFields);
1227:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1228:   DMPlexGetAnchors(dm,&aSec,&aIS);
1229:   ISGetIndices(aIS,&anchors);
1230:   PetscSectionGetChart(cSec,&conStart,&conEnd);
1231:   DMGetDimension(dm,&spdim);
1232:   PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);

1234:   for (f = 0; f < numFields; f++) {
1235:     PetscObject disc;
1236:     PetscClassId id;
1237:     PetscSpace     bspace;
1238:     PetscDualSpace dspace;
1239:     PetscInt i, j, k, nPoints, Nc, offset;
1240:     PetscInt fSize, maxDof;
1241:     PetscReal   *weights, *pointsRef, *pointsReal, *work;
1242:     PetscScalar *scwork, *X;
1243:     PetscInt  *sizes, *workIndRow, *workIndCol;
1244:     Mat Amat, Bmat, Xmat;
1245:     const PetscInt    *numDof  = NULL;
1246:     const PetscInt    ***perms = NULL;
1247:     const PetscScalar ***flips = NULL;

1249:     PetscDSGetDiscretization(ds,f,&disc);
1250:     PetscObjectGetClassId(disc,&id);
1251:     if (id == PETSCFE_CLASSID) {
1252:       PetscFE fe = (PetscFE) disc;

1254:       PetscFEGetBasisSpace(fe,&bspace);
1255:       PetscFEGetDualSpace(fe,&dspace);
1256:       PetscDualSpaceGetDimension(dspace,&fSize);
1257:       PetscFEGetNumComponents(fe,&Nc);
1258:     }
1259:     else if (id == PETSCFV_CLASSID) {
1260:       PetscFV fv = (PetscFV) disc;

1262:       PetscFVGetNumComponents(fv,&Nc);
1263:       PetscSpaceCreate(PetscObjectComm((PetscObject)fv),&bspace);
1264:       PetscSpaceSetType(bspace,PETSCSPACEPOLYNOMIAL);
1265:       PetscSpaceSetOrder(bspace,0);
1266:       PetscSpaceSetNumComponents(bspace,Nc);
1267:       PetscSpacePolynomialSetNumVariables(bspace,spdim);
1268:       PetscSpaceSetUp(bspace);
1269:       PetscFVGetDualSpace(fv,&dspace);
1270:       PetscDualSpaceGetDimension(dspace,&fSize);
1271:     }
1272:     else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1273:     PetscDualSpaceGetNumDof(dspace,&numDof);
1274:     for (i = 0, maxDof = 0; i <= spdim; i++) {maxDof = PetscMax(maxDof,numDof[i]);}
1275:     PetscDualSpaceGetSymmetries(dspace,&perms,&flips);

1277:     MatCreate(PETSC_COMM_SELF,&Amat);
1278:     MatSetSizes(Amat,fSize,fSize,fSize,fSize);
1279:     MatSetType(Amat,MATSEQDENSE);
1280:     MatSetUp(Amat);
1281:     MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);
1282:     MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);
1283:     nPoints = 0;
1284:     for (i = 0; i < fSize; i++) {
1285:       PetscInt        qPoints, thisNc;
1286:       PetscQuadrature quad;

1288:       PetscDualSpaceGetFunctional(dspace,i,&quad);
1289:       PetscQuadratureGetData(quad,NULL,&thisNc,&qPoints,NULL,NULL);
1290:       if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc);
1291:       nPoints += qPoints;
1292:     }
1293:     PetscMalloc7(fSize,&sizes,nPoints*Nc,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal,nPoints*fSize*Nc,&work,maxDof,&workIndRow,maxDof,&workIndCol);
1294:     PetscMalloc1(maxDof * maxDof,&scwork);
1295:     offset = 0;
1296:     for (i = 0; i < fSize; i++) {
1297:       PetscInt        qPoints;
1298:       const PetscReal    *p, *w;
1299:       PetscQuadrature quad;

1301:       PetscDualSpaceGetFunctional(dspace,i,&quad);
1302:       PetscQuadratureGetData(quad,NULL,NULL,&qPoints,&p,&w);
1303:       PetscMemcpy(weights+Nc*offset,w,Nc*qPoints*sizeof(*w));
1304:       PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));
1305:       sizes[i] = qPoints;
1306:       offset  += qPoints;
1307:     }
1308:     EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsRef,weights,work,Amat);
1309:     MatLUFactor(Amat,NULL,NULL,NULL);
1310:     for (c = cStart; c < cEnd; c++) {
1311:       PetscInt        parent;
1312:       PetscInt        closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1313:       PetscInt        *childOffsets, *parentOffsets;

1315:       DMPlexGetTreeParent(dm,c,&parent,NULL);
1316:       if (parent == c) continue;
1317:       DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1318:       for (i = 0; i < closureSize; i++) {
1319:         PetscInt p = closure[2*i];
1320:         PetscInt conDof;

1322:         if (p < conStart || p >= conEnd) continue;
1323:         if (numFields) {
1324:           PetscSectionGetFieldDof(cSec,p,f,&conDof);
1325:         }
1326:         else {
1327:           PetscSectionGetDof(cSec,p,&conDof);
1328:         }
1329:         if (conDof) break;
1330:       }
1331:       if (i == closureSize) {
1332:         DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1333:         continue;
1334:       }

1336:       DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ);
1337:       DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent);
1338:       for (i = 0; i < nPoints; i++) {
1339:         CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp);
1340:         CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);
1341:       }
1342:       EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsReal,weights,work,Bmat);
1343:       MatMatSolve(Amat,Bmat,Xmat);
1344:       MatDenseGetArray(Xmat,&X);
1345:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1346:       PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);
1347:       childOffsets[0] = 0;
1348:       for (i = 0; i < closureSize; i++) {
1349:         PetscInt p = closure[2*i];
1350:         PetscInt dof;

1352:         if (numFields) {
1353:           PetscSectionGetFieldDof(section,p,f,&dof);
1354:         }
1355:         else {
1356:           PetscSectionGetDof(section,p,&dof);
1357:         }
1358:         childOffsets[i+1]=childOffsets[i]+dof;
1359:       }
1360:       parentOffsets[0] = 0;
1361:       for (i = 0; i < closureSizeP; i++) {
1362:         PetscInt p = closureP[2*i];
1363:         PetscInt dof;

1365:         if (numFields) {
1366:           PetscSectionGetFieldDof(section,p,f,&dof);
1367:         }
1368:         else {
1369:           PetscSectionGetDof(section,p,&dof);
1370:         }
1371:         parentOffsets[i+1]=parentOffsets[i]+dof;
1372:       }
1373:       for (i = 0; i < closureSize; i++) {
1374:         PetscInt conDof, conOff, aDof, aOff, nWork;
1375:         PetscInt p = closure[2*i];
1376:         PetscInt o = closure[2*i+1];
1377:         const PetscInt    *perm;
1378:         const PetscScalar *flip;

1380:         if (p < conStart || p >= conEnd) continue;
1381:         if (numFields) {
1382:           PetscSectionGetFieldDof(cSec,p,f,&conDof);
1383:           PetscSectionGetFieldOffset(cSec,p,f,&conOff);
1384:         }
1385:         else {
1386:           PetscSectionGetDof(cSec,p,&conDof);
1387:           PetscSectionGetOffset(cSec,p,&conOff);
1388:         }
1389:         if (!conDof) continue;
1390:         perm  = (perms && perms[i]) ? perms[i][o] : NULL;
1391:         flip  = (flips && flips[i]) ? flips[i][o] : NULL;
1392:         PetscSectionGetDof(aSec,p,&aDof);
1393:         PetscSectionGetOffset(aSec,p,&aOff);
1394:         nWork = childOffsets[i+1]-childOffsets[i];
1395:         for (k = 0; k < aDof; k++) {
1396:           PetscInt a = anchors[aOff + k];
1397:           PetscInt aSecDof, aSecOff;

1399:           if (numFields) {
1400:             PetscSectionGetFieldDof(section,a,f,&aSecDof);
1401:             PetscSectionGetFieldOffset(section,a,f,&aSecOff);
1402:           }
1403:           else {
1404:             PetscSectionGetDof(section,a,&aSecDof);
1405:             PetscSectionGetOffset(section,a,&aSecOff);
1406:           }
1407:           if (!aSecDof) continue;

1409:           for (j = 0; j < closureSizeP; j++) {
1410:             PetscInt q = closureP[2*j];
1411:             PetscInt oq = closureP[2*j+1];

1413:             if (q == a) {
1414:               PetscInt           r, s, nWorkP;
1415:               const PetscInt    *permP;
1416:               const PetscScalar *flipP;

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

1470:   return(0);
1471: }

1473: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1474: {
1475:   Mat               refCmat;
1476:   PetscDS           ds;
1477:   PetscInt          numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1478:   PetscScalar       ***refPointFieldMats;
1479:   PetscSection      refConSec, refAnSec, refSection;
1480:   IS                refAnIS;
1481:   const PetscInt    *refAnchors;
1482:   const PetscInt    **perms;
1483:   const PetscScalar **flips;
1484:   PetscErrorCode    ierr;

1487:   DMGetDS(refTree,&ds);
1488:   PetscDSGetNumFields(ds,&numFields);
1489:   maxFields = PetscMax(1,numFields);
1490:   DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1491:   DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1492:   ISGetIndices(refAnIS,&refAnchors);
1493:   DMGetDefaultSection(refTree,&refSection);
1494:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1495:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
1496:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);
1497:   PetscSectionGetMaxDof(refConSec,&maxDof);
1498:   PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1499:   PetscMalloc1(maxDof,&rows);
1500:   PetscMalloc1(maxDof*maxAnDof,&cols);
1501:   for (p = pRefStart; p < pRefEnd; p++) {
1502:     PetscInt parent, closureSize, *closure = NULL, pDof;

1504:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
1505:     PetscSectionGetDof(refConSec,p,&pDof);
1506:     if (!pDof || parent == p) continue;

1508:     PetscMalloc1(maxFields,&refPointFieldMats[p-pRefStart]);
1509:     PetscCalloc1(maxFields,&refPointFieldN[p-pRefStart]);
1510:     DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1511:     for (f = 0; f < maxFields; f++) {
1512:       PetscInt cDof, cOff, numCols, r, i;

1514:       if (f < numFields) {
1515:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1516:         PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
1517:         PetscSectionGetFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);
1518:       } else {
1519:         PetscSectionGetDof(refConSec,p,&cDof);
1520:         PetscSectionGetOffset(refConSec,p,&cOff);
1521:         PetscSectionGetPointSyms(refSection,closureSize,closure,&perms,&flips);
1522:       }

1524:       for (r = 0; r < cDof; r++) {
1525:         rows[r] = cOff + r;
1526:       }
1527:       numCols = 0;
1528:       for (i = 0; i < closureSize; i++) {
1529:         PetscInt          q = closure[2*i];
1530:         PetscInt          aDof, aOff, j;
1531:         const PetscInt    *perm = perms ? perms[i] : NULL;

1533:         if (numFields) {
1534:           PetscSectionGetFieldDof(refSection,q,f,&aDof);
1535:           PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1536:         }
1537:         else {
1538:           PetscSectionGetDof(refSection,q,&aDof);
1539:           PetscSectionGetOffset(refSection,q,&aOff);
1540:         }

1542:         for (j = 0; j < aDof; j++) {
1543:           cols[numCols++] = aOff + (perm ? perm[j] : j);
1544:         }
1545:       }
1546:       refPointFieldN[p-pRefStart][f] = numCols;
1547:       PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
1548:       MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);
1549:       if (flips) {
1550:         PetscInt colOff = 0;

1552:         for (i = 0; i < closureSize; i++) {
1553:           PetscInt          q = closure[2*i];
1554:           PetscInt          aDof, aOff, j;
1555:           const PetscScalar *flip = flips ? flips[i] : NULL;

1557:           if (numFields) {
1558:             PetscSectionGetFieldDof(refSection,q,f,&aDof);
1559:             PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1560:           }
1561:           else {
1562:             PetscSectionGetDof(refSection,q,&aDof);
1563:             PetscSectionGetOffset(refSection,q,&aOff);
1564:           }
1565:           if (flip) {
1566:             PetscInt k;
1567:             for (k = 0; k < cDof; k++) {
1568:               for (j = 0; j < aDof; j++) {
1569:                 refPointFieldMats[p-pRefStart][f][k * numCols + colOff + j] *= flip[j];
1570:               }
1571:             }
1572:           }
1573:           colOff += aDof;
1574:         }
1575:       }
1576:       if (numFields) {
1577:         PetscSectionRestoreFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);
1578:       } else {
1579:         PetscSectionRestorePointSyms(refSection,closureSize,closure,&perms,&flips);
1580:       }
1581:     }
1582:     DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1583:   }
1584:   *childrenMats = refPointFieldMats;
1585:   *childrenN = refPointFieldN;
1586:   ISRestoreIndices(refAnIS,&refAnchors);
1587:   PetscFree(rows);
1588:   PetscFree(cols);
1589:   return(0);
1590: }

1592: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1593: {
1594:   PetscDS        ds;
1595:   PetscInt       **refPointFieldN;
1596:   PetscScalar    ***refPointFieldMats;
1597:   PetscInt       numFields, maxFields, pRefStart, pRefEnd, p, f;
1598:   PetscSection   refConSec;

1602:   refPointFieldN = *childrenN;
1603:   *childrenN = NULL;
1604:   refPointFieldMats = *childrenMats;
1605:   *childrenMats = NULL;
1606:   DMGetDS(refTree,&ds);
1607:   PetscDSGetNumFields(ds,&numFields);
1608:   maxFields = PetscMax(1,numFields);
1609:   DMGetDefaultConstraints(refTree,&refConSec,NULL);
1610:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1611:   for (p = pRefStart; p < pRefEnd; p++) {
1612:     PetscInt parent, pDof;

1614:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
1615:     PetscSectionGetDof(refConSec,p,&pDof);
1616:     if (!pDof || parent == p) continue;

1618:     for (f = 0; f < maxFields; f++) {
1619:       PetscInt cDof;

1621:       if (numFields) {
1622:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1623:       }
1624:       else {
1625:         PetscSectionGetDof(refConSec,p,&cDof);
1626:       }

1628:       PetscFree(refPointFieldMats[p - pRefStart][f]);
1629:     }
1630:     PetscFree(refPointFieldMats[p - pRefStart]);
1631:     PetscFree(refPointFieldN[p - pRefStart]);
1632:   }
1633:   PetscFree(refPointFieldMats);
1634:   PetscFree(refPointFieldN);
1635:   return(0);
1636: }

1638: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1639: {
1640:   DM             refTree;
1641:   PetscDS        ds;
1642:   Mat            refCmat;
1643:   PetscInt       numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1644:   PetscScalar ***refPointFieldMats, *pointWork;
1645:   PetscSection   refConSec, refAnSec, anSec;
1646:   IS             refAnIS, anIS;
1647:   const PetscInt *anchors;

1652:   DMGetDS(dm,&ds);
1653:   PetscDSGetNumFields(ds,&numFields);
1654:   maxFields = PetscMax(1,numFields);
1655:   DMPlexGetReferenceTree(dm,&refTree);
1656:   DMSetDS(refTree,ds);
1657:   DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1658:   DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1659:   DMPlexGetAnchors(dm,&anSec,&anIS);
1660:   ISGetIndices(anIS,&anchors);
1661:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1662:   PetscSectionGetChart(conSec,&conStart,&conEnd);
1663:   PetscSectionGetMaxDof(refConSec,&maxDof);
1664:   PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1665:   PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);

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

1670:   /* step 2: compute the preorder */
1671:   DMPlexGetChart(dm,&pStart,&pEnd);
1672:   PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);
1673:   for (p = pStart; p < pEnd; p++) {
1674:     perm[p - pStart] = p;
1675:     iperm[p - pStart] = p-pStart;
1676:   }
1677:   for (p = 0; p < pEnd - pStart;) {
1678:     PetscInt point = perm[p];
1679:     PetscInt parent;

1681:     DMPlexGetTreeParent(dm,point,&parent,NULL);
1682:     if (parent == point) {
1683:       p++;
1684:     }
1685:     else {
1686:       PetscInt size, closureSize, *closure = NULL, i;

1688:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1689:       for (i = 0; i < closureSize; i++) {
1690:         PetscInt q = closure[2*i];
1691:         if (iperm[q-pStart] > iperm[point-pStart]) {
1692:           /* swap */
1693:           perm[p]               = q;
1694:           perm[iperm[q-pStart]] = point;
1695:           iperm[point-pStart]   = iperm[q-pStart];
1696:           iperm[q-pStart]       = p;
1697:           break;
1698:         }
1699:       }
1700:       size = closureSize;
1701:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1702:       if (i == size) {
1703:         p++;
1704:       }
1705:     }
1706:   }

1708:   /* step 3: fill the constraint matrix */
1709:   /* we are going to use a preorder progressive fill strategy.  Mat doesn't
1710:    * allow progressive fill without assembly, so we are going to set up the
1711:    * values outside of the Mat first.
1712:    */
1713:   {
1714:     PetscInt nRows, row, nnz;
1715:     PetscBool done;
1716:     const PetscInt *ia, *ja;
1717:     PetscScalar *vals;

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

1729:       DMPlexGetTreeParent(dm,point,&parent,&childid);
1730:       if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1731:       PetscSectionGetDof(conSec,point,&pointDof);
1732:       if (!pointDof) continue;
1733:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1734:       for (f = 0; f < maxFields; f++) {
1735:         PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
1736:         PetscScalar *pointMat;
1737:         const PetscInt    **perms;
1738:         const PetscScalar **flips;

1740:         if (numFields) {
1741:           PetscSectionGetFieldDof(conSec,point,f,&cDof);
1742:           PetscSectionGetFieldOffset(conSec,point,f,&cOff);
1743:         }
1744:         else {
1745:           PetscSectionGetDof(conSec,point,&cDof);
1746:           PetscSectionGetOffset(conSec,point,&cOff);
1747:         }
1748:         if (!cDof) continue;
1749:         if (numFields) {PetscSectionGetFieldPointSyms(section,f,closureSize,closure,&perms,&flips);}
1750:         else           {PetscSectionGetPointSyms(section,closureSize,closure,&perms,&flips);}

1752:         /* make sure that every row for this point is the same size */
1753: #if defined(PETSC_USE_DEBUG)
1754:         for (r = 0; r < cDof; r++) {
1755:           if (cDof > 1 && r) {
1756:             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]));
1757:           }
1758:         }
1759: #endif
1760:         /* zero rows */
1761:         for (i = ia[cOff] ; i< ia[cOff+cDof];i++) {
1762:           vals[i] = 0.;
1763:         }
1764:         matOffset = ia[cOff];
1765:         numFillCols = ia[cOff+1] - matOffset;
1766:         pointMat = refPointFieldMats[childid-pRefStart][f];
1767:         numCols = refPointFieldN[childid-pRefStart][f];
1768:         offset = 0;
1769:         for (i = 0; i < closureSize; i++) {
1770:           PetscInt q = closure[2*i];
1771:           PetscInt aDof, aOff, j, k, qConDof, qConOff;
1772:           const PetscInt    *perm = perms ? perms[i] : NULL;
1773:           const PetscScalar *flip = flips ? flips[i] : NULL;

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

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

1836:                 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
1837:               }
1838:             }
1839:           }
1840:           offset += aDof;
1841:         }
1842:         if (numFields) {
1843:           PetscSectionRestoreFieldPointSyms(section,f,closureSize,closure,&perms,&flips);
1844:         } else {
1845:           PetscSectionRestorePointSyms(section,closureSize,closure,&perms,&flips);
1846:         }
1847:       }
1848:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1849:     }
1850:     for (row = 0; row < nRows; row++) {
1851:       MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);
1852:     }
1853:     MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1854:     if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix");
1855:     MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1856:     MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1857:     PetscFree(vals);
1858:   }

1860:   /* clean up */
1861:   ISRestoreIndices(anIS,&anchors);
1862:   PetscFree2(perm,iperm);
1863:   PetscFree(pointWork);
1864:   DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
1865:   return(0);
1866: }

1868: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1869:  * a non-conforming mesh.  Local refinement comes later */
1870: PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm)
1871: {
1872:   DM K;
1873:   PetscMPIInt rank;
1874:   PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1875:   PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1876:   PetscInt *Kembedding;
1877:   PetscInt *cellClosure=NULL, nc;
1878:   PetscScalar *newVertexCoords;
1879:   PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1880:   PetscSection parentSection;

1884:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1885:   DMGetDimension(dm,&dim);
1886:   DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);
1887:   DMSetDimension(*ncdm,dim);

1889:   DMPlexGetChart(dm, &pStart, &pEnd);
1890:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);
1891:   DMPlexGetReferenceTree(dm,&K);
1892:   if (!rank) {
1893:     /* compute the new charts */
1894:     PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);
1895:     offset = 0;
1896:     for (d = 0; d <= dim; d++) {
1897:       PetscInt pOldCount, kStart, kEnd, k;

1899:       pNewStart[d] = offset;
1900:       DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);
1901:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1902:       pOldCount = pOldEnd[d] - pOldStart[d];
1903:       /* adding the new points */
1904:       pNewCount[d] = pOldCount + kEnd - kStart;
1905:       if (!d) {
1906:         /* removing the cell */
1907:         pNewCount[d]--;
1908:       }
1909:       for (k = kStart; k < kEnd; k++) {
1910:         PetscInt parent;
1911:         DMPlexGetTreeParent(K,k,&parent,NULL);
1912:         if (parent == k) {
1913:           /* avoid double counting points that won't actually be new */
1914:           pNewCount[d]--;
1915:         }
1916:       }
1917:       pNewEnd[d] = pNewStart[d] + pNewCount[d];
1918:       offset = pNewEnd[d];

1920:     }
1921:     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]);
1922:     /* get the current closure of the cell that we are removing */
1923:     DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);

1925:     PetscMalloc1(pNewEnd[dim],&newConeSizes);
1926:     {
1927:       PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;

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

1932:       for (k = kStart; k < kEnd; k++) {
1933:         perm[k - kStart] = k;
1934:         iperm [k - kStart] = k - kStart;
1935:         preOrient[k - kStart] = 0;
1936:       }

1938:       DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1939:       for (j = 1; j < closureSizeK; j++) {
1940:         PetscInt parentOrientA = closureK[2*j+1];
1941:         PetscInt parentOrientB = cellClosure[2*j+1];
1942:         PetscInt p, q;

1944:         p = closureK[2*j];
1945:         q = cellClosure[2*j];
1946:         for (d = 0; d <= dim; d++) {
1947:           if (q >= pOldStart[d] && q < pOldEnd[d]) {
1948:             Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1949:           }
1950:         }
1951:         if (parentOrientA != parentOrientB) {
1952:           PetscInt numChildren, i;
1953:           const PetscInt *children;

1955:           DMPlexGetTreeChildren(K,p,&numChildren,&children);
1956:           for (i = 0; i < numChildren; i++) {
1957:             PetscInt kPerm, oPerm;

1959:             k    = children[i];
1960:             DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);
1961:             /* perm = what refTree position I'm in */
1962:             perm[kPerm-kStart]      = k;
1963:             /* iperm = who is at this position */
1964:             iperm[k-kStart]         = kPerm-kStart;
1965:             preOrient[kPerm-kStart] = oPerm;
1966:           }
1967:         }
1968:       }
1969:       DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1970:     }
1971:     PetscSectionSetChart(parentSection,0,pNewEnd[dim]);
1972:     offset = 0;
1973:     numNewCones = 0;
1974:     for (d = 0; d <= dim; d++) {
1975:       PetscInt kStart, kEnd, k;
1976:       PetscInt p;
1977:       PetscInt size;

1979:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1980:         /* skip cell 0 */
1981:         if (p == cell) continue;
1982:         /* old cones to new cones */
1983:         DMPlexGetConeSize(dm,p,&size);
1984:         newConeSizes[offset++] = size;
1985:         numNewCones += size;
1986:       }

1988:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1989:       for (k = kStart; k < kEnd; k++) {
1990:         PetscInt kParent;

1992:         DMPlexGetTreeParent(K,k,&kParent,NULL);
1993:         if (kParent != k) {
1994:           Kembedding[k] = offset;
1995:           DMPlexGetConeSize(K,k,&size);
1996:           newConeSizes[offset++] = size;
1997:           numNewCones += size;
1998:           if (kParent != 0) {
1999:             PetscSectionSetDof(parentSection,Kembedding[k],1);
2000:           }
2001:         }
2002:       }
2003:     }

2005:     PetscSectionSetUp(parentSection);
2006:     PetscSectionGetStorageSize(parentSection,&numPointsWithParents);
2007:     PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);
2008:     PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);

2010:     /* fill new cones */
2011:     offset = 0;
2012:     for (d = 0; d <= dim; d++) {
2013:       PetscInt kStart, kEnd, k, l;
2014:       PetscInt p;
2015:       PetscInt size;
2016:       const PetscInt *cone, *orientation;

2018:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
2019:         /* skip cell 0 */
2020:         if (p == cell) continue;
2021:         /* old cones to new cones */
2022:         DMPlexGetConeSize(dm,p,&size);
2023:         DMPlexGetCone(dm,p,&cone);
2024:         DMPlexGetConeOrientation(dm,p,&orientation);
2025:         for (l = 0; l < size; l++) {
2026:           newCones[offset]          = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
2027:           newOrientations[offset++] = orientation[l];
2028:         }
2029:       }

2031:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
2032:       for (k = kStart; k < kEnd; k++) {
2033:         PetscInt kPerm = perm[k], kParent;
2034:         PetscInt preO  = preOrient[k];

2036:         DMPlexGetTreeParent(K,k,&kParent,NULL);
2037:         if (kParent != k) {
2038:           /* embed new cones */
2039:           DMPlexGetConeSize(K,k,&size);
2040:           DMPlexGetCone(K,kPerm,&cone);
2041:           DMPlexGetConeOrientation(K,kPerm,&orientation);
2042:           for (l = 0; l < size; l++) {
2043:             PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size);
2044:             PetscInt newO, lSize, oTrue;

2046:             q                         = iperm[cone[m]];
2047:             newCones[offset]          = Kembedding[q];
2048:             DMPlexGetConeSize(K,q,&lSize);
2049:             oTrue                     = orientation[m];
2050:             oTrue                     = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
2051:             newO                      = DihedralCompose(lSize,oTrue,preOrient[q]);
2052:             newOrientations[offset++] = newO;
2053:           }
2054:           if (kParent != 0) {
2055:             PetscInt newPoint = Kembedding[kParent];
2056:             PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);
2057:             parents[pOffset]  = newPoint;
2058:             childIDs[pOffset] = k;
2059:           }
2060:         }
2061:       }
2062:     }

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

2066:     /* fill coordinates */
2067:     offset = 0;
2068:     {
2069:       PetscInt kStart, kEnd, l;
2070:       PetscSection vSection;
2071:       PetscInt v;
2072:       Vec coords;
2073:       PetscScalar *coordvals;
2074:       PetscInt dof, off;
2075:       PetscReal v0[3], J[9], detJ;

2077: #if defined(PETSC_USE_DEBUG)
2078:       {
2079:         PetscInt k;
2080:         DMPlexGetHeightStratum(K,0,&kStart,&kEnd);
2081:         for (k = kStart; k < kEnd; k++) {
2082:           DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);
2083:           if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k);
2084:         }
2085:       }
2086: #endif
2087:       DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);
2088:       DMGetCoordinateSection(dm,&vSection);
2089:       DMGetCoordinatesLocal(dm,&coords);
2090:       VecGetArray(coords,&coordvals);
2091:       for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {

2093:         PetscSectionGetDof(vSection,v,&dof);
2094:         PetscSectionGetOffset(vSection,v,&off);
2095:         for (l = 0; l < dof; l++) {
2096:           newVertexCoords[offset++] = coordvals[off + l];
2097:         }
2098:       }
2099:       VecRestoreArray(coords,&coordvals);

2101:       DMGetCoordinateSection(K,&vSection);
2102:       DMGetCoordinatesLocal(K,&coords);
2103:       VecGetArray(coords,&coordvals);
2104:       DMPlexGetDepthStratum(K,0,&kStart,&kEnd);
2105:       for (v = kStart; v < kEnd; v++) {
2106:         PetscReal coord[3], newCoord[3];
2107:         PetscInt  vPerm = perm[v];
2108:         PetscInt  kParent;

2110:         DMPlexGetTreeParent(K,v,&kParent,NULL);
2111:         if (kParent != v) {
2112:           /* this is a new vertex */
2113:           PetscSectionGetOffset(vSection,vPerm,&off);
2114:           for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]);
2115:           CoordinatesRefToReal(dim, dim, v0, J, coord, newCoord);
2116:           for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l];
2117:           offset += dim;
2118:         }
2119:       }
2120:       VecRestoreArray(coords,&coordvals);
2121:     }

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

2127:       tmp = pNewCount[d];
2128:       pNewCount[d] = pNewCount[dim - d];
2129:       pNewCount[dim - d] = tmp;
2130:     }

2132:     DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);
2133:     DMPlexSetReferenceTree(*ncdm,K);
2134:     DMPlexSetTree(*ncdm,parentSection,parents,childIDs);

2136:     /* clean up */
2137:     DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
2138:     PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);
2139:     PetscFree(newConeSizes);
2140:     PetscFree2(newCones,newOrientations);
2141:     PetscFree(newVertexCoords);
2142:     PetscFree2(parents,childIDs);
2143:     PetscFree4(Kembedding,perm,iperm,preOrient);
2144:   }
2145:   else {
2146:     PetscInt    p, counts[4];
2147:     PetscInt    *coneSizes, *cones, *orientations;
2148:     Vec         coordVec;
2149:     PetscScalar *coords;

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

2154:       DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
2155:       counts[d] = dEnd - dStart;
2156:     }
2157:     PetscMalloc1(pEnd-pStart,&coneSizes);
2158:     for (p = pStart; p < pEnd; p++) {
2159:       DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);
2160:     }
2161:     DMPlexGetCones(dm, &cones);
2162:     DMPlexGetConeOrientations(dm, &orientations);
2163:     DMGetCoordinatesLocal(dm,&coordVec);
2164:     VecGetArray(coordVec,&coords);

2166:     PetscSectionSetChart(parentSection,pStart,pEnd);
2167:     PetscSectionSetUp(parentSection);
2168:     DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);
2169:     DMPlexSetReferenceTree(*ncdm,K);
2170:     DMPlexSetTree(*ncdm,parentSection,NULL,NULL);
2171:     VecRestoreArray(coordVec,&coords);
2172:   }
2173:   PetscSectionDestroy(&parentSection);

2175:   return(0);
2176: }

2178: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2179: {
2180:   PetscSF           coarseToFineEmbedded;
2181:   PetscSection      globalCoarse, globalFine;
2182:   PetscSection      localCoarse, localFine;
2183:   PetscSection      aSec, cSec;
2184:   PetscSection      rootIndicesSec, rootMatricesSec;
2185:   PetscSection      leafIndicesSec, leafMatricesSec;
2186:   PetscInt          *rootIndices, *leafIndices;
2187:   PetscScalar       *rootMatrices, *leafMatrices;
2188:   IS                aIS;
2189:   const PetscInt    *anchors;
2190:   Mat               cMat;
2191:   PetscInt          numFields, maxFields;
2192:   PetscInt          pStartC, pEndC, pStartF, pEndF, p;
2193:   PetscInt          aStart, aEnd, cStart, cEnd;
2194:   PetscInt          *maxChildIds;
2195:   PetscInt          *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2196:   const PetscInt    ***perms;
2197:   const PetscScalar ***flips;
2198:   PetscErrorCode    ierr;

2201:   DMPlexGetChart(coarse,&pStartC,&pEndC);
2202:   DMPlexGetChart(fine,&pStartF,&pEndF);
2203:   DMGetDefaultGlobalSection(fine,&globalFine);
2204:   { /* winnow fine points that don't have global dofs out of the sf */
2205:     PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2206:     const PetscInt *leaves;

2208:     PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);
2209:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2210:       p = leaves ? leaves[l] : l;
2211:       PetscSectionGetDof(globalFine,p,&dof);
2212:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
2213:       if ((dof - cdof) > 0) {
2214:         numPointsWithDofs++;
2215:       }
2216:     }
2217:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
2218:     for (l = 0, offset = 0; l < nleaves; l++) {
2219:       p = leaves ? leaves[l] : l;
2220:       PetscSectionGetDof(globalFine,p,&dof);
2221:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
2222:       if ((dof - cdof) > 0) {
2223:         pointsWithDofs[offset++] = l;
2224:       }
2225:     }
2226:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
2227:     PetscFree(pointsWithDofs);
2228:   }
2229:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2230:   PetscMalloc1(pEndC-pStartC,&maxChildIds);
2231:   for (p = pStartC; p < pEndC; p++) {
2232:     maxChildIds[p - pStartC] = -2;
2233:   }
2234:   PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);
2235:   PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);

2237:   DMGetDefaultSection(coarse,&localCoarse);
2238:   DMGetDefaultGlobalSection(coarse,&globalCoarse);

2240:   DMPlexGetAnchors(coarse,&aSec,&aIS);
2241:   ISGetIndices(aIS,&anchors);
2242:   PetscSectionGetChart(aSec,&aStart,&aEnd);

2244:   DMGetDefaultConstraints(coarse,&cSec,&cMat);
2245:   PetscSectionGetChart(cSec,&cStart,&cEnd);

2247:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2248:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
2249:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);
2250:   PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);
2251:   PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);
2252:   PetscSectionGetNumFields(localCoarse,&numFields);
2253:   maxFields = PetscMax(1,numFields);
2254:   PetscMalloc7(maxFields+1,&offsets,maxFields+1,&offsetsCopy,maxFields+1,&newOffsets,maxFields+1,&newOffsetsCopy,maxFields+1,&rowOffsets,maxFields+1,&numD,maxFields+1,&numO);
2255:   PetscMalloc2(maxFields+1,(PetscInt****)&perms,maxFields+1,(PetscScalar****)&flips);
2256:   PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **));
2257:   PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **));

2259:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2260:     PetscInt dof, matSize   = 0;
2261:     PetscInt aDof           = 0;
2262:     PetscInt cDof           = 0;
2263:     PetscInt maxChildId     = maxChildIds[p - pStartC];
2264:     PetscInt numRowIndices  = 0;
2265:     PetscInt numColIndices  = 0;
2266:     PetscInt f;

2268:     PetscSectionGetDof(globalCoarse,p,&dof);
2269:     if (dof < 0) {
2270:       dof = -(dof + 1);
2271:     }
2272:     if (p >= aStart && p < aEnd) {
2273:       PetscSectionGetDof(aSec,p,&aDof);
2274:     }
2275:     if (p >= cStart && p < cEnd) {
2276:       PetscSectionGetDof(cSec,p,&cDof);
2277:     }
2278:     for (f = 0; f <= numFields; f++) offsets[f] = 0;
2279:     for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2280:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2281:       PetscInt *closure = NULL, closureSize, cl;

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

2287:         PetscSectionGetDof(localCoarse,c,&clDof);
2288:         numRowIndices += clDof;
2289:         for (f = 0; f < numFields; f++) {
2290:           PetscSectionGetFieldDof(localCoarse,c,f,&clDof);
2291:           offsets[f + 1] += clDof;
2292:         }
2293:       }
2294:       for (f = 0; f < numFields; f++) {
2295:         offsets[f + 1]   += offsets[f];
2296:         newOffsets[f + 1] = offsets[f + 1];
2297:       }
2298:       /* get the number of indices needed and their field offsets */
2299:       DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);
2300:       DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
2301:       if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2302:         numColIndices = numRowIndices;
2303:         matSize = 0;
2304:       }
2305:       else if (numFields) { /* we send one submat for each field: sum their sizes */
2306:         matSize = 0;
2307:         for (f = 0; f < numFields; f++) {
2308:           PetscInt numRow, numCol;

2310:           numRow = offsets[f + 1] - offsets[f];
2311:           numCol = newOffsets[f + 1] - newOffsets[f];
2312:           matSize += numRow * numCol;
2313:         }
2314:       }
2315:       else {
2316:         matSize = numRowIndices * numColIndices;
2317:       }
2318:     } else if (maxChildId == -1) {
2319:       if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2320:         PetscInt aOff, a;

2322:         PetscSectionGetOffset(aSec,p,&aOff);
2323:         for (f = 0; f < numFields; f++) {
2324:           PetscInt fDof;

2326:           PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2327:           offsets[f+1] = fDof;
2328:         }
2329:         for (a = 0; a < aDof; a++) {
2330:           PetscInt anchor = anchors[a + aOff], aLocalDof;

2332:           PetscSectionGetDof(localCoarse,anchor,&aLocalDof);
2333:           numColIndices += aLocalDof;
2334:           for (f = 0; f < numFields; f++) {
2335:             PetscInt fDof;

2337:             PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2338:             newOffsets[f+1] += fDof;
2339:           }
2340:         }
2341:         if (numFields) {
2342:           matSize = 0;
2343:           for (f = 0; f < numFields; f++) {
2344:             matSize += offsets[f+1] * newOffsets[f+1];
2345:           }
2346:         }
2347:         else {
2348:           matSize = numColIndices * dof;
2349:         }
2350:       }
2351:       else { /* no children, and no constraints on dofs: just get the global indices */
2352:         numColIndices = dof;
2353:         matSize       = 0;
2354:       }
2355:     }
2356:     /* we will pack the column indices with the field offsets */
2357:     PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);
2358:     PetscSectionSetDof(rootMatricesSec,p,matSize);
2359:   }
2360:   PetscSectionSetUp(rootIndicesSec);
2361:   PetscSectionSetUp(rootMatricesSec);
2362:   {
2363:     PetscInt numRootIndices, numRootMatrices;

2365:     PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
2366:     PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);
2367:     PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);
2368:     for (p = pStartC; p < pEndC; p++) {
2369:       PetscInt    numRowIndices, numColIndices, matSize, dof;
2370:       PetscInt    pIndOff, pMatOff, f;
2371:       PetscInt    *pInd;
2372:       PetscInt    maxChildId = maxChildIds[p - pStartC];
2373:       PetscScalar *pMat = NULL;

2375:       PetscSectionGetDof(rootIndicesSec,p,&numColIndices);
2376:       if (!numColIndices) {
2377:         continue;
2378:       }
2379:       for (f = 0; f <= numFields; f++) {
2380:         offsets[f]        = 0;
2381:         newOffsets[f]     = 0;
2382:         offsetsCopy[f]    = 0;
2383:         newOffsetsCopy[f] = 0;
2384:       }
2385:       numColIndices -= 2 * numFields;
2386:       PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);
2387:       pInd = &(rootIndices[pIndOff]);
2388:       PetscSectionGetDof(rootMatricesSec,p,&matSize);
2389:       if (matSize) {
2390:         PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);
2391:         pMat = &rootMatrices[pMatOff];
2392:       }
2393:       PetscSectionGetDof(globalCoarse,p,&dof);
2394:       if (dof < 0) {
2395:         dof = -(dof + 1);
2396:       }
2397:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2398:         PetscInt i, j;
2399:         PetscInt numRowIndices = matSize / numColIndices;

2401:         if (!numRowIndices) { /* don't need to calculate the mat, just the indices */
2402:           PetscInt numIndices, *indices;
2403:           DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2404:           if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations");
2405:           for (i = 0; i < numColIndices; i++) {
2406:             pInd[i] = indices[i];
2407:           }
2408:           for (i = 0; i < numFields; i++) {
2409:             pInd[numColIndices + i]             = offsets[i+1];
2410:             pInd[numColIndices + numFields + i] = offsets[i+1];
2411:           }
2412:           DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2413:         }
2414:         else {
2415:           PetscInt closureSize, *closure = NULL, cl;
2416:           PetscScalar *pMatIn, *pMatModified;
2417:           PetscInt numPoints,*points;

2419:           DMGetWorkArray(coarse,numRowIndices * numRowIndices,MPIU_SCALAR,&pMatIn);
2420:           for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2421:             for (j = 0; j < numRowIndices; j++) {
2422:               pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2423:             }
2424:           }
2425:           DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2426:           for (f = 0; f < maxFields; f++) {
2427:             if (numFields) {PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);}
2428:             else           {PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);}
2429:           }
2430:           if (numFields) {
2431:             for (cl = 0; cl < closureSize; cl++) {
2432:               PetscInt c = closure[2 * cl];

2434:               for (f = 0; f < numFields; f++) {
2435:                 PetscInt fDof;

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

2490:               PetscSectionGetOffset(globalCoarse, c, &globalOff);
2491:               DMPlexGetIndicesPoint_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, pInd);
2492:             }
2493:           }
2494:           for (f = 0; f < maxFields; f++) {
2495:             if (numFields) {PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);}
2496:             else           {PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);}
2497:           }
2498:           DMRestoreWorkArray(coarse,numPoints,MPIU_SCALAR,&points);
2499:         }
2500:       }
2501:       else if (matSize) {
2502:         PetscInt cOff;
2503:         PetscInt *rowIndices, *colIndices, a, aDof, aOff;

2505:         numRowIndices = matSize / numColIndices;
2506:         if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs");
2507:         DMGetWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);
2508:         DMGetWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);
2509:         PetscSectionGetOffset(cSec,p,&cOff);
2510:         PetscSectionGetDof(aSec,p,&aDof);
2511:         PetscSectionGetOffset(aSec,p,&aOff);
2512:         if (numFields) {
2513:           for (f = 0; f < numFields; f++) {
2514:             PetscInt fDof;

2516:             PetscSectionGetFieldDof(cSec,p,f,&fDof);
2517:             offsets[f + 1] = fDof;
2518:             for (a = 0; a < aDof; a++) {
2519:               PetscInt anchor = anchors[a + aOff];
2520:               PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2521:               newOffsets[f + 1] += fDof;
2522:             }
2523:           }
2524:           for (f = 0; f < numFields; f++) {
2525:             offsets[f + 1]       += offsets[f];
2526:             offsetsCopy[f + 1]    = offsets[f + 1];
2527:             newOffsets[f + 1]    += newOffsets[f];
2528:             newOffsetsCopy[f + 1] = newOffsets[f + 1];
2529:           }
2530:           DMPlexGetIndicesPointFields_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1,rowIndices);
2531:           for (a = 0; a < aDof; a++) {
2532:             PetscInt anchor = anchors[a + aOff], lOff;
2533:             PetscSectionGetOffset(localCoarse,anchor,&lOff);
2534:             DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1,colIndices);
2535:           }
2536:         }
2537:         else {
2538:           DMPlexGetIndicesPoint_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,rowIndices);
2539:           for (a = 0; a < aDof; a++) {
2540:             PetscInt anchor = anchors[a + aOff], lOff;
2541:             PetscSectionGetOffset(localCoarse,anchor,&lOff);
2542:             DMPlexGetIndicesPoint_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,colIndices);
2543:           }
2544:         }
2545:         if (numFields) {
2546:           PetscInt count, a;

2548:           for (f = 0, count = 0; f < numFields; f++) {
2549:             PetscInt iSize = offsets[f + 1] - offsets[f];
2550:             PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2551:             MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);
2552:             count += iSize * jSize;
2553:             pInd[numColIndices + f]             = offsets[f+1];
2554:             pInd[numColIndices + numFields + f] = newOffsets[f+1];
2555:           }
2556:           for (a = 0; a < aDof; a++) {
2557:             PetscInt anchor = anchors[a + aOff];
2558:             PetscInt gOff;
2559:             PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2560:             DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1,pInd);
2561:           }
2562:         }
2563:         else {
2564:           PetscInt a;
2565:           MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);
2566:           for (a = 0; a < aDof; a++) {
2567:             PetscInt anchor = anchors[a + aOff];
2568:             PetscInt gOff;
2569:             PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2570:             DMPlexGetIndicesPoint_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,pInd);
2571:           }
2572:         }
2573:         DMRestoreWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);
2574:         DMRestoreWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);
2575:       }
2576:       else {
2577:         PetscInt gOff;

2579:         PetscSectionGetOffset(globalCoarse,p,&gOff);
2580:         if (numFields) {
2581:           for (f = 0; f < numFields; f++) {
2582:             PetscInt fDof;
2583:             PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2584:             offsets[f + 1] = fDof + offsets[f];
2585:           }
2586:           for (f = 0; f < numFields; f++) {
2587:             pInd[numColIndices + f]             = offsets[f+1];
2588:             pInd[numColIndices + numFields + f] = offsets[f+1];
2589:           }
2590:           DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);
2591:         }
2592:         else {
2593:           DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);
2594:         }
2595:       }
2596:     }
2597:     PetscFree(maxChildIds);
2598:   }
2599:   {
2600:     PetscSF  indicesSF, matricesSF;
2601:     PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;

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

2641:     PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);
2642:     PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);
2643:     MatGetLayouts(mat,&rowMap,&colMap);
2644:     PetscLayoutSetUp(rowMap);
2645:     PetscLayoutSetUp(colMap);
2646:     PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
2647:     PetscLayoutGetRange(colMap,&colStart,&colEnd);
2648:     PetscSectionGetMaxDof(globalFine,&maxDof);
2649:     PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd);
2650:     DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
2651:     for (p = leafStart; p < leafEnd; p++) {
2652:       PetscInt    gDof, gcDof, gOff;
2653:       PetscInt    numColIndices, pIndOff, *pInd;
2654:       PetscInt    matSize;
2655:       PetscInt    i;

2657:       PetscSectionGetDof(globalFine,p,&gDof);
2658:       PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2659:       if ((gDof - gcDof) <= 0) {
2660:         continue;
2661:       }
2662:       PetscSectionGetOffset(globalFine,p,&gOff);
2663:       if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset");
2664:       if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs");
2665:       PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2666:       PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2667:       numColIndices -= 2 * numFields;
2668:       if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from");
2669:       pInd = &leafIndices[pIndOff];
2670:       offsets[0]        = 0;
2671:       offsetsCopy[0]    = 0;
2672:       newOffsets[0]     = 0;
2673:       newOffsetsCopy[0] = 0;
2674:       if (numFields) {
2675:         PetscInt f;
2676:         for (f = 0; f < numFields; f++) {
2677:           PetscInt rowDof;

2679:           PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2680:           offsets[f + 1]        = offsets[f] + rowDof;
2681:           offsetsCopy[f + 1]    = offsets[f + 1];
2682:           newOffsets[f + 1]     = pInd[numColIndices + numFields + f];
2683:           numD[f] = 0;
2684:           numO[f] = 0;
2685:         }
2686:         DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);
2687:         for (f = 0; f < numFields; f++) {
2688:           PetscInt colOffset    = newOffsets[f];
2689:           PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];

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

2694:             if (gInd >= colStart && gInd < colEnd) {
2695:               numD[f]++;
2696:             }
2697:             else if (gInd >= 0) { /* negative means non-entry */
2698:               numO[f]++;
2699:             }
2700:           }
2701:         }
2702:       }
2703:       else {
2704:         DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);
2705:         numD[0] = 0;
2706:         numO[0] = 0;
2707:         for (i = 0; i < numColIndices; i++) {
2708:           PetscInt gInd = pInd[i];

2710:           if (gInd >= colStart && gInd < colEnd) {
2711:             numD[0]++;
2712:           }
2713:           else if (gInd >= 0) { /* negative means non-entry */
2714:             numO[0]++;
2715:           }
2716:         }
2717:       }
2718:       PetscSectionGetDof(leafMatricesSec,p,&matSize);
2719:       if (!matSize) { /* incoming matrix is identity */
2720:         PetscInt childId;

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

2823:     DMPlexGetReferenceTree(fine,&refTree);
2824:     DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
2825:     DMGetDefaultConstraints(refTree,&refConSec,NULL);
2826:     DMPlexGetAnchors(refTree,&refAnSec,NULL);
2827:     PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
2828:     PetscSectionGetMaxDof(refConSec,&maxConDof);
2829:     PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);
2830:     PetscMalloc1(maxConDof*maxColumns,&pointWork);
2831:     for (p = leafStart; p < leafEnd; p++) {
2832:       PetscInt    gDof, gcDof, gOff;
2833:       PetscInt    numColIndices, pIndOff, *pInd;
2834:       PetscInt    matSize;
2835:       PetscInt    childId;


2838:       PetscSectionGetDof(globalFine,p,&gDof);
2839:       PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2840:       if ((gDof - gcDof) <= 0) {
2841:         continue;
2842:       }
2843:       childId = childIds[p-pStartF];
2844:       PetscSectionGetOffset(globalFine,p,&gOff);
2845:       PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2846:       PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2847:       numColIndices -= 2 * numFields;
2848:       pInd = &leafIndices[pIndOff];
2849:       offsets[0]        = 0;
2850:       offsetsCopy[0]    = 0;
2851:       newOffsets[0]     = 0;
2852:       newOffsetsCopy[0] = 0;
2853:       rowOffsets[0]     = 0;
2854:       if (numFields) {
2855:         PetscInt f;
2856:         for (f = 0; f < numFields; f++) {
2857:           PetscInt rowDof;

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

2907:         PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);
2908:         pMat = &leafMatrices[pMatOff];
2909:         if (childId < 0) { /* copy the incoming matrix */
2910:           if (numFields) {
2911:             PetscInt f, count;
2912:             for (f = 0, count = 0; f < numFields; f++) {
2913:               PetscInt numRows = offsets[f+1]-offsets[f];
2914:               PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2915:               PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2916:               PetscScalar *inMat = &pMat[count];

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

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

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

3029:     if (canonical) {
3030:       PetscInt pCanonical;
3031:       DMLabelGetValue(canonical,p,&pCanonical);
3032:       if (p != pCanonical) continue;
3033:     }
3034:     DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3035:     if (!numChildren) continue;
3036:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
3037:       PetscInt child = children[i];
3038:       PetscInt dof;

3040:       PetscSectionGetDof(section,child,&dof);
3041:       numChildDof += dof;
3042:     }
3043:     PetscSectionGetDof(section,p,&numSelfDof);
3044:     if (!numChildDof || !numSelfDof) continue;
3045:     for (f = 0; f < numFields; f++) {
3046:       PetscInt selfOff;

3048:       if (numSecFields) { /* count the dofs for just this field */
3049:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
3050:           PetscInt child = children[i];
3051:           PetscInt dof;

3053:           PetscSectionGetFieldDof(section,child,f,&dof);
3054:           numChildDof += dof;
3055:         }
3056:         PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3057:         PetscSectionGetFieldOffset(section,p,f,&selfOff);
3058:       }
3059:       else {
3060:         PetscSectionGetOffset(section,p,&selfOff);
3061:       }
3062:       for (i = 0; i < numSelfDof; i++) {
3063:         nnz[selfOff + i] = numChildDof;
3064:       }
3065:     }
3066:   }
3067:   MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);
3068:   PetscFree(nnz);
3069:   /* Setp 2: compute entries */
3070:   for (p = pStart; p < pEnd; p++) {
3071:     const PetscInt *children;
3072:     PetscInt numChildren;
3073:     PetscInt i, numChildDof, numSelfDof;

3075:     /* same conditions about when entries occur */
3076:     if (canonical) {
3077:       PetscInt pCanonical;
3078:       DMLabelGetValue(canonical,p,&pCanonical);
3079:       if (p != pCanonical) continue;
3080:     }
3081:     DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3082:     if (!numChildren) continue;
3083:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
3084:       PetscInt child = children[i];
3085:       PetscInt dof;

3087:       PetscSectionGetDof(section,child,&dof);
3088:       numChildDof += dof;
3089:     }
3090:     PetscSectionGetDof(section,p,&numSelfDof);
3091:     if (!numChildDof || !numSelfDof) continue;

3093:     for (f = 0; f < numFields; f++) {
3094:       PetscInt       selfOff, Nc, parentCell;
3095:       PetscInt       cellShapeOff;
3096:       PetscObject    disc;
3097:       PetscDualSpace dsp;
3098:       PetscClassId   classId;
3099:       PetscScalar    *pointMat;
3100:       PetscInt       *matRows, *matCols;
3101:       PetscInt       pO = PETSC_MIN_INT;
3102:       const PetscInt *depthNumDof;

3104:       if (numSecFields) {
3105:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
3106:           PetscInt child = children[i];
3107:           PetscInt dof;

3109:           PetscSectionGetFieldDof(section,child,f,&dof);
3110:           numChildDof += dof;
3111:         }
3112:         PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3113:         PetscSectionGetFieldOffset(section,p,f,&selfOff);
3114:       }
3115:       else {
3116:         PetscSectionGetOffset(section,p,&selfOff);
3117:       }

3119:       /* find a cell whose closure contains p */
3120:       if (p >= cStart && p < cEnd) {
3121:         parentCell = p;
3122:       }
3123:       else {
3124:         PetscInt *star = NULL;
3125:         PetscInt numStar;

3127:         parentCell = -1;
3128:         DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3129:         for (i = numStar - 1; i >= 0; i--) {
3130:           PetscInt c = star[2 * i];

3132:           if (c >= cStart && c < cEnd) {
3133:             parentCell = c;
3134:             break;
3135:           }
3136:         }
3137:         DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3138:       }
3139:       /* determine the offset of p's shape functions withing parentCell's shape functions */
3140:       PetscDSGetDiscretization(ds,f,&disc);
3141:       PetscObjectGetClassId(disc,&classId);
3142:       if (classId == PETSCFE_CLASSID) {
3143:         PetscFEGetDualSpace((PetscFE)disc,&dsp);
3144:       }
3145:       else if (classId == PETSCFV_CLASSID) {
3146:         PetscFVGetDualSpace((PetscFV)disc,&dsp);
3147:       }
3148:       else {
3149:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object");
3150:       }
3151:       PetscDualSpaceGetNumDof(dsp,&depthNumDof);
3152:       PetscDualSpaceGetNumComponents(dsp,&Nc);
3153:       {
3154:         PetscInt *closure = NULL;
3155:         PetscInt numClosure;

3157:         DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);
3158:         for (i = 0, cellShapeOff = 0; i < numClosure; i++) {
3159:           PetscInt point = closure[2 * i], pointDepth;

3161:           pO = closure[2 * i + 1];
3162:           if (point == p) break;
3163:           DMLabelGetValue(depth,point,&pointDepth);
3164:           cellShapeOff += depthNumDof[pointDepth];
3165:         }
3166:         DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);
3167:       }

3169:       DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);
3170:       DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);
3171:       matCols = matRows + numSelfDof;
3172:       for (i = 0; i < numSelfDof; i++) {
3173:         matRows[i] = selfOff + i;
3174:       }
3175:       for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
3176:       {
3177:         PetscInt colOff = 0;

3179:         for (i = 0; i < numChildren; i++) {
3180:           PetscInt child = children[i];
3181:           PetscInt dof, off, j;

3183:           if (numSecFields) {
3184:             PetscSectionGetFieldDof(cSection,child,f,&dof);
3185:             PetscSectionGetFieldOffset(cSection,child,f,&off);
3186:           }
3187:           else {
3188:             PetscSectionGetDof(cSection,child,&dof);
3189:             PetscSectionGetOffset(cSection,child,&off);
3190:           }

3192:           for (j = 0; j < dof; j++) {
3193:             matCols[colOff++] = off + j;
3194:           }
3195:         }
3196:       }
3197:       if (classId == PETSCFE_CLASSID) {
3198:         PetscFE        fe = (PetscFE) disc;
3199:         PetscInt       fSize;

3201:         PetscFEGetDualSpace(fe,&dsp);
3202:         PetscDualSpaceGetDimension(dsp,&fSize);
3203:         for (i = 0; i < numSelfDof; i++) { /* for every shape function */
3204:           PetscQuadrature q;
3205:           PetscInt        dim, thisNc, numPoints, j, k;
3206:           const PetscReal *points;
3207:           const PetscReal *weights;
3208:           PetscInt        *closure = NULL;
3209:           PetscInt        numClosure;
3210:           PetscInt        parentCellShapeDof = cellShapeOff + (pO < 0 ? (numSelfDof - 1 - i) : i);
3211:           PetscReal       *Bparent;

3213:           PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);
3214:           PetscQuadratureGetData(q,&dim,&thisNc,&numPoints,&points,&weights);
3215:           if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc);
3216:           PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3217:           for (j = 0; j < numPoints; j++) {
3218:             PetscInt          childCell = -1;
3219:             PetscReal         *parentValAtPoint;
3220:             const PetscReal   *pointReal = &points[dim * j];
3221:             const PetscScalar *point;
3222:             PetscReal         *Bchild;
3223:             PetscInt          childCellShapeOff, pointMatOff;
3224: #if defined(PETSC_USE_COMPLEX)
3225:             PetscInt          d;

3227:             for (d = 0; d < dim; d++) {
3228:               pointScalar[d] = points[dim * j + d];
3229:             }
3230:             point = pointScalar;
3231: #else
3232:             point = pointReal;
3233: #endif

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

3237:             for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3238:               PetscInt child = children[k];
3239:               PetscInt *star = NULL;
3240:               PetscInt numStar, s;

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

3246:                 if (c < cStart || c >= cEnd) continue;
3247:                 DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);
3248:                 if (childCell >= 0) break;
3249:               }
3250:               DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);
3251:               if (childCell >= 0) break;
3252:             }
3253:             if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point");
3254:             DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);
3255:             DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);
3256:             CoordinatesRefToReal(dim, dim, v0parent, Jparent, pointReal, vtmp);
3257:             CoordinatesRealToRef(dim, dim, v0, invJ, vtmp, pointRef);

3259:             PetscFEGetTabulation(fe,1,pointRef,&Bchild,NULL,NULL);
3260:             DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3261:             for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3262:               PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT;
3263:               PetscInt l;

3265:               DMLabelGetValue(depth,child,&childDepth);
3266:               childDof = depthNumDof[childDepth];
3267:               for (l = 0, childCellShapeOff = 0; l < numClosure; l++) {
3268:                 PetscInt point = closure[2 * l];
3269:                 PetscInt pointDepth;

3271:                 childO = closure[2 * l + 1];
3272:                 if (point == child) break;
3273:                 DMLabelGetValue(depth,point,&pointDepth);
3274:                 childCellShapeOff += depthNumDof[pointDepth];
3275:               }
3276:               if (l == numClosure) {
3277:                 pointMatOff += childDof;
3278:                 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3279:               }
3280:               for (l = 0; l < childDof; l++) {
3281:                 PetscInt    childCellDof = childCellShapeOff + (childO ? (childDof - 1 - l) : l), m;
3282:                 PetscReal   *childValAtPoint;
3283:                 PetscReal   val = 0.;

3285:                 childValAtPoint = &Bchild[childCellDof * Nc];
3286:                 for (m = 0; m < Nc; m++) {
3287:                   val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];
3288:                 }

3290:                 pointMat[i * numChildDof + pointMatOff + l] += val;
3291:               }
3292:               pointMatOff += childDof;
3293:             }
3294:             DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3295:             PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);
3296:           }
3297:           PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);
3298:         }
3299:       }
3300:       else { /* just the volume-weighted averages of the children */
3301:         PetscReal parentVol;
3302:         PetscInt  childCell;

3304:         DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);
3305:         for (i = 0, childCell = 0; i < numChildren; i++) {
3306:           PetscInt  child = children[i], j;
3307:           PetscReal childVol;

3309:           if (child < cStart || child >= cEnd) continue;
3310:           DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);
3311:           for (j = 0; j < Nc; j++) {
3312:             pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3313:           }
3314:           childCell++;
3315:         }
3316:       }
3317:       /* Insert pointMat into mat */
3318:       MatSetValues(mat,numSelfDof,matRows,numChildDof,matCols,pointMat,INSERT_VALUES);
3319:       DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);
3320:       DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);
3321:     }
3322:   }
3323:   PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);
3324:   PetscFree2(pointScalar,pointRef);
3325:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3326:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3327:   *inj = mat;
3328:   return(0);
3329: }

3331: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3332: {
3333:   PetscDS        ds;
3334:   PetscInt       numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3335:   PetscScalar    ***refPointFieldMats;
3336:   PetscSection   refConSec, refSection;

3340:   DMGetDS(refTree,&ds);
3341:   PetscDSGetNumFields(ds,&numFields);
3342:   DMGetDefaultConstraints(refTree,&refConSec,NULL);
3343:   DMGetDefaultSection(refTree,&refSection);
3344:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3345:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
3346:   PetscSectionGetMaxDof(refConSec,&maxDof);
3347:   PetscMalloc1(maxDof,&rows);
3348:   PetscMalloc1(maxDof*maxDof,&cols);
3349:   for (p = pRefStart; p < pRefEnd; p++) {
3350:     PetscInt parent, pDof, parentDof;

3352:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
3353:     PetscSectionGetDof(refConSec,p,&pDof);
3354:     PetscSectionGetDof(refSection,parent,&parentDof);
3355:     if (!pDof || !parentDof || parent == p) continue;

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

3361:       if (numFields > 1) {
3362:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3363:         PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
3364:       }
3365:       else {
3366:         PetscSectionGetDof(refConSec,p,&cDof);
3367:         PetscSectionGetOffset(refConSec,p,&cOff);
3368:       }

3370:       for (r = 0; r < cDof; r++) {
3371:         rows[r] = cOff + r;
3372:       }
3373:       numCols = 0;
3374:       {
3375:         PetscInt aDof, aOff, j;

3377:         if (numFields > 1) {
3378:           PetscSectionGetFieldDof(refSection,parent,f,&aDof);
3379:           PetscSectionGetFieldOffset(refSection,parent,f,&aOff);
3380:         }
3381:         else {
3382:           PetscSectionGetDof(refSection,parent,&aDof);
3383:           PetscSectionGetOffset(refSection,parent,&aOff);
3384:         }

3386:         for (j = 0; j < aDof; j++) {
3387:           cols[numCols++] = aOff + j;
3388:         }
3389:       }
3390:       PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
3391:       /* transpose of constraint matrix */
3392:       MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);
3393:     }
3394:   }
3395:   *childrenMats = refPointFieldMats;
3396:   PetscFree(rows);
3397:   PetscFree(cols);
3398:   return(0);
3399: }

3401: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3402: {
3403:   PetscDS        ds;
3404:   PetscScalar    ***refPointFieldMats;
3405:   PetscInt       numFields, pRefStart, pRefEnd, p, f;
3406:   PetscSection   refConSec, refSection;

3410:   refPointFieldMats = *childrenMats;
3411:   *childrenMats = NULL;
3412:   DMGetDS(refTree,&ds);
3413:   DMGetDefaultSection(refTree,&refSection);
3414:   PetscDSGetNumFields(ds,&numFields);
3415:   DMGetDefaultConstraints(refTree,&refConSec,NULL);
3416:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3417:   for (p = pRefStart; p < pRefEnd; p++) {
3418:     PetscInt parent, pDof, parentDof;

3420:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
3421:     PetscSectionGetDof(refConSec,p,&pDof);
3422:     PetscSectionGetDof(refSection,parent,&parentDof);
3423:     if (!pDof || !parentDof || parent == p) continue;

3425:     for (f = 0; f < numFields; f++) {
3426:       PetscInt cDof;

3428:       if (numFields > 1) {
3429:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3430:       }
3431:       else {
3432:         PetscSectionGetDof(refConSec,p,&cDof);
3433:       }

3435:       PetscFree(refPointFieldMats[p - pRefStart][f]);
3436:     }
3437:     PetscFree(refPointFieldMats[p - pRefStart]);
3438:   }
3439:   PetscFree(refPointFieldMats);
3440:   return(0);
3441: }

3443: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef)
3444: {
3445:   Mat            cMatRef;
3446:   PetscObject    injRefObj;

3450:   DMGetDefaultConstraints(refTree,NULL,&cMatRef);
3451:   PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);
3452:   *injRef = (Mat) injRefObj;
3453:   if (!*injRef) {
3454:     DMPlexComputeInjectorReferenceTree(refTree,injRef);
3455:     PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);
3456:     /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3457:     PetscObjectDereference((PetscObject)*injRef);
3458:   }
3459:   return(0);
3460: }

3462: 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)
3463: {
3464:   PetscInt       pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3465:   PetscSection   globalCoarse, globalFine;
3466:   PetscSection   localCoarse, localFine, leafIndicesSec;
3467:   PetscSection   multiRootSec, rootIndicesSec;
3468:   PetscInt       *leafInds, *rootInds = NULL;
3469:   const PetscInt *rootDegrees;
3470:   PetscScalar    *leafVals = NULL, *rootVals = NULL;
3471:   PetscSF        coarseToFineEmbedded;

3475:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3476:   DMPlexGetChart(fine,&pStartF,&pEndF);
3477:   DMGetDefaultSection(fine,&localFine);
3478:   DMGetDefaultGlobalSection(fine,&globalFine);
3479:   PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
3480:   PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);
3481:   PetscSectionGetMaxDof(localFine,&maxDof);
3482:   { /* winnow fine points that don't have global dofs out of the sf */
3483:     PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3484:     const PetscInt *leaves;

3486:     PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);
3487:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3488:       p    = leaves ? leaves[l] : l;
3489:       PetscSectionGetDof(globalFine,p,&dof);
3490:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3491:       if ((dof - cdof) > 0) {
3492:         numPointsWithDofs++;

3494:         PetscSectionGetDof(localFine,p,&dof);
3495:         PetscSectionSetDof(leafIndicesSec,p,dof + 1);
3496:       }
3497:     }
3498:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3499:     PetscSectionSetUp(leafIndicesSec);
3500:     PetscSectionGetStorageSize(leafIndicesSec,&numIndices);
3501:     PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);
3502:     if (gatheredValues)  {PetscMalloc1(numIndices,&leafVals);}
3503:     for (l = 0, offset = 0; l < nleaves; l++) {
3504:       p    = leaves ? leaves[l] : l;
3505:       PetscSectionGetDof(globalFine,p,&dof);
3506:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3507:       if ((dof - cdof) > 0) {
3508:         PetscInt    off, gOff;
3509:         PetscInt    *pInd;
3510:         PetscScalar *pVal = NULL;

3512:         pointsWithDofs[offset++] = l;

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

3516:         pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3517:         if (gatheredValues) {
3518:           PetscInt i;

3520:           pVal = &leafVals[off + 1];
3521:           for (i = 0; i < dof; i++) pVal[i] = 0.;
3522:         }
3523:         PetscSectionGetOffset(globalFine,p,&gOff);

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

3529:           for (f = 0; f < numFields; f++) {
3530:             PetscInt fDof;
3531:             PetscSectionGetFieldDof(localFine,p,f,&fDof);
3532:             offsets[f + 1] = fDof + offsets[f];
3533:           }
3534:           DMPlexGetIndicesPointFields_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);
3535:         }
3536:         else {
3537:           DMPlexGetIndicesPoint_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);
3538:         }
3539:         if (gatheredValues) {VecGetValues(fineVec,dof,pInd,pVal);}
3540:       }
3541:     }
3542:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
3543:     PetscFree(pointsWithDofs);
3544:   }

3546:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3547:   DMGetDefaultSection(coarse,&localCoarse);
3548:   DMGetDefaultGlobalSection(coarse,&globalCoarse);

3550:   { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3551:     MPI_Datatype threeInt;
3552:     PetscMPIInt  rank;
3553:     PetscInt     (*parentNodeAndIdCoarse)[3];
3554:     PetscInt     (*parentNodeAndIdFine)[3];
3555:     PetscInt     p, nleaves, nleavesToParents;
3556:     PetscSF      pointSF, sfToParents;
3557:     const PetscInt *ilocal;
3558:     const PetscSFNode *iremote;
3559:     PetscSFNode  *iremoteToParents;
3560:     PetscInt     *ilocalToParents;

3562:     MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
3563:     MPI_Type_contiguous(3,MPIU_INT,&threeInt);
3564:     MPI_Type_commit(&threeInt);
3565:     PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);
3566:     DMGetPointSF(coarse,&pointSF);
3567:     PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);
3568:     for (p = pStartC; p < pEndC; p++) {
3569:       PetscInt parent, childId;
3570:       DMPlexGetTreeParent(coarse,p,&parent,&childId);
3571:       parentNodeAndIdCoarse[p - pStartC][0] = rank;
3572:       parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3573:       parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3574:       if (nleaves > 0) {
3575:         PetscInt leaf = -1;

3577:         if (ilocal) {
3578:           PetscFindInt(parent,nleaves,ilocal,&leaf);
3579:         }
3580:         else {
3581:           leaf = p - pStartC;
3582:         }
3583:         if (leaf >= 0) {
3584:           parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3585:           parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3586:         }
3587:       }
3588:     }
3589:     for (p = pStartF; p < pEndF; p++) {
3590:       parentNodeAndIdFine[p - pStartF][0] = -1;
3591:       parentNodeAndIdFine[p - pStartF][1] = -1;
3592:       parentNodeAndIdFine[p - pStartF][2] = -1;
3593:     }
3594:     PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3595:     PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3596:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3597:       PetscInt dof;

3599:       PetscSectionGetDof(leafIndicesSec,p,&dof);
3600:       if (dof) {
3601:         PetscInt off;

3603:         PetscSectionGetOffset(leafIndicesSec,p,&off);
3604:         if (gatheredIndices) {
3605:           leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3606:         } else if (gatheredValues) {
3607:           leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3608:         }
3609:       }
3610:       if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3611:         nleavesToParents++;
3612:       }
3613:     }
3614:     PetscMalloc1(nleavesToParents,&ilocalToParents);
3615:     PetscMalloc1(nleavesToParents,&iremoteToParents);
3616:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3617:       if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3618:         ilocalToParents[nleavesToParents] = p - pStartF;
3619:         iremoteToParents[nleavesToParents].rank  = parentNodeAndIdFine[p-pStartF][0];
3620:         iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1];
3621:         nleavesToParents++;
3622:       }
3623:     }
3624:     PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);
3625:     PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);
3626:     PetscSFDestroy(&coarseToFineEmbedded);

3628:     coarseToFineEmbedded = sfToParents;

3630:     PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);
3631:     MPI_Type_free(&threeInt);
3632:   }

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

3638:     for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3639:       PetscSectionGetDof(globalCoarse,p,&dof);
3640:       PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3641:       if ((dof - cdof) > 0) {
3642:         numPointsWithDofs++;
3643:       }
3644:     }
3645:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3646:     for (p = pStartC, offset = 0; p < pEndC; p++) {
3647:       PetscSectionGetDof(globalCoarse,p,&dof);
3648:       PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3649:       if ((dof - cdof) > 0) {
3650:         pointsWithDofs[offset++] = p - pStartC;
3651:       }
3652:     }
3653:     PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);
3654:     PetscSFDestroy(&coarseToFineEmbedded);
3655:     PetscFree(pointsWithDofs);
3656:     coarseToFineEmbedded = sfDofsOnly;
3657:   }

3659:   /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3660:   PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);
3661:   PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);
3662:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);
3663:   PetscSectionSetChart(multiRootSec,pStartC,pEndC);
3664:   for (p = pStartC; p < pEndC; p++) {
3665:     PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);
3666:   }
3667:   PetscSectionSetUp(multiRootSec);
3668:   PetscSectionGetStorageSize(multiRootSec,&numMulti);
3669:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
3670:   { /* distribute the leaf section */
3671:     PetscSF multi, multiInv, indicesSF;
3672:     PetscInt *remoteOffsets, numRootIndices;

3674:     PetscSFGetMultiSF(coarseToFineEmbedded,&multi);
3675:     PetscSFCreateInverseSF(multi,&multiInv);
3676:     PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);
3677:     PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);
3678:     PetscFree(remoteOffsets);
3679:     PetscSFDestroy(&multiInv);
3680:     PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
3681:     if (gatheredIndices) {
3682:       PetscMalloc1(numRootIndices,&rootInds);
3683:       PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds);
3684:       PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds);
3685:     }
3686:     if (gatheredValues) {
3687:       PetscMalloc1(numRootIndices,&rootVals);
3688:       PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals);
3689:       PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals);
3690:     }
3691:     PetscSFDestroy(&indicesSF);
3692:   }
3693:   PetscSectionDestroy(&leafIndicesSec);
3694:   PetscFree(leafInds);
3695:   PetscFree(leafVals);
3696:   PetscSFDestroy(&coarseToFineEmbedded);
3697:   *rootMultiSec = multiRootSec;
3698:   *multiLeafSec = rootIndicesSec;
3699:   if (gatheredIndices) *gatheredIndices = rootInds;
3700:   if (gatheredValues)  *gatheredValues  = rootVals;
3701:   return(0);
3702: }

3704: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3705: {
3706:   DM             refTree;
3707:   PetscSection   multiRootSec, rootIndicesSec;
3708:   PetscSection   globalCoarse, globalFine;
3709:   PetscSection   localCoarse, localFine;
3710:   PetscSection   cSecRef;
3711:   PetscInt       *rootIndices, *parentIndices, pRefStart, pRefEnd;
3712:   Mat            injRef;
3713:   PetscInt       numFields, maxDof;
3714:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
3715:   PetscInt       *offsets, *offsetsCopy, *rowOffsets;
3716:   PetscLayout    rowMap, colMap;
3717:   PetscInt       rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3718:   PetscScalar    ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */


3723:   /* get the templates for the fine-to-coarse injection from the reference tree */
3724:   DMPlexGetReferenceTree(coarse,&refTree);
3725:   DMGetDefaultConstraints(refTree,&cSecRef,NULL);
3726:   PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
3727:   DMPlexReferenceTreeGetInjector(refTree,&injRef);

3729:   DMPlexGetChart(fine,&pStartF,&pEndF);
3730:   DMGetDefaultSection(fine,&localFine);
3731:   DMGetDefaultGlobalSection(fine,&globalFine);
3732:   PetscSectionGetNumFields(localFine,&numFields);
3733:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3734:   DMGetDefaultSection(coarse,&localCoarse);
3735:   DMGetDefaultGlobalSection(coarse,&globalCoarse);
3736:   PetscSectionGetMaxDof(localCoarse,&maxDof);
3737:   {
3738:     PetscInt maxFields = PetscMax(1,numFields) + 1;
3739:     PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
3740:   }

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

3744:   PetscMalloc1(maxDof,&parentIndices);

3746:   /* count indices */
3747:   MatGetLayouts(mat,&rowMap,&colMap);
3748:   PetscLayoutSetUp(rowMap);
3749:   PetscLayoutSetUp(colMap);
3750:   PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
3751:   PetscLayoutGetRange(colMap,&colStart,&colEnd);
3752:   PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);
3753:   for (p = pStartC; p < pEndC; p++) {
3754:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3756:     PetscSectionGetDof(globalCoarse,p,&dof);
3757:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3758:     if ((dof - cdof) <= 0) continue;
3759:     PetscSectionGetOffset(globalCoarse,p,&gOff);

3761:     rowOffsets[0] = 0;
3762:     offsetsCopy[0] = 0;
3763:     if (numFields) {
3764:       PetscInt f;

3766:       for (f = 0; f < numFields; f++) {
3767:         PetscInt fDof;
3768:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3769:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3770:       }
3771:       DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);
3772:     }
3773:     else {
3774:       DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);
3775:       rowOffsets[1] = offsetsCopy[0];
3776:     }

3778:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
3779:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
3780:     leafEnd = leafStart + numLeaves;
3781:     for (l = leafStart; l < leafEnd; l++) {
3782:       PetscInt numIndices, childId, offset;
3783:       const PetscInt *childIndices;

3785:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3786:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
3787:       childId = rootIndices[offset++];
3788:       childIndices = &rootIndices[offset];
3789:       numIndices--;

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

3794:         for (i = 0; i < numIndices; i++) {
3795:           PetscInt colIndex = childIndices[i];
3796:           PetscInt rowIndex = parentIndices[i];
3797:           if (rowIndex < 0) continue;
3798:           if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse");
3799:           if (colIndex >= colStart && colIndex < colEnd) {
3800:             nnzD[rowIndex - rowStart] = 1;
3801:           }
3802:           else {
3803:             nnzO[rowIndex - rowStart] = 1;
3804:           }
3805:         }
3806:       }
3807:       else {
3808:         PetscInt parentId, f, lim;

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

3812:         lim = PetscMax(1,numFields);
3813:         offsets[0] = 0;
3814:         if (numFields) {
3815:           PetscInt f;

3817:           for (f = 0; f < numFields; f++) {
3818:             PetscInt fDof;
3819:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

3821:             offsets[f + 1] = fDof + offsets[f];
3822:           }
3823:         }
3824:         else {
3825:           PetscInt cDof;

3827:           PetscSectionGetDof(cSecRef,childId,&cDof);
3828:           offsets[1] = cDof;
3829:         }
3830:         for (f = 0; f < lim; f++) {
3831:           PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3832:           PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3833:           PetscInt i, numD = 0, numO = 0;

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

3838:             if (colIndex < 0) continue;
3839:             if (colIndex >= colStart && colIndex < colEnd) {
3840:               numD++;
3841:             }
3842:             else {
3843:               numO++;
3844:             }
3845:           }
3846:           for (i = parentStart; i < parentEnd; i++) {
3847:             PetscInt rowIndex = parentIndices[i];

3849:             if (rowIndex < 0) continue;
3850:             nnzD[rowIndex - rowStart] += numD;
3851:             nnzO[rowIndex - rowStart] += numO;
3852:           }
3853:         }
3854:       }
3855:     }
3856:   }
3857:   /* preallocate */
3858:   MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);
3859:   PetscFree2(nnzD,nnzO);
3860:   /* insert values */
3861:   DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3862:   for (p = pStartC; p < pEndC; p++) {
3863:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3865:     PetscSectionGetDof(globalCoarse,p,&dof);
3866:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3867:     if ((dof - cdof) <= 0) continue;
3868:     PetscSectionGetOffset(globalCoarse,p,&gOff);

3870:     rowOffsets[0] = 0;
3871:     offsetsCopy[0] = 0;
3872:     if (numFields) {
3873:       PetscInt f;

3875:       for (f = 0; f < numFields; f++) {
3876:         PetscInt fDof;
3877:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3878:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3879:       }
3880:       DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);
3881:     }
3882:     else {
3883:       DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);
3884:       rowOffsets[1] = offsetsCopy[0];
3885:     }

3887:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
3888:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
3889:     leafEnd = leafStart + numLeaves;
3890:     for (l = leafStart; l < leafEnd; l++) {
3891:       PetscInt numIndices, childId, offset;
3892:       const PetscInt *childIndices;

3894:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3895:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
3896:       childId = rootIndices[offset++];
3897:       childIndices = &rootIndices[offset];
3898:       numIndices--;

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

3903:         for (i = 0; i < numIndices; i++) {
3904:           MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);
3905:         }
3906:       }
3907:       else {
3908:         PetscInt parentId, f, lim;

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

3912:         lim = PetscMax(1,numFields);
3913:         offsets[0] = 0;
3914:         if (numFields) {
3915:           PetscInt f;

3917:           for (f = 0; f < numFields; f++) {
3918:             PetscInt fDof;
3919:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

3921:             offsets[f + 1] = fDof + offsets[f];
3922:           }
3923:         }
3924:         else {
3925:           PetscInt cDof;

3927:           PetscSectionGetDof(cSecRef,childId,&cDof);
3928:           offsets[1] = cDof;
3929:         }
3930:         for (f = 0; f < lim; f++) {
3931:           PetscScalar    *childMat   = &childrenMats[childId - pRefStart][f][0];
3932:           PetscInt       *rowIndices = &parentIndices[rowOffsets[f]];
3933:           const PetscInt *colIndices = &childIndices[offsets[f]];

3935:           MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);
3936:         }
3937:       }
3938:     }
3939:   }
3940:   PetscSectionDestroy(&multiRootSec);
3941:   PetscSectionDestroy(&rootIndicesSec);
3942:   PetscFree(parentIndices);
3943:   DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3944:   PetscFree(rootIndices);
3945:   PetscFree3(offsets,offsetsCopy,rowOffsets);

3947:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3948:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3949:   return(0);
3950: }

3952: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3953: {
3954:   PetscDS           ds;
3955:   PetscSF           coarseToFineEmbedded;
3956:   PetscSection      globalCoarse, globalFine;
3957:   PetscSection      localCoarse, localFine;
3958:   PetscSection      aSec, cSec;
3959:   PetscSection      rootValuesSec;
3960:   PetscSection      leafValuesSec;
3961:   PetscScalar       *rootValues, *leafValues;
3962:   IS                aIS;
3963:   const PetscInt    *anchors;
3964:   Mat               cMat;
3965:   PetscInt          numFields;
3966:   PetscInt          pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd, cellEndInterior;
3967:   PetscInt          aStart, aEnd, cStart, cEnd;
3968:   PetscInt          *maxChildIds;
3969:   PetscInt          *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3970:   PetscFV           fv = NULL;
3971:   PetscInt          dim, numFVcomps = -1, fvField = -1;
3972:   DM                cellDM = NULL, gradDM = NULL;
3973:   const PetscScalar *cellGeomArray = NULL;
3974:   const PetscScalar *gradArray = NULL;
3975:   PetscErrorCode    ierr;

3978:   VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
3979:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3980:   DMPlexGetHeightStratum(coarse,0,&cellStart,&cellEnd);
3981:   DMPlexGetHybridBounds(coarse,&cellEndInterior,NULL,NULL,NULL);
3982:   cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior;
3983:   DMPlexGetChart(fine,&pStartF,&pEndF);
3984:   DMGetDefaultGlobalSection(fine,&globalFine);
3985:   DMGetCoordinateDim(coarse,&dim);
3986:   { /* winnow fine points that don't have global dofs out of the sf */
3987:     PetscInt       nleaves, l;
3988:     const PetscInt *leaves;
3989:     PetscInt       dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;

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

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

3996:       PetscSectionGetDof(globalFine,p,&dof);
3997:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3998:       if ((dof - cdof) > 0) {
3999:         numPointsWithDofs++;
4000:       }
4001:     }
4002:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
4003:     for (l = 0, offset = 0; l < nleaves; l++) {
4004:       PetscInt p = leaves ? leaves[l] : l;

4006:       PetscSectionGetDof(globalFine,p,&dof);
4007:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
4008:       if ((dof - cdof) > 0) {
4009:         pointsWithDofs[offset++] = l;
4010:       }
4011:     }
4012:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
4013:     PetscFree(pointsWithDofs);
4014:   }
4015:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
4016:   PetscMalloc1(pEndC-pStartC,&maxChildIds);
4017:   for (p = pStartC; p < pEndC; p++) {
4018:     maxChildIds[p - pStartC] = -2;
4019:   }
4020:   PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);
4021:   PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);

4023:   DMGetDefaultSection(coarse,&localCoarse);
4024:   DMGetDefaultGlobalSection(coarse,&globalCoarse);

4026:   DMPlexGetAnchors(coarse,&aSec,&aIS);
4027:   ISGetIndices(aIS,&anchors);
4028:   PetscSectionGetChart(aSec,&aStart,&aEnd);

4030:   DMGetDefaultConstraints(coarse,&cSec,&cMat);
4031:   PetscSectionGetChart(cSec,&cStart,&cEnd);

4033:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
4034:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);
4035:   PetscSectionSetChart(rootValuesSec,pStartC,pEndC);
4036:   PetscSectionGetNumFields(localCoarse,&numFields);
4037:   {
4038:     PetscInt maxFields = PetscMax(1,numFields) + 1;
4039:     PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);
4040:   }
4041:   if (grad) {
4042:     PetscInt i;

4044:     VecGetDM(cellGeom,&cellDM);
4045:     VecGetArrayRead(cellGeom,&cellGeomArray);
4046:     VecGetDM(grad,&gradDM);
4047:     VecGetArrayRead(grad,&gradArray);
4048:     for (i = 0; i < PetscMax(1,numFields); i++) {
4049:       PetscObject  obj;
4050:       PetscClassId id;

4052:       DMGetField(coarse, i, &obj);
4053:       PetscObjectGetClassId(obj,&id);
4054:       if (id == PETSCFV_CLASSID) {
4055:         fv      = (PetscFV) obj;
4056:         PetscFVGetNumComponents(fv,&numFVcomps);
4057:         fvField = i;
4058:         break;
4059:       }
4060:     }
4061:   }

4063:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
4064:     PetscInt dof;
4065:     PetscInt maxChildId     = maxChildIds[p - pStartC];
4066:     PetscInt numValues      = 0;

4068:     PetscSectionGetDof(globalCoarse,p,&dof);
4069:     if (dof < 0) {
4070:       dof = -(dof + 1);
4071:     }
4072:     offsets[0]    = 0;
4073:     newOffsets[0] = 0;
4074:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
4075:       PetscInt *closure = NULL, closureSize, cl;

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

4081:         PetscSectionGetDof(localCoarse,c,&clDof);
4082:         numValues += clDof;
4083:       }
4084:       DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
4085:     }
4086:     else if (maxChildId == -1) {
4087:       PetscSectionGetDof(localCoarse,p,&numValues);
4088:     }
4089:     /* we will pack the column indices with the field offsets */
4090:     if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
4091:       /* also send the centroid, and the gradient */
4092:       numValues += dim * (1 + numFVcomps);
4093:     }
4094:     PetscSectionSetDof(rootValuesSec,p,numValues);
4095:   }
4096:   PetscSectionSetUp(rootValuesSec);
4097:   {
4098:     PetscInt          numRootValues;
4099:     const PetscScalar *coarseArray;

4101:     PetscSectionGetStorageSize(rootValuesSec,&numRootValues);
4102:     PetscMalloc1(numRootValues,&rootValues);
4103:     VecGetArrayRead(vecCoarseLocal,&coarseArray);
4104:     for (p = pStartC; p < pEndC; p++) {
4105:       PetscInt    numValues;
4106:       PetscInt    pValOff;
4107:       PetscScalar *pVal;
4108:       PetscInt    maxChildId = maxChildIds[p - pStartC];

4110:       PetscSectionGetDof(rootValuesSec,p,&numValues);
4111:       if (!numValues) {
4112:         continue;
4113:       }
4114:       PetscSectionGetOffset(rootValuesSec,p,&pValOff);
4115:       pVal = &(rootValues[pValOff]);
4116:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
4117:         PetscInt closureSize = numValues;
4118:         DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);
4119:         if (grad && p >= cellStart && p < cellEnd) {
4120:           PetscFVCellGeom *cg;
4121:           PetscScalar     *gradVals = NULL;
4122:           PetscInt        i;

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

4126:           DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);
4127:           for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
4128:           pVal += dim;
4129:           DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);
4130:           for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
4131:         }
4132:       }
4133:       else if (maxChildId == -1) {
4134:         PetscInt lDof, lOff, i;

4136:         PetscSectionGetDof(localCoarse,p,&lDof);
4137:         PetscSectionGetOffset(localCoarse,p,&lOff);
4138:         for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
4139:       }
4140:     }
4141:     VecRestoreArrayRead(vecCoarseLocal,&coarseArray);
4142:     PetscFree(maxChildIds);
4143:   }
4144:   {
4145:     PetscSF  valuesSF;
4146:     PetscInt *remoteOffsetsValues, numLeafValues;

4148:     PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);
4149:     PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);
4150:     PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);
4151:     PetscSFDestroy(&coarseToFineEmbedded);
4152:     PetscFree(remoteOffsetsValues);
4153:     PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);
4154:     PetscMalloc1(numLeafValues,&leafValues);
4155:     PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues);
4156:     PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues);
4157:     PetscSFDestroy(&valuesSF);
4158:     PetscFree(rootValues);
4159:     PetscSectionDestroy(&rootValuesSec);
4160:   }
4161:   DMGetDefaultSection(fine,&localFine);
4162:   {
4163:     PetscInt    maxDof;
4164:     PetscInt    *rowIndices;
4165:     DM           refTree;
4166:     PetscInt     **refPointFieldN;
4167:     PetscScalar  ***refPointFieldMats;
4168:     PetscSection refConSec, refAnSec;
4169:     PetscInt     pRefStart,pRefEnd,leafStart,leafEnd;
4170:     PetscScalar  *pointWork;

4172:     PetscSectionGetMaxDof(globalFine,&maxDof);
4173:     DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
4174:     DMGetWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);
4175:     DMPlexGetReferenceTree(fine,&refTree);
4176:     DMGetDS(fine,&ds);
4177:     DMSetDS(refTree,ds);
4178:     DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
4179:     DMGetDefaultConstraints(refTree,&refConSec,NULL);
4180:     DMPlexGetAnchors(refTree,&refAnSec,NULL);
4181:     PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
4182:     PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);
4183:     DMPlexGetHeightStratum(fine,0,&cellStart,&cellEnd);
4184:     DMPlexGetHybridBounds(fine,&cellEndInterior,NULL,NULL,NULL);
4185:     cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior;
4186:     for (p = leafStart; p < leafEnd; p++) {
4187:       PetscInt          gDof, gcDof, gOff, lDof;
4188:       PetscInt          numValues, pValOff;
4189:       PetscInt          childId;
4190:       const PetscScalar *pVal;
4191:       const PetscScalar *fvGradData = NULL;

4193:       PetscSectionGetDof(globalFine,p,&gDof);
4194:       PetscSectionGetDof(localFine,p,&lDof);
4195:       PetscSectionGetConstraintDof(globalFine,p,&gcDof);
4196:       if ((gDof - gcDof) <= 0) {
4197:         continue;
4198:       }
4199:       PetscSectionGetOffset(globalFine,p,&gOff);
4200:       PetscSectionGetDof(leafValuesSec,p,&numValues);
4201:       if (!numValues) continue;
4202:       PetscSectionGetOffset(leafValuesSec,p,&pValOff);
4203:       pVal = &leafValues[pValOff];
4204:       offsets[0]        = 0;
4205:       offsetsCopy[0]    = 0;
4206:       newOffsets[0]     = 0;
4207:       newOffsetsCopy[0] = 0;
4208:       childId           = cids[p - pStartF];
4209:       if (numFields) {
4210:         PetscInt f;
4211:         for (f = 0; f < numFields; f++) {
4212:           PetscInt rowDof;

4214:           PetscSectionGetFieldDof(localFine,p,f,&rowDof);
4215:           offsets[f + 1]        = offsets[f] + rowDof;
4216:           offsetsCopy[f + 1]    = offsets[f + 1];
4217:           /* TODO: closure indices */
4218:           newOffsets[f + 1]     = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
4219:         }
4220:         DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);
4221:       }
4222:       else {
4223:         offsets[0]    = 0;
4224:         offsets[1]    = lDof;
4225:         newOffsets[0] = 0;
4226:         newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
4227:         DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);
4228:       }
4229:       if (childId == -1) { /* no child interpolation: one nnz per */
4230:         VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);
4231:       } else {
4232:         PetscInt f;

4234:         if (grad && p >= cellStart && p < cellEnd) {
4235:           numValues -= (dim * (1 + numFVcomps));
4236:           fvGradData = &pVal[numValues];
4237:         }
4238:         for (f = 0; f < PetscMax(1,numFields); f++) {
4239:           const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
4240:           PetscInt numRows = offsets[f+1] - offsets[f];
4241:           PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
4242:           const PetscScalar *cVal = &pVal[newOffsets[f]];
4243:           PetscScalar *rVal = &pointWork[offsets[f]];
4244:           PetscInt i, j;

4246: #if 0
4247:           PetscInfo6(coarse,"field %D, numFields %D, childId %D, numRows %D, numCols %D, refPointFieldN %D\n",f,numFields,childId,numRows,numCols,refPointFieldN[childId - pRefStart][f]);
4248: #endif
4249:           for (i = 0; i < numRows; i++) {
4250:             PetscScalar val = 0.;
4251:             for (j = 0; j < numCols; j++) {
4252:               val += childMat[i * numCols + j] * cVal[j];
4253:             }
4254:             rVal[i] = val;
4255:           }
4256:           if (f == fvField && p >= cellStart && p < cellEnd) {
4257:             PetscReal   centroid[3];
4258:             PetscScalar diff[3];
4259:             const PetscScalar *parentCentroid = &fvGradData[0];
4260:             const PetscScalar *gradient       = &fvGradData[dim];

4262:             DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);
4263:             for (i = 0; i < dim; i++) {
4264:               diff[i] = centroid[i] - parentCentroid[i];
4265:             }
4266:             for (i = 0; i < numFVcomps; i++) {
4267:               PetscScalar val = 0.;

4269:               for (j = 0; j < dim; j++) {
4270:                 val += gradient[dim * i + j] * diff[j];
4271:               }
4272:               rVal[i] += val;
4273:             }
4274:           }
4275:           VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);
4276:         }
4277:       }
4278:     }
4279:     DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
4280:     DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
4281:     DMRestoreWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);
4282:   }
4283:   PetscFree(leafValues);
4284:   PetscSectionDestroy(&leafValuesSec);
4285:   PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);
4286:   ISRestoreIndices(aIS,&anchors);
4287:   return(0);
4288: }

4290: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4291: {
4292:   PetscDS        ds;
4293:   DM             refTree;
4294:   PetscSection   multiRootSec, rootIndicesSec;
4295:   PetscSection   globalCoarse, globalFine;
4296:   PetscSection   localCoarse, localFine;
4297:   PetscSection   cSecRef;
4298:   PetscInt       *parentIndices, pRefStart, pRefEnd;
4299:   PetscScalar    *rootValues, *parentValues;
4300:   Mat            injRef;
4301:   PetscInt       numFields, maxDof;
4302:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
4303:   PetscInt       *offsets, *offsetsCopy, *rowOffsets;
4304:   PetscLayout    rowMap, colMap;
4305:   PetscInt       rowStart, rowEnd, colStart, colEnd;
4306:   PetscScalar    ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */


4311:   /* get the templates for the fine-to-coarse injection from the reference tree */
4312:   VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
4313:   VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
4314:   DMPlexGetReferenceTree(coarse,&refTree);
4315:   DMGetDS(coarse,&ds);
4316:   DMSetDS(refTree,ds);
4317:   DMGetDefaultConstraints(refTree,&cSecRef,NULL);
4318:   PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
4319:   DMPlexReferenceTreeGetInjector(refTree,&injRef);

4321:   DMPlexGetChart(fine,&pStartF,&pEndF);
4322:   DMGetDefaultSection(fine,&localFine);
4323:   DMGetDefaultGlobalSection(fine,&globalFine);
4324:   PetscSectionGetNumFields(localFine,&numFields);
4325:   DMPlexGetChart(coarse,&pStartC,&pEndC);
4326:   DMGetDefaultSection(coarse,&localCoarse);
4327:   DMGetDefaultGlobalSection(coarse,&globalCoarse);
4328:   PetscSectionGetMaxDof(localCoarse,&maxDof);
4329:   {
4330:     PetscInt maxFields = PetscMax(1,numFields) + 1;
4331:     PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
4332:   }

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

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

4338:   /* count indices */
4339:   VecGetLayout(vecFine,&colMap);
4340:   VecGetLayout(vecCoarse,&rowMap);
4341:   PetscLayoutSetUp(rowMap);
4342:   PetscLayoutSetUp(colMap);
4343:   PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
4344:   PetscLayoutGetRange(colMap,&colStart,&colEnd);
4345:   /* insert values */
4346:   DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
4347:   for (p = pStartC; p < pEndC; p++) {
4348:     PetscInt  numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4349:     PetscBool contribute = PETSC_FALSE;

4351:     PetscSectionGetDof(globalCoarse,p,&dof);
4352:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
4353:     if ((dof - cdof) <= 0) continue;
4354:     PetscSectionGetDof(localCoarse,p,&dof);
4355:     PetscSectionGetOffset(globalCoarse,p,&gOff);

4357:     rowOffsets[0] = 0;
4358:     offsetsCopy[0] = 0;
4359:     if (numFields) {
4360:       PetscInt f;

4362:       for (f = 0; f < numFields; f++) {
4363:         PetscInt fDof;
4364:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
4365:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
4366:       }
4367:       DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);
4368:     }
4369:     else {
4370:       DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);
4371:       rowOffsets[1] = offsetsCopy[0];
4372:     }

4374:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
4375:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
4376:     leafEnd = leafStart + numLeaves;
4377:     for (l = 0; l < dof; l++) parentValues[l] = 0.;
4378:     for (l = leafStart; l < leafEnd; l++) {
4379:       PetscInt numIndices, childId, offset;
4380:       const PetscScalar *childValues;

4382:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
4383:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
4384:       childId = (PetscInt) PetscRealPart(rootValues[offset++]);
4385:       childValues = &rootValues[offset];
4386:       numIndices--;

4388:       if (childId == -2) { /* skip */
4389:         continue;
4390:       } else if (childId == -1) { /* equivalent points: scatter */
4391:         PetscInt m;

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

4398:         contribute = PETSC_TRUE;
4399:         DMPlexGetTreeParent(refTree,childId,&parentId,NULL);

4401:         lim = PetscMax(1,numFields);
4402:         offsets[0] = 0;
4403:         if (numFields) {
4404:           PetscInt f;

4406:           for (f = 0; f < numFields; f++) {
4407:             PetscInt fDof;
4408:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

4410:             offsets[f + 1] = fDof + offsets[f];
4411:           }
4412:         }
4413:         else {
4414:           PetscInt cDof;

4416:           PetscSectionGetDof(cSecRef,childId,&cDof);
4417:           offsets[1] = cDof;
4418:         }
4419:         for (f = 0; f < lim; f++) {
4420:           PetscScalar       *childMat   = &childrenMats[childId - pRefStart][f][0];
4421:           PetscInt          n           = offsets[f+1]-offsets[f];
4422:           PetscInt          m           = rowOffsets[f+1]-rowOffsets[f];
4423:           PetscInt          i, j;
4424:           const PetscScalar *colValues  = &childValues[offsets[f]];

4426:           for (i = 0; i < m; i++) {
4427:             PetscScalar val = 0.;
4428:             for (j = 0; j < n; j++) {
4429:               val += childMat[n * i + j] * colValues[j];
4430:             }
4431:             parentValues[rowOffsets[f] + i] += val;
4432:           }
4433:         }
4434:       }
4435:     }
4436:     if (contribute) {VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);}
4437:   }
4438:   PetscSectionDestroy(&multiRootSec);
4439:   PetscSectionDestroy(&rootIndicesSec);
4440:   PetscFree2(parentIndices,parentValues);
4441:   DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
4442:   PetscFree(rootValues);
4443:   PetscFree3(offsets,offsetsCopy,rowOffsets);
4444:   return(0);
4445: }

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

4452:   collective

4454:   Input Parameters:
4455: + dmIn        - The DMPlex mesh for the input vector
4456: . vecIn       - The input vector
4457: . sfRefine    - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in
4458:                 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn
4459: . sfCoarsen   - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in
4460:                 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn
4461: . cidsRefine  - The childIds of the points in dmOut.  These childIds relate back to the reference tree: childid[j] = k implies
4462:                 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference
4463:                 tree was refined from its parent.  childid[j] = -1 indicates that the point j in dmOut is exactly
4464:                 equivalent to its root in dmIn, so no interpolation is necessary.  childid[j] = -2 indicates that this
4465:                 point j in dmOut is not a leaf of sfRefine.
4466: . cidsCoarsen - The childIds of the points in dmIn.  These childIds relate back to the reference tree: childid[j] = k implies
4467:                 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference
4468:                 tree coarsens to its parent.  childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen.
4469: . useBCs      - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer.
4470: - time        - Used if boundary values are time dependent.

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

4478:   Level: developer

4480: .seealso(): DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients()
4481: @*/
4482: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4483: {

4487:   VecSet(vecOut,0.0);
4488:   if (sfRefine) {
4489:     Vec vecInLocal;
4490:     DM  dmGrad = NULL;
4491:     Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;

4493:     DMGetLocalVector(dmIn,&vecInLocal);
4494:     VecSet(vecInLocal,0.0);
4495:     {
4496:       PetscInt  numFields, i;

4498:       DMGetNumFields(dmIn, &numFields);
4499:       for (i = 0; i < numFields; i++) {
4500:         PetscObject  obj;
4501:         PetscClassId classid;

4503:         DMGetField(dmIn, i, &obj);
4504:         PetscObjectGetClassId(obj, &classid);
4505:         if (classid == PETSCFV_CLASSID) {
4506:           DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);
4507:           break;
4508:         }
4509:       }
4510:     }
4511:     if (useBCs) {
4512:       DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);
4513:     }
4514:     DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);
4515:     DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);
4516:     if (dmGrad) {
4517:       DMGetGlobalVector(dmGrad,&grad);
4518:       DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);
4519:     }
4520:     DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);
4521:     DMRestoreLocalVector(dmIn,&vecInLocal);
4522:     if (dmGrad) {
4523:       DMRestoreGlobalVector(dmGrad,&grad);
4524:     }
4525:   }
4526:   if (sfCoarsen) {
4527:     DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);
4528:   }
4529:   VecAssemblyBegin(vecOut);
4530:   VecAssemblyEnd(vecOut);
4531:   return(0);
4532: }