Actual source code: plextree.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/isimpl.h>
  3: #include <petsc/private/petscfeimpl.h>
  4: #include <petscsf.h>
  5: #include <petscds.h>

  7: /* hierarchy routines */

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

 12:   Not collective

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

 18:   Level: intermediate

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

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

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

 37:   Not collective

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

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

 45:   Level: intermediate

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

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

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

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

 81:     DMPlexGetSupportSize(dm,childA,&size);
 82:     DMPlexGetSupport(dm,childA,&supp);

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

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

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

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

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

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

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

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

191:   Level: developer

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

201:   mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);
202:   return 0;
203: }

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

207: PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
208: {
209:   DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);
210:   return 0;
211: }

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

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: .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree()
450: @*/
451: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
452: {
453:   DM_Plex       *mesh;
454:   DM             K, Kref;
455:   PetscInt       p, pStart, pEnd;
456:   DMLabel        identity;

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

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

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

490:   PetscSectionDestroy(&mesh->childSection);
491:   PetscFree(mesh->children);
492:   pSec = mesh->parentSection;
493:   if (!pSec) return 0;
494:   PetscSectionGetStorageSize(pSec,&pSize);
495:   for (p = 0; p < pSize; p++) {
496:     PetscInt par = mesh->parents[p];

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

510:     PetscSectionAddDof(childSec,par,1);
511:   }
512:   PetscSectionSetUp(childSec);
513:   PetscSectionGetStorageSize(childSec,&cSize);
514:   PetscMalloc1(cSize,&children);
515:   PetscCalloc1(parMax-parMin,&offsets);
516:   PetscSectionGetChart(pSec,&pStart,&pEnd);
517:   for (p = pStart; p < pEnd; p++) {
518:     PetscInt dof, off, i;

520:     PetscSectionGetDof(pSec,p,&dof);
521:     PetscSectionGetOffset(pSec,p,&off);
522:     for (i = 0; i < dof; i++) {
523:       PetscInt par = mesh->parents[off + i], cOff;

525:       PetscSectionGetOffset(childSec,par,&cOff);
526:       children[cOff + offsets[par-parMin]++] = p;
527:     }
528:   }
529:   mesh->childSection = childSec;
530:   mesh->children = children;
531:   PetscFree(offsets);
532:   return 0;
533: }

535: static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
536: {
537:   PetscInt       pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
538:   const PetscInt *vals;
539:   PetscSection   secNew;
540:   PetscBool      anyNew, globalAnyNew;
541:   PetscBool      compress;

543:   PetscSectionGetChart(section,&pStart,&pEnd);
544:   ISGetLocalSize(is,&size);
545:   ISGetIndices(is,&vals);
546:   PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);
547:   PetscSectionSetChart(secNew,pStart,pEnd);
548:   for (i = 0; i < size; i++) {
549:     PetscInt dof;

551:     p = vals[i];
552:     if (p < pStart || p >= pEnd) continue;
553:     PetscSectionGetDof(section, p, &dof);
554:     if (dof) break;
555:   }
556:   if (i == size) {
557:     PetscSectionSetUp(secNew);
558:     anyNew   = PETSC_FALSE;
559:     compress = PETSC_FALSE;
560:     sizeNew  = 0;
561:   }
562:   else {
563:     anyNew = PETSC_TRUE;
564:     for (p = pStart; p < pEnd; p++) {
565:       PetscInt dof, off;

567:       PetscSectionGetDof(section, p, &dof);
568:       PetscSectionGetOffset(section, p, &off);
569:       for (i = 0; i < dof; i++) {
570:         PetscInt q = vals[off + i], qDof = 0;

572:         if (q >= pStart && q < pEnd) {
573:           PetscSectionGetDof(section, q, &qDof);
574:         }
575:         if (qDof) {
576:           PetscSectionAddDof(secNew, p, qDof);
577:         }
578:         else {
579:           PetscSectionAddDof(secNew, p, 1);
580:         }
581:       }
582:     }
583:     PetscSectionSetUp(secNew);
584:     PetscSectionGetStorageSize(secNew,&sizeNew);
585:     PetscMalloc1(sizeNew,&valsNew);
586:     compress = PETSC_FALSE;
587:     for (p = pStart; p < pEnd; p++) {
588:       PetscInt dof, off, count, offNew, dofNew;

590:       PetscSectionGetDof(section, p, &dof);
591:       PetscSectionGetOffset(section, p, &off);
592:       PetscSectionGetDof(secNew, p, &dofNew);
593:       PetscSectionGetOffset(secNew, p, &offNew);
594:       count = 0;
595:       for (i = 0; i < dof; i++) {
596:         PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;

598:         if (q >= pStart && q < pEnd) {
599:           PetscSectionGetDof(section, q, &qDof);
600:           PetscSectionGetOffset(section, q, &qOff);
601:         }
602:         if (qDof) {
603:           PetscInt oldCount = count;

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

608:             for (k = 0; k < oldCount; k++) {
609:               if (valsNew[offNew + k] == r) {
610:                 break;
611:               }
612:             }
613:             if (k == oldCount) {
614:               valsNew[offNew + count++] = r;
615:             }
616:           }
617:         }
618:         else {
619:           PetscInt k, oldCount = count;

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

647:     MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
648:     if (compress) {
649:       PetscSection secComp;
650:       PetscInt *valsComp = NULL;

652:       PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);
653:       PetscSectionSetChart(secComp,pStart,pEnd);
654:       for (p = pStart; p < pEnd; p++) {
655:         PetscInt dof;

657:         PetscSectionGetDof(secNew, p, &dof);
658:         PetscSectionSetDof(secComp, p, dof);
659:       }
660:       PetscSectionSetUp(secComp);
661:       PetscSectionGetStorageSize(secComp,&sizeNew);
662:       PetscMalloc1(sizeNew,&valsComp);
663:       for (p = pStart; p < pEnd; p++) {
664:         PetscInt dof, off, offNew, j;

666:         PetscSectionGetDof(secNew, p, &dof);
667:         PetscSectionGetOffset(secNew, p, &off);
668:         PetscSectionGetOffset(secComp, p, &offNew);
669:         for (j = 0; j < dof; j++) {
670:           valsComp[offNew + j] = valsNew[off + j];
671:         }
672:       }
673:       PetscSectionDestroy(&secNew);
674:       secNew  = secComp;
675:       PetscFree(valsNew);
676:       valsNew = valsComp;
677:     }
678:     ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);
679:   }
680:   return 0;
681: }

683: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
684: {
685:   PetscInt       p, pStart, pEnd, *anchors, size;
686:   PetscInt       aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
687:   PetscSection   aSec;
688:   DMLabel        canonLabel;
689:   IS             aIS;

692:   DMPlexGetChart(dm,&pStart,&pEnd);
693:   DMGetLabel(dm,"canonical",&canonLabel);
694:   for (p = pStart; p < pEnd; p++) {
695:     PetscInt parent;

697:     if (canonLabel) {
698:       PetscInt canon;

700:       DMLabelGetValue(canonLabel,p,&canon);
701:       if (p != canon) continue;
702:     }
703:     DMPlexGetTreeParent(dm,p,&parent,NULL);
704:     if (parent != p) {
705:       aMin = PetscMin(aMin,p);
706:       aMax = PetscMax(aMax,p+1);
707:     }
708:   }
709:   if (aMin > aMax) {
710:     aMin = -1;
711:     aMax = -1;
712:   }
713:   PetscSectionCreate(PETSC_COMM_SELF,&aSec);
714:   PetscSectionSetChart(aSec,aMin,aMax);
715:   for (p = aMin; p < aMax; p++) {
716:     PetscInt parent, ancestor = p;

718:     if (canonLabel) {
719:       PetscInt canon;

721:       DMLabelGetValue(canonLabel,p,&canon);
722:       if (p != canon) continue;
723:     }
724:     DMPlexGetTreeParent(dm,p,&parent,NULL);
725:     while (parent != ancestor) {
726:       ancestor = parent;
727:       DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
728:     }
729:     if (ancestor != p) {
730:       PetscInt closureSize, *closure = NULL;

732:       DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
733:       PetscSectionSetDof(aSec,p,closureSize);
734:       DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
735:     }
736:   }
737:   PetscSectionSetUp(aSec);
738:   PetscSectionGetStorageSize(aSec,&size);
739:   PetscMalloc1(size,&anchors);
740:   for (p = aMin; p < aMax; p++) {
741:     PetscInt parent, ancestor = p;

743:     if (canonLabel) {
744:       PetscInt canon;

746:       DMLabelGetValue(canonLabel,p,&canon);
747:       if (p != canon) continue;
748:     }
749:     DMPlexGetTreeParent(dm,p,&parent,NULL);
750:     while (parent != ancestor) {
751:       ancestor = parent;
752:       DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
753:     }
754:     if (ancestor != p) {
755:       PetscInt j, closureSize, *closure = NULL, aOff;

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

759:       DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
760:       for (j = 0; j < closureSize; j++) {
761:         anchors[aOff + j] = closure[2*j];
762:       }
763:       DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
764:     }
765:   }
766:   ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);
767:   {
768:     PetscSection aSecNew = aSec;
769:     IS           aISNew  = aIS;

771:     PetscObjectReference((PetscObject)aSec);
772:     PetscObjectReference((PetscObject)aIS);
773:     while (aSecNew) {
774:       PetscSectionDestroy(&aSec);
775:       ISDestroy(&aIS);
776:       aSec    = aSecNew;
777:       aIS     = aISNew;
778:       aSecNew = NULL;
779:       aISNew  = NULL;
780:       AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);
781:     }
782:   }
783:   DMPlexSetAnchors(dm,aSec,aIS);
784:   PetscSectionDestroy(&aSec);
785:   ISDestroy(&aIS);
786:   return 0;
787: }

789: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp)
790: {
791:   if (numTrueSupp[p] == -1) {
792:     PetscInt i, alldof;
793:     const PetscInt *supp;
794:     PetscInt count = 0;

796:     DMPlexGetSupportSize(dm,p,&alldof);
797:     DMPlexGetSupport(dm,p,&supp);
798:     for (i = 0; i < alldof; i++) {
799:       PetscInt q = supp[i], numCones, j;
800:       const PetscInt *cone;

802:       DMPlexGetConeSize(dm,q,&numCones);
803:       DMPlexGetCone(dm,q,&cone);
804:       for (j = 0; j < numCones; j++) {
805:         if (cone[j] == p) break;
806:       }
807:       if (j < numCones) count++;
808:     }
809:     numTrueSupp[p] = count;
810:   }
811:   *dof = numTrueSupp[p];
812:   return 0;
813: }

815: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
816: {
817:   DM_Plex        *mesh = (DM_Plex *)dm->data;
818:   PetscSection   newSupportSection;
819:   PetscInt       newSize, *newSupports, pStart, pEnd, p, d, depth;
820:   PetscInt       *numTrueSupp;
821:   PetscInt       *offsets;

824:   /* symmetrize the hierarchy */
825:   DMPlexGetDepth(dm,&depth);
826:   PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);
827:   DMPlexGetChart(dm,&pStart,&pEnd);
828:   PetscSectionSetChart(newSupportSection,pStart,pEnd);
829:   PetscCalloc1(pEnd,&offsets);
830:   PetscMalloc1(pEnd,&numTrueSupp);
831:   for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
832:   /* if a point is in the (true) support of q, it should be in the support of
833:    * parent(q) */
834:   for (d = 0; d <= depth; d++) {
835:     DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
836:     for (p = pStart; p < pEnd; ++p) {
837:       PetscInt dof, q, qdof, parent;

839:       DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);
840:       PetscSectionAddDof(newSupportSection, p, dof);
841:       q    = p;
842:       DMPlexGetTreeParent(dm,q,&parent,NULL);
843:       while (parent != q && parent >= pStart && parent < pEnd) {
844:         q = parent;

846:         DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);
847:         PetscSectionAddDof(newSupportSection,p,qdof);
848:         PetscSectionAddDof(newSupportSection,q,dof);
849:         DMPlexGetTreeParent(dm,q,&parent,NULL);
850:       }
851:     }
852:   }
853:   PetscSectionSetUp(newSupportSection);
854:   PetscSectionGetStorageSize(newSupportSection,&newSize);
855:   PetscMalloc1(newSize,&newSupports);
856:   for (d = 0; d <= depth; d++) {
857:     DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
858:     for (p = pStart; p < pEnd; p++) {
859:       PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;

861:       PetscSectionGetDof(mesh->supportSection, p, &dof);
862:       PetscSectionGetOffset(mesh->supportSection, p, &off);
863:       PetscSectionGetDof(newSupportSection, p, &newDof);
864:       PetscSectionGetOffset(newSupportSection, p, &newOff);
865:       for (i = 0; i < dof; i++) {
866:         PetscInt numCones, j;
867:         const PetscInt *cone;
868:         PetscInt q = mesh->supports[off + i];

870:         DMPlexGetConeSize(dm,q,&numCones);
871:         DMPlexGetCone(dm,q,&cone);
872:         for (j = 0; j < numCones; j++) {
873:           if (cone[j] == p) break;
874:         }
875:         if (j < numCones) newSupports[newOff+offsets[p]++] = q;
876:       }

878:       q    = p;
879:       DMPlexGetTreeParent(dm,q,&parent,NULL);
880:       while (parent != q && parent >= pStart && parent < pEnd) {
881:         q = parent;
882:         PetscSectionGetDof(mesh->supportSection, q, &qdof);
883:         PetscSectionGetOffset(mesh->supportSection, q, &qoff);
884:         PetscSectionGetOffset(newSupportSection, q, &newqOff);
885:         for (i = 0; i < qdof; i++) {
886:           PetscInt numCones, j;
887:           const PetscInt *cone;
888:           PetscInt r = mesh->supports[qoff + i];

890:           DMPlexGetConeSize(dm,r,&numCones);
891:           DMPlexGetCone(dm,r,&cone);
892:           for (j = 0; j < numCones; j++) {
893:             if (cone[j] == q) break;
894:           }
895:           if (j < numCones) newSupports[newOff+offsets[p]++] = r;
896:         }
897:         for (i = 0; i < dof; i++) {
898:           PetscInt numCones, j;
899:           const PetscInt *cone;
900:           PetscInt r = mesh->supports[off + i];

902:           DMPlexGetConeSize(dm,r,&numCones);
903:           DMPlexGetCone(dm,r,&cone);
904:           for (j = 0; j < numCones; j++) {
905:             if (cone[j] == p) break;
906:           }
907:           if (j < numCones) newSupports[newqOff+offsets[q]++] = r;
908:         }
909:         DMPlexGetTreeParent(dm,q,&parent,NULL);
910:       }
911:     }
912:   }
913:   PetscSectionDestroy(&mesh->supportSection);
914:   mesh->supportSection = newSupportSection;
915:   PetscFree(mesh->supports);
916:   mesh->supports = newSupports;
917:   PetscFree(offsets);
918:   PetscFree(numTrueSupp);

920:   return 0;
921: }

923: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat);
924: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat);

926: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
927: {
928:   DM_Plex       *mesh = (DM_Plex *)dm->data;
929:   DM             refTree;
930:   PetscInt       size;

934:   PetscObjectReference((PetscObject)parentSection);
935:   PetscSectionDestroy(&mesh->parentSection);
936:   mesh->parentSection = parentSection;
937:   PetscSectionGetStorageSize(parentSection,&size);
938:   if (parents != mesh->parents) {
939:     PetscFree(mesh->parents);
940:     PetscMalloc1(size,&mesh->parents);
941:     PetscArraycpy(mesh->parents, parents, size);
942:   }
943:   if (childIDs != mesh->childIDs) {
944:     PetscFree(mesh->childIDs);
945:     PetscMalloc1(size,&mesh->childIDs);
946:     PetscArraycpy(mesh->childIDs, childIDs, size);
947:   }
948:   DMPlexGetReferenceTree(dm,&refTree);
949:   if (refTree) {
950:     DMLabel canonLabel;

952:     DMGetLabel(refTree,"canonical",&canonLabel);
953:     if (canonLabel) {
954:       PetscInt i;

956:       for (i = 0; i < size; i++) {
957:         PetscInt canon;
958:         DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);
959:         if (canon >= 0) {
960:           mesh->childIDs[i] = canon;
961:         }
962:       }
963:     }
964:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
965:   } else {
966:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
967:   }
968:   DMPlexTreeSymmetrize(dm);
969:   if (computeCanonical) {
970:     PetscInt d, dim;

972:     /* add the canonical label */
973:     DMGetDimension(dm,&dim);
974:     DMCreateLabel(dm,"canonical");
975:     for (d = 0; d <= dim; d++) {
976:       PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
977:       const PetscInt *cChildren;

979:       DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
980:       for (p = dStart; p < dEnd; p++) {
981:         DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);
982:         if (cNumChildren) {
983:           canon = p;
984:           break;
985:         }
986:       }
987:       if (canon == -1) continue;
988:       for (p = dStart; p < dEnd; p++) {
989:         PetscInt numChildren, i;
990:         const PetscInt *children;

992:         DMPlexGetTreeChildren(dm,p,&numChildren,&children);
993:         if (numChildren) {
995:           DMSetLabelValue(dm,"canonical",p,canon);
996:           for (i = 0; i < numChildren; i++) {
997:             DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);
998:           }
999:         }
1000:       }
1001:     }
1002:   }
1003:   if (exchangeSupports) {
1004:     DMPlexTreeExchangeSupports(dm);
1005:   }
1006:   mesh->createanchors = DMPlexCreateAnchors_Tree;
1007:   /* reset anchors */
1008:   DMPlexSetAnchors(dm,NULL,NULL);
1009:   return 0;
1010: }

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

1017:   Collective on dm

1019:   Input Parameters:
1020: + dm - the DMPlex object
1021: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1022:                   offset indexes the parent and childID list; the reference count of parentSection is incremented
1023: . parents - a list of the point parents; copied, can be destroyed
1024: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1025:              the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed

1027:   Level: intermediate

1029: .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1030: @*/
1031: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
1032: {
1033:   DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);
1034:   return 0;
1035: }

1037: /*@
1038:   DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1039:   Collective on dm

1041:   Input Parameter:
1042: . dm - the DMPlex object

1044:   Output Parameters:
1045: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1046:                   offset indexes the parent and childID list
1047: . parents - a list of the point parents
1048: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1049:              the child corresponds to the point in the reference tree with index childID
1050: . childSection - the inverse of the parent section
1051: - children - a list of the point children

1053:   Level: intermediate

1055: .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1056: @*/
1057: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1058: {
1059:   DM_Plex        *mesh = (DM_Plex *)dm->data;

1062:   if (parentSection) *parentSection = mesh->parentSection;
1063:   if (parents)       *parents       = mesh->parents;
1064:   if (childIDs)      *childIDs      = mesh->childIDs;
1065:   if (childSection)  *childSection  = mesh->childSection;
1066:   if (children)      *children      = mesh->children;
1067:   return 0;
1068: }

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

1073:   Input Parameters:
1074: + dm - the DMPlex object
1075: - point - the query point

1077:   Output Parameters:
1078: + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent
1079: - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point
1080:             does not have a parent

1082:   Level: intermediate

1084: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren()
1085: @*/
1086: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1087: {
1088:   DM_Plex       *mesh = (DM_Plex *)dm->data;
1089:   PetscSection   pSec;

1092:   pSec = mesh->parentSection;
1093:   if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1094:     PetscInt dof;

1096:     PetscSectionGetDof (pSec, point, &dof);
1097:     if (dof) {
1098:       PetscInt off;

1100:       PetscSectionGetOffset (pSec, point, &off);
1101:       if (parent)  *parent = mesh->parents[off];
1102:       if (childID) *childID = mesh->childIDs[off];
1103:       return 0;
1104:     }
1105:   }
1106:   if (parent) {
1107:     *parent = point;
1108:   }
1109:   if (childID) {
1110:     *childID = 0;
1111:   }
1112:   return 0;
1113: }

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

1118:   Input Parameters:
1119: + dm - the DMPlex object
1120: - point - the query point

1122:   Output Parameters:
1123: + numChildren - if not NULL, set to the number of children
1124: - children - if not NULL, set to a list children, or set to NULL if the point has no children

1126:   Level: intermediate

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

1132: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent()
1133: @*/
1134: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1135: {
1136:   DM_Plex       *mesh = (DM_Plex *)dm->data;
1137:   PetscSection   childSec;
1138:   PetscInt       dof = 0;

1141:   childSec = mesh->childSection;
1142:   if (childSec && point >= childSec->pStart && point < childSec->pEnd) {
1143:     PetscSectionGetDof (childSec, point, &dof);
1144:   }
1145:   if (numChildren) *numChildren = dof;
1146:   if (children) {
1147:     if (dof) {
1148:       PetscInt off;

1150:       PetscSectionGetOffset (childSec, point, &off);
1151:       *children = &mesh->children[off];
1152:     }
1153:     else {
1154:       *children = NULL;
1155:     }
1156:   }
1157:   return 0;
1158: }

1160: 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)
1161: {
1162:   PetscInt       f, b, p, c, offset, qPoints;

1164:   PetscSpaceEvaluate(space,nPoints,points,work,NULL,NULL);
1165:   for (f = 0, offset = 0; f < nFunctionals; f++) {
1166:     qPoints = pointsPerFn[f];
1167:     for (b = 0; b < nBasis; b++) {
1168:       PetscScalar val = 0.;

1170:       for (p = 0; p < qPoints; p++) {
1171:         for (c = 0; c < nComps; c++) {
1172:           val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c];
1173:         }
1174:       }
1175:       MatSetValue(basisAtPoints,b,f,val,INSERT_VALUES);
1176:     }
1177:     offset += qPoints;
1178:   }
1179:   MatAssemblyBegin(basisAtPoints,MAT_FINAL_ASSEMBLY);
1180:   MatAssemblyEnd(basisAtPoints,MAT_FINAL_ASSEMBLY);
1181:   return 0;
1182: }

1184: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
1185: {
1186:   PetscDS        ds;
1187:   PetscInt       spdim;
1188:   PetscInt       numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
1189:   const PetscInt *anchors;
1190:   PetscSection   aSec;
1191:   PetscReal      *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
1192:   IS             aIS;

1194:   DMPlexGetChart(dm,&pStart,&pEnd);
1195:   DMGetDS(dm,&ds);
1196:   PetscDSGetNumFields(ds,&numFields);
1197:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1198:   DMPlexGetAnchors(dm,&aSec,&aIS);
1199:   ISGetIndices(aIS,&anchors);
1200:   PetscSectionGetChart(cSec,&conStart,&conEnd);
1201:   DMGetDimension(dm,&spdim);
1202:   PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);

1204:   for (f = 0; f < numFields; f++) {
1205:     PetscObject       disc;
1206:     PetscClassId      id;
1207:     PetscSpace        bspace;
1208:     PetscDualSpace    dspace;
1209:     PetscInt          i, j, k, nPoints, Nc, offset;
1210:     PetscInt          fSize, maxDof;
1211:     PetscReal         *weights, *pointsRef, *pointsReal, *work;
1212:     PetscScalar       *scwork;
1213:     const PetscScalar *X;
1214:     PetscInt          *sizes, *workIndRow, *workIndCol;
1215:     Mat               Amat, Bmat, Xmat;
1216:     const PetscInt    *numDof  = NULL;
1217:     const PetscInt    ***perms = NULL;
1218:     const PetscScalar ***flips = NULL;

1220:     PetscDSGetDiscretization(ds,f,&disc);
1221:     PetscObjectGetClassId(disc,&id);
1222:     if (id == PETSCFE_CLASSID) {
1223:       PetscFE fe = (PetscFE) disc;

1225:       PetscFEGetBasisSpace(fe,&bspace);
1226:       PetscFEGetDualSpace(fe,&dspace);
1227:       PetscDualSpaceGetDimension(dspace,&fSize);
1228:       PetscFEGetNumComponents(fe,&Nc);
1229:     }
1230:     else if (id == PETSCFV_CLASSID) {
1231:       PetscFV fv = (PetscFV) disc;

1233:       PetscFVGetNumComponents(fv,&Nc);
1234:       PetscSpaceCreate(PetscObjectComm((PetscObject)fv),&bspace);
1235:       PetscSpaceSetType(bspace,PETSCSPACEPOLYNOMIAL);
1236:       PetscSpaceSetDegree(bspace,0,PETSC_DETERMINE);
1237:       PetscSpaceSetNumComponents(bspace,Nc);
1238:       PetscSpaceSetNumVariables(bspace,spdim);
1239:       PetscSpaceSetUp(bspace);
1240:       PetscFVGetDualSpace(fv,&dspace);
1241:       PetscDualSpaceGetDimension(dspace,&fSize);
1242:     }
1243:     else SETERRQ(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1244:     PetscDualSpaceGetNumDof(dspace,&numDof);
1245:     for (i = 0, maxDof = 0; i <= spdim; i++) {maxDof = PetscMax(maxDof,numDof[i]);}
1246:     PetscDualSpaceGetSymmetries(dspace,&perms,&flips);

1248:     MatCreate(PETSC_COMM_SELF,&Amat);
1249:     MatSetSizes(Amat,fSize,fSize,fSize,fSize);
1250:     MatSetType(Amat,MATSEQDENSE);
1251:     MatSetUp(Amat);
1252:     MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);
1253:     MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);
1254:     nPoints = 0;
1255:     for (i = 0; i < fSize; i++) {
1256:       PetscInt        qPoints, thisNc;
1257:       PetscQuadrature quad;

1259:       PetscDualSpaceGetFunctional(dspace,i,&quad);
1260:       PetscQuadratureGetData(quad,NULL,&thisNc,&qPoints,NULL,NULL);
1262:       nPoints += qPoints;
1263:     }
1264:     PetscMalloc7(fSize,&sizes,nPoints*Nc,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal,nPoints*fSize*Nc,&work,maxDof,&workIndRow,maxDof,&workIndCol);
1265:     PetscMalloc1(maxDof * maxDof,&scwork);
1266:     offset = 0;
1267:     for (i = 0; i < fSize; i++) {
1268:       PetscInt        qPoints;
1269:       const PetscReal    *p, *w;
1270:       PetscQuadrature quad;

1272:       PetscDualSpaceGetFunctional(dspace,i,&quad);
1273:       PetscQuadratureGetData(quad,NULL,NULL,&qPoints,&p,&w);
1274:       PetscArraycpy(weights+Nc*offset,w,Nc*qPoints);
1275:       PetscArraycpy(pointsRef+spdim*offset,p,spdim*qPoints);
1276:       sizes[i] = qPoints;
1277:       offset  += qPoints;
1278:     }
1279:     EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsRef,weights,work,Amat);
1280:     MatLUFactor(Amat,NULL,NULL,NULL);
1281:     for (c = cStart; c < cEnd; c++) {
1282:       PetscInt        parent;
1283:       PetscInt        closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1284:       PetscInt        *childOffsets, *parentOffsets;

1286:       DMPlexGetTreeParent(dm,c,&parent,NULL);
1287:       if (parent == c) continue;
1288:       DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1289:       for (i = 0; i < closureSize; i++) {
1290:         PetscInt p = closure[2*i];
1291:         PetscInt conDof;

1293:         if (p < conStart || p >= conEnd) continue;
1294:         if (numFields) {
1295:           PetscSectionGetFieldDof(cSec,p,f,&conDof);
1296:         }
1297:         else {
1298:           PetscSectionGetDof(cSec,p,&conDof);
1299:         }
1300:         if (conDof) break;
1301:       }
1302:       if (i == closureSize) {
1303:         DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1304:         continue;
1305:       }

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

1312:         CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i*spdim],vtmp);
1313:         CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);
1314:       }
1315:       EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsReal,weights,work,Bmat);
1316:       MatMatSolve(Amat,Bmat,Xmat);
1317:       MatDenseGetArrayRead(Xmat,&X);
1318:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1319:       PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);
1320:       childOffsets[0] = 0;
1321:       for (i = 0; i < closureSize; i++) {
1322:         PetscInt p = closure[2*i];
1323:         PetscInt dof;

1325:         if (numFields) {
1326:           PetscSectionGetFieldDof(section,p,f,&dof);
1327:         }
1328:         else {
1329:           PetscSectionGetDof(section,p,&dof);
1330:         }
1331:         childOffsets[i+1]=childOffsets[i]+dof;
1332:       }
1333:       parentOffsets[0] = 0;
1334:       for (i = 0; i < closureSizeP; i++) {
1335:         PetscInt p = closureP[2*i];
1336:         PetscInt dof;

1338:         if (numFields) {
1339:           PetscSectionGetFieldDof(section,p,f,&dof);
1340:         }
1341:         else {
1342:           PetscSectionGetDof(section,p,&dof);
1343:         }
1344:         parentOffsets[i+1]=parentOffsets[i]+dof;
1345:       }
1346:       for (i = 0; i < closureSize; i++) {
1347:         PetscInt conDof, conOff, aDof, aOff, nWork;
1348:         PetscInt p = closure[2*i];
1349:         PetscInt o = closure[2*i+1];
1350:         const PetscInt    *perm;
1351:         const PetscScalar *flip;

1353:         if (p < conStart || p >= conEnd) continue;
1354:         if (numFields) {
1355:           PetscSectionGetFieldDof(cSec,p,f,&conDof);
1356:           PetscSectionGetFieldOffset(cSec,p,f,&conOff);
1357:         }
1358:         else {
1359:           PetscSectionGetDof(cSec,p,&conDof);
1360:           PetscSectionGetOffset(cSec,p,&conOff);
1361:         }
1362:         if (!conDof) continue;
1363:         perm  = (perms && perms[i]) ? perms[i][o] : NULL;
1364:         flip  = (flips && flips[i]) ? flips[i][o] : NULL;
1365:         PetscSectionGetDof(aSec,p,&aDof);
1366:         PetscSectionGetOffset(aSec,p,&aOff);
1367:         nWork = childOffsets[i+1]-childOffsets[i];
1368:         for (k = 0; k < aDof; k++) {
1369:           PetscInt a = anchors[aOff + k];
1370:           PetscInt aSecDof, aSecOff;

1372:           if (numFields) {
1373:             PetscSectionGetFieldDof(section,a,f,&aSecDof);
1374:             PetscSectionGetFieldOffset(section,a,f,&aSecOff);
1375:           }
1376:           else {
1377:             PetscSectionGetDof(section,a,&aSecDof);
1378:             PetscSectionGetOffset(section,a,&aSecOff);
1379:           }
1380:           if (!aSecDof) continue;

1382:           for (j = 0; j < closureSizeP; j++) {
1383:             PetscInt q = closureP[2*j];
1384:             PetscInt oq = closureP[2*j+1];

1386:             if (q == a) {
1387:               PetscInt           r, s, nWorkP;
1388:               const PetscInt    *permP;
1389:               const PetscScalar *flipP;

1391:               permP  = (perms && perms[j]) ? perms[j][oq] : NULL;
1392:               flipP  = (flips && flips[j]) ? flips[j][oq] : NULL;
1393:               nWorkP = parentOffsets[j+1]-parentOffsets[j];
1394:               /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the
1395:                * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArrayRead() is
1396:                * column-major, so transpose-transpose = do nothing */
1397:               for (r = 0; r < nWork; r++) {
1398:                 for (s = 0; s < nWorkP; s++) {
1399:                   scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])];
1400:                 }
1401:               }
1402:               for (r = 0; r < nWork; r++)  {workIndRow[perm  ? perm[r]  : r] = conOff  + r;}
1403:               for (s = 0; s < nWorkP; s++) {workIndCol[permP ? permP[s] : s] = aSecOff + s;}
1404:               if (flip) {
1405:                 for (r = 0; r < nWork; r++) {
1406:                   for (s = 0; s < nWorkP; s++) {
1407:                     scwork[r * nWorkP + s] *= flip[r];
1408:                   }
1409:                 }
1410:               }
1411:               if (flipP) {
1412:                 for (r = 0; r < nWork; r++) {
1413:                   for (s = 0; s < nWorkP; s++) {
1414:                     scwork[r * nWorkP + s] *= flipP[s];
1415:                   }
1416:                 }
1417:               }
1418:               MatSetValues(cMat,nWork,workIndRow,nWorkP,workIndCol,scwork,INSERT_VALUES);
1419:               break;
1420:             }
1421:           }
1422:         }
1423:       }
1424:       MatDenseRestoreArrayRead(Xmat,&X);
1425:       PetscFree2(childOffsets,parentOffsets);
1426:       DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1427:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1428:     }
1429:     MatDestroy(&Amat);
1430:     MatDestroy(&Bmat);
1431:     MatDestroy(&Xmat);
1432:     PetscFree(scwork);
1433:     PetscFree7(sizes,weights,pointsRef,pointsReal,work,workIndRow,workIndCol);
1434:     if (id == PETSCFV_CLASSID) {
1435:       PetscSpaceDestroy(&bspace);
1436:     }
1437:   }
1438:   MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1439:   MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1440:   PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);
1441:   ISRestoreIndices(aIS,&anchors);

1443:   return 0;
1444: }

1446: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1447: {
1448:   Mat               refCmat;
1449:   PetscDS           ds;
1450:   PetscInt          numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1451:   PetscScalar       ***refPointFieldMats;
1452:   PetscSection      refConSec, refAnSec, refSection;
1453:   IS                refAnIS;
1454:   const PetscInt    *refAnchors;
1455:   const PetscInt    **perms;
1456:   const PetscScalar **flips;

1458:   DMGetDS(refTree,&ds);
1459:   PetscDSGetNumFields(ds,&numFields);
1460:   maxFields = PetscMax(1,numFields);
1461:   DMGetDefaultConstraints(refTree,&refConSec,&refCmat,NULL);
1462:   DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1463:   ISGetIndices(refAnIS,&refAnchors);
1464:   DMGetLocalSection(refTree,&refSection);
1465:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1466:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
1467:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);
1468:   PetscSectionGetMaxDof(refConSec,&maxDof);
1469:   PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1470:   PetscMalloc1(maxDof,&rows);
1471:   PetscMalloc1(maxDof*maxAnDof,&cols);
1472:   for (p = pRefStart; p < pRefEnd; p++) {
1473:     PetscInt parent, closureSize, *closure = NULL, pDof;

1475:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
1476:     PetscSectionGetDof(refConSec,p,&pDof);
1477:     if (!pDof || parent == p) continue;

1479:     PetscMalloc1(maxFields,&refPointFieldMats[p-pRefStart]);
1480:     PetscCalloc1(maxFields,&refPointFieldN[p-pRefStart]);
1481:     DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1482:     for (f = 0; f < maxFields; f++) {
1483:       PetscInt cDof, cOff, numCols, r, i;

1485:       if (f < numFields) {
1486:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1487:         PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
1488:         PetscSectionGetFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);
1489:       } else {
1490:         PetscSectionGetDof(refConSec,p,&cDof);
1491:         PetscSectionGetOffset(refConSec,p,&cOff);
1492:         PetscSectionGetPointSyms(refSection,closureSize,closure,&perms,&flips);
1493:       }

1495:       for (r = 0; r < cDof; r++) {
1496:         rows[r] = cOff + r;
1497:       }
1498:       numCols = 0;
1499:       for (i = 0; i < closureSize; i++) {
1500:         PetscInt          q = closure[2*i];
1501:         PetscInt          aDof, aOff, j;
1502:         const PetscInt    *perm = perms ? perms[i] : NULL;

1504:         if (numFields) {
1505:           PetscSectionGetFieldDof(refSection,q,f,&aDof);
1506:           PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1507:         }
1508:         else {
1509:           PetscSectionGetDof(refSection,q,&aDof);
1510:           PetscSectionGetOffset(refSection,q,&aOff);
1511:         }

1513:         for (j = 0; j < aDof; j++) {
1514:           cols[numCols++] = aOff + (perm ? perm[j] : j);
1515:         }
1516:       }
1517:       refPointFieldN[p-pRefStart][f] = numCols;
1518:       PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
1519:       MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);
1520:       if (flips) {
1521:         PetscInt colOff = 0;

1523:         for (i = 0; i < closureSize; i++) {
1524:           PetscInt          q = closure[2*i];
1525:           PetscInt          aDof, aOff, j;
1526:           const PetscScalar *flip = flips ? flips[i] : NULL;

1528:           if (numFields) {
1529:             PetscSectionGetFieldDof(refSection,q,f,&aDof);
1530:             PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1531:           }
1532:           else {
1533:             PetscSectionGetDof(refSection,q,&aDof);
1534:             PetscSectionGetOffset(refSection,q,&aOff);
1535:           }
1536:           if (flip) {
1537:             PetscInt k;
1538:             for (k = 0; k < cDof; k++) {
1539:               for (j = 0; j < aDof; j++) {
1540:                 refPointFieldMats[p-pRefStart][f][k * numCols + colOff + j] *= flip[j];
1541:               }
1542:             }
1543:           }
1544:           colOff += aDof;
1545:         }
1546:       }
1547:       if (numFields) {
1548:         PetscSectionRestoreFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);
1549:       } else {
1550:         PetscSectionRestorePointSyms(refSection,closureSize,closure,&perms,&flips);
1551:       }
1552:     }
1553:     DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1554:   }
1555:   *childrenMats = refPointFieldMats;
1556:   *childrenN = refPointFieldN;
1557:   ISRestoreIndices(refAnIS,&refAnchors);
1558:   PetscFree(rows);
1559:   PetscFree(cols);
1560:   return 0;
1561: }

1563: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1564: {
1565:   PetscDS        ds;
1566:   PetscInt       **refPointFieldN;
1567:   PetscScalar    ***refPointFieldMats;
1568:   PetscInt       numFields, maxFields, pRefStart, pRefEnd, p, f;
1569:   PetscSection   refConSec;

1571:   refPointFieldN = *childrenN;
1572:   *childrenN = NULL;
1573:   refPointFieldMats = *childrenMats;
1574:   *childrenMats = NULL;
1575:   DMGetDS(refTree,&ds);
1576:   PetscDSGetNumFields(ds,&numFields);
1577:   maxFields = PetscMax(1,numFields);
1578:   DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL);
1579:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1580:   for (p = pRefStart; p < pRefEnd; p++) {
1581:     PetscInt parent, pDof;

1583:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
1584:     PetscSectionGetDof(refConSec,p,&pDof);
1585:     if (!pDof || parent == p) continue;

1587:     for (f = 0; f < maxFields; f++) {
1588:       PetscInt cDof;

1590:       if (numFields) {
1591:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1592:       }
1593:       else {
1594:         PetscSectionGetDof(refConSec,p,&cDof);
1595:       }

1597:       PetscFree(refPointFieldMats[p - pRefStart][f]);
1598:     }
1599:     PetscFree(refPointFieldMats[p - pRefStart]);
1600:     PetscFree(refPointFieldN[p - pRefStart]);
1601:   }
1602:   PetscFree(refPointFieldMats);
1603:   PetscFree(refPointFieldN);
1604:   return 0;
1605: }

1607: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1608: {
1609:   DM             refTree;
1610:   PetscDS        ds;
1611:   Mat            refCmat;
1612:   PetscInt       numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1613:   PetscScalar ***refPointFieldMats, *pointWork;
1614:   PetscSection   refConSec, refAnSec, anSec;
1615:   IS             refAnIS, anIS;
1616:   const PetscInt *anchors;

1619:   DMGetDS(dm,&ds);
1620:   PetscDSGetNumFields(ds,&numFields);
1621:   maxFields = PetscMax(1,numFields);
1622:   DMPlexGetReferenceTree(dm,&refTree);
1623:   DMCopyDisc(dm,refTree);
1624:   DMGetDefaultConstraints(refTree,&refConSec,&refCmat,NULL);
1625:   DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1626:   DMPlexGetAnchors(dm,&anSec,&anIS);
1627:   ISGetIndices(anIS,&anchors);
1628:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1629:   PetscSectionGetChart(conSec,&conStart,&conEnd);
1630:   PetscSectionGetMaxDof(refConSec,&maxDof);
1631:   PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1632:   PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);

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

1637:   /* step 2: compute the preorder */
1638:   DMPlexGetChart(dm,&pStart,&pEnd);
1639:   PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);
1640:   for (p = pStart; p < pEnd; p++) {
1641:     perm[p - pStart] = p;
1642:     iperm[p - pStart] = p-pStart;
1643:   }
1644:   for (p = 0; p < pEnd - pStart;) {
1645:     PetscInt point = perm[p];
1646:     PetscInt parent;

1648:     DMPlexGetTreeParent(dm,point,&parent,NULL);
1649:     if (parent == point) {
1650:       p++;
1651:     }
1652:     else {
1653:       PetscInt size, closureSize, *closure = NULL, i;

1655:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1656:       for (i = 0; i < closureSize; i++) {
1657:         PetscInt q = closure[2*i];
1658:         if (iperm[q-pStart] > iperm[point-pStart]) {
1659:           /* swap */
1660:           perm[p]               = q;
1661:           perm[iperm[q-pStart]] = point;
1662:           iperm[point-pStart]   = iperm[q-pStart];
1663:           iperm[q-pStart]       = p;
1664:           break;
1665:         }
1666:       }
1667:       size = closureSize;
1668:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1669:       if (i == size) {
1670:         p++;
1671:       }
1672:     }
1673:   }

1675:   /* step 3: fill the constraint matrix */
1676:   /* we are going to use a preorder progressive fill strategy.  Mat doesn't
1677:    * allow progressive fill without assembly, so we are going to set up the
1678:    * values outside of the Mat first.
1679:    */
1680:   {
1681:     PetscInt nRows, row, nnz;
1682:     PetscBool done;
1683:     const PetscInt *ia, *ja;
1684:     PetscScalar *vals;

1686:     MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1688:     nnz  = ia[nRows];
1689:     /* malloc and then zero rows right before we fill them: this way valgrind
1690:      * can tell if we are doing progressive fill in the wrong order */
1691:     PetscMalloc1(nnz,&vals);
1692:     for (p = 0; p < pEnd - pStart; p++) {
1693:       PetscInt        parent, childid, closureSize, *closure = NULL;
1694:       PetscInt        point = perm[p], pointDof;

1696:       DMPlexGetTreeParent(dm,point,&parent,&childid);
1697:       if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1698:       PetscSectionGetDof(conSec,point,&pointDof);
1699:       if (!pointDof) continue;
1700:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1701:       for (f = 0; f < maxFields; f++) {
1702:         PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
1703:         PetscScalar *pointMat;
1704:         const PetscInt    **perms;
1705:         const PetscScalar **flips;

1707:         if (numFields) {
1708:           PetscSectionGetFieldDof(conSec,point,f,&cDof);
1709:           PetscSectionGetFieldOffset(conSec,point,f,&cOff);
1710:         }
1711:         else {
1712:           PetscSectionGetDof(conSec,point,&cDof);
1713:           PetscSectionGetOffset(conSec,point,&cOff);
1714:         }
1715:         if (!cDof) continue;
1716:         if (numFields) PetscSectionGetFieldPointSyms(section,f,closureSize,closure,&perms,&flips);
1717:         else           PetscSectionGetPointSyms(section,closureSize,closure,&perms,&flips);

1719:         /* make sure that every row for this point is the same size */
1720:         if (PetscDefined(USE_DEBUG)) {
1721:           for (r = 0; r < cDof; r++) {
1722:             if (cDof > 1 && r) {
1724:             }
1725:           }
1726:         }
1727:         /* zero rows */
1728:         for (i = ia[cOff] ; i< ia[cOff+cDof];i++) {
1729:           vals[i] = 0.;
1730:         }
1731:         matOffset = ia[cOff];
1732:         numFillCols = ia[cOff+1] - matOffset;
1733:         pointMat = refPointFieldMats[childid-pRefStart][f];
1734:         numCols = refPointFieldN[childid-pRefStart][f];
1735:         offset = 0;
1736:         for (i = 0; i < closureSize; i++) {
1737:           PetscInt q = closure[2*i];
1738:           PetscInt aDof, aOff, j, k, qConDof, qConOff;
1739:           const PetscInt    *perm = perms ? perms[i] : NULL;
1740:           const PetscScalar *flip = flips ? flips[i] : NULL;

1742:           qConDof = qConOff = 0;
1743:           if (numFields) {
1744:             PetscSectionGetFieldDof(section,q,f,&aDof);
1745:             PetscSectionGetFieldOffset(section,q,f,&aOff);
1746:             if (q >= conStart && q < conEnd) {
1747:               PetscSectionGetFieldDof(conSec,q,f,&qConDof);
1748:               PetscSectionGetFieldOffset(conSec,q,f,&qConOff);
1749:             }
1750:           }
1751:           else {
1752:             PetscSectionGetDof(section,q,&aDof);
1753:             PetscSectionGetOffset(section,q,&aOff);
1754:             if (q >= conStart && q < conEnd) {
1755:               PetscSectionGetDof(conSec,q,&qConDof);
1756:               PetscSectionGetOffset(conSec,q,&qConOff);
1757:             }
1758:           }
1759:           if (!aDof) continue;
1760:           if (qConDof) {
1761:             /* this point has anchors: its rows of the matrix should already
1762:              * be filled, thanks to preordering */
1763:             /* first multiply into pointWork, then set in matrix */
1764:             PetscInt aMatOffset = ia[qConOff];
1765:             PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1766:             for (r = 0; r < cDof; r++) {
1767:               for (j = 0; j < aNumFillCols; j++) {
1768:                 PetscScalar inVal = 0;
1769:                 for (k = 0; k < aDof; k++) {
1770:                   PetscInt col = perm ? perm[k] : k;

1772:                   inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.);
1773:                 }
1774:                 pointWork[r * aNumFillCols + j] = inVal;
1775:               }
1776:             }
1777:             /* assume that the columns are sorted, spend less time searching */
1778:             for (j = 0, k = 0; j < aNumFillCols; j++) {
1779:               PetscInt col = ja[aMatOffset + j];
1780:               for (;k < numFillCols; k++) {
1781:                 if (ja[matOffset + k] == col) {
1782:                   break;
1783:                 }
1784:               }
1786:               for (r = 0; r < cDof; r++) {
1787:                 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1788:               }
1789:             }
1790:           }
1791:           else {
1792:             /* find where to put this portion of pointMat into the matrix */
1793:             for (k = 0; k < numFillCols; k++) {
1794:               if (ja[matOffset + k] == aOff) {
1795:                 break;
1796:               }
1797:             }
1799:             for (r = 0; r < cDof; r++) {
1800:               for (j = 0; j < aDof; j++) {
1801:                 PetscInt col = perm ? perm[j] : j;

1803:                 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
1804:               }
1805:             }
1806:           }
1807:           offset += aDof;
1808:         }
1809:         if (numFields) {
1810:           PetscSectionRestoreFieldPointSyms(section,f,closureSize,closure,&perms,&flips);
1811:         } else {
1812:           PetscSectionRestorePointSyms(section,closureSize,closure,&perms,&flips);
1813:         }
1814:       }
1815:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1816:     }
1817:     for (row = 0; row < nRows; row++) {
1818:       MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);
1819:     }
1820:     MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1822:     MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1823:     MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1824:     PetscFree(vals);
1825:   }

1827:   /* clean up */
1828:   ISRestoreIndices(anIS,&anchors);
1829:   PetscFree2(perm,iperm);
1830:   PetscFree(pointWork);
1831:   DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
1832:   return 0;
1833: }

1835: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1836:  * a non-conforming mesh.  Local refinement comes later */
1837: PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm)
1838: {
1839:   DM K;
1840:   PetscMPIInt rank;
1841:   PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1842:   PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1843:   PetscInt *Kembedding;
1844:   PetscInt *cellClosure=NULL, nc;
1845:   PetscScalar *newVertexCoords;
1846:   PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1847:   PetscSection parentSection;

1849:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1850:   DMGetDimension(dm,&dim);
1851:   DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);
1852:   DMSetDimension(*ncdm,dim);

1854:   DMPlexGetChart(dm, &pStart, &pEnd);
1855:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);
1856:   DMPlexGetReferenceTree(dm,&K);
1857:   if (rank == 0) {
1858:     /* compute the new charts */
1859:     PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);
1860:     offset = 0;
1861:     for (d = 0; d <= dim; d++) {
1862:       PetscInt pOldCount, kStart, kEnd, k;

1864:       pNewStart[d] = offset;
1865:       DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);
1866:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1867:       pOldCount = pOldEnd[d] - pOldStart[d];
1868:       /* adding the new points */
1869:       pNewCount[d] = pOldCount + kEnd - kStart;
1870:       if (!d) {
1871:         /* removing the cell */
1872:         pNewCount[d]--;
1873:       }
1874:       for (k = kStart; k < kEnd; k++) {
1875:         PetscInt parent;
1876:         DMPlexGetTreeParent(K,k,&parent,NULL);
1877:         if (parent == k) {
1878:           /* avoid double counting points that won't actually be new */
1879:           pNewCount[d]--;
1880:         }
1881:       }
1882:       pNewEnd[d] = pNewStart[d] + pNewCount[d];
1883:       offset = pNewEnd[d];

1885:     }
1887:     /* get the current closure of the cell that we are removing */
1888:     DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);

1890:     PetscMalloc1(pNewEnd[dim],&newConeSizes);
1891:     {
1892:       DMPolytopeType pct, qct;
1893:       PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;

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

1898:       for (k = kStart; k < kEnd; k++) {
1899:         perm[k - kStart] = k;
1900:         iperm [k - kStart] = k - kStart;
1901:         preOrient[k - kStart] = 0;
1902:       }

1904:       DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1905:       for (j = 1; j < closureSizeK; j++) {
1906:         PetscInt parentOrientA = closureK[2*j+1];
1907:         PetscInt parentOrientB = cellClosure[2*j+1];
1908:         PetscInt p, q;

1910:         p = closureK[2*j];
1911:         q = cellClosure[2*j];
1912:         DMPlexGetCellType(K, p, &pct);
1913:         DMPlexGetCellType(dm, q, &qct);
1914:         for (d = 0; d <= dim; d++) {
1915:           if (q >= pOldStart[d] && q < pOldEnd[d]) {
1916:             Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1917:           }
1918:         }
1919:         parentOrientA = DMPolytopeConvertNewOrientation_Internal(pct, parentOrientA);
1920:         parentOrientB = DMPolytopeConvertNewOrientation_Internal(qct, parentOrientB);
1921:         if (parentOrientA != parentOrientB) {
1922:           PetscInt numChildren, i;
1923:           const PetscInt *children;

1925:           DMPlexGetTreeChildren(K,p,&numChildren,&children);
1926:           for (i = 0; i < numChildren; i++) {
1927:             PetscInt kPerm, oPerm;

1929:             k    = children[i];
1930:             DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);
1931:             /* perm = what refTree position I'm in */
1932:             perm[kPerm-kStart]      = k;
1933:             /* iperm = who is at this position */
1934:             iperm[k-kStart]         = kPerm-kStart;
1935:             preOrient[kPerm-kStart] = oPerm;
1936:           }
1937:         }
1938:       }
1939:       DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1940:     }
1941:     PetscSectionSetChart(parentSection,0,pNewEnd[dim]);
1942:     offset = 0;
1943:     numNewCones = 0;
1944:     for (d = 0; d <= dim; d++) {
1945:       PetscInt kStart, kEnd, k;
1946:       PetscInt p;
1947:       PetscInt size;

1949:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1950:         /* skip cell 0 */
1951:         if (p == cell) continue;
1952:         /* old cones to new cones */
1953:         DMPlexGetConeSize(dm,p,&size);
1954:         newConeSizes[offset++] = size;
1955:         numNewCones += size;
1956:       }

1958:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1959:       for (k = kStart; k < kEnd; k++) {
1960:         PetscInt kParent;

1962:         DMPlexGetTreeParent(K,k,&kParent,NULL);
1963:         if (kParent != k) {
1964:           Kembedding[k] = offset;
1965:           DMPlexGetConeSize(K,k,&size);
1966:           newConeSizes[offset++] = size;
1967:           numNewCones += size;
1968:           if (kParent != 0) {
1969:             PetscSectionSetDof(parentSection,Kembedding[k],1);
1970:           }
1971:         }
1972:       }
1973:     }

1975:     PetscSectionSetUp(parentSection);
1976:     PetscSectionGetStorageSize(parentSection,&numPointsWithParents);
1977:     PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);
1978:     PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);

1980:     /* fill new cones */
1981:     offset = 0;
1982:     for (d = 0; d <= dim; d++) {
1983:       PetscInt kStart, kEnd, k, l;
1984:       PetscInt p;
1985:       PetscInt size;
1986:       const PetscInt *cone, *orientation;

1988:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1989:         /* skip cell 0 */
1990:         if (p == cell) continue;
1991:         /* old cones to new cones */
1992:         DMPlexGetConeSize(dm,p,&size);
1993:         DMPlexGetCone(dm,p,&cone);
1994:         DMPlexGetConeOrientation(dm,p,&orientation);
1995:         for (l = 0; l < size; l++) {
1996:           newCones[offset]          = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
1997:           newOrientations[offset++] = orientation[l];
1998:         }
1999:       }

2001:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
2002:       for (k = kStart; k < kEnd; k++) {
2003:         PetscInt kPerm = perm[k], kParent;
2004:         PetscInt preO  = preOrient[k];

2006:         DMPlexGetTreeParent(K,k,&kParent,NULL);
2007:         if (kParent != k) {
2008:           /* embed new cones */
2009:           DMPlexGetConeSize(K,k,&size);
2010:           DMPlexGetCone(K,kPerm,&cone);
2011:           DMPlexGetConeOrientation(K,kPerm,&orientation);
2012:           for (l = 0; l < size; l++) {
2013:             PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size);
2014:             PetscInt newO, lSize, oTrue;
2015:             DMPolytopeType ct = DM_NUM_POLYTOPES;

2017:             q                         = iperm[cone[m]];
2018:             newCones[offset]          = Kembedding[q];
2019:             DMPlexGetConeSize(K,q,&lSize);
2020:             if (lSize == 2) ct = DM_POLYTOPE_SEGMENT;
2021:             else if (lSize == 4) ct = DM_POLYTOPE_QUADRILATERAL;
2022:             oTrue                     = DMPolytopeConvertNewOrientation_Internal(ct, orientation[m]);
2023:             oTrue                     = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
2024:             newO                      = DihedralCompose(lSize,oTrue,preOrient[q]);
2025:             newOrientations[offset++] = DMPolytopeConvertOldOrientation_Internal(ct, newO);
2026:           }
2027:           if (kParent != 0) {
2028:             PetscInt newPoint = Kembedding[kParent];
2029:             PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);
2030:             parents[pOffset]  = newPoint;
2031:             childIDs[pOffset] = k;
2032:           }
2033:         }
2034:       }
2035:     }

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

2039:     /* fill coordinates */
2040:     offset = 0;
2041:     {
2042:       PetscInt kStart, kEnd, l;
2043:       PetscSection vSection;
2044:       PetscInt v;
2045:       Vec coords;
2046:       PetscScalar *coordvals;
2047:       PetscInt dof, off;
2048:       PetscReal v0[3], J[9], detJ;

2050:       if (PetscDefined(USE_DEBUG)) {
2051:         PetscInt k;
2052:         DMPlexGetHeightStratum(K,0,&kStart,&kEnd);
2053:         for (k = kStart; k < kEnd; k++) {
2054:           DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);
2056:         }
2057:       }
2058:       DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);
2059:       DMGetCoordinateSection(dm,&vSection);
2060:       DMGetCoordinatesLocal(dm,&coords);
2061:       VecGetArray(coords,&coordvals);
2062:       for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {

2064:         PetscSectionGetDof(vSection,v,&dof);
2065:         PetscSectionGetOffset(vSection,v,&off);
2066:         for (l = 0; l < dof; l++) {
2067:           newVertexCoords[offset++] = coordvals[off + l];
2068:         }
2069:       }
2070:       VecRestoreArray(coords,&coordvals);

2072:       DMGetCoordinateSection(K,&vSection);
2073:       DMGetCoordinatesLocal(K,&coords);
2074:       VecGetArray(coords,&coordvals);
2075:       DMPlexGetDepthStratum(K,0,&kStart,&kEnd);
2076:       for (v = kStart; v < kEnd; v++) {
2077:         PetscReal coord[3], newCoord[3];
2078:         PetscInt  vPerm = perm[v];
2079:         PetscInt  kParent;
2080:         const PetscReal xi0[3] = {-1.,-1.,-1.};

2082:         DMPlexGetTreeParent(K,v,&kParent,NULL);
2083:         if (kParent != v) {
2084:           /* this is a new vertex */
2085:           PetscSectionGetOffset(vSection,vPerm,&off);
2086:           for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]);
2087:           CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
2088:           for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l];
2089:           offset += dim;
2090:         }
2091:       }
2092:       VecRestoreArray(coords,&coordvals);
2093:     }

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

2099:       tmp = pNewCount[d];
2100:       pNewCount[d] = pNewCount[dim - d];
2101:       pNewCount[dim - d] = tmp;
2102:     }

2104:     DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);
2105:     DMPlexSetReferenceTree(*ncdm,K);
2106:     DMPlexSetTree(*ncdm,parentSection,parents,childIDs);

2108:     /* clean up */
2109:     DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
2110:     PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);
2111:     PetscFree(newConeSizes);
2112:     PetscFree2(newCones,newOrientations);
2113:     PetscFree(newVertexCoords);
2114:     PetscFree2(parents,childIDs);
2115:     PetscFree4(Kembedding,perm,iperm,preOrient);
2116:   }
2117:   else {
2118:     PetscInt    p, counts[4];
2119:     PetscInt    *coneSizes, *cones, *orientations;
2120:     Vec         coordVec;
2121:     PetscScalar *coords;

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

2126:       DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
2127:       counts[d] = dEnd - dStart;
2128:     }
2129:     PetscMalloc1(pEnd-pStart,&coneSizes);
2130:     for (p = pStart; p < pEnd; p++) {
2131:       DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);
2132:     }
2133:     DMPlexGetCones(dm, &cones);
2134:     DMPlexGetConeOrientations(dm, &orientations);
2135:     DMGetCoordinatesLocal(dm,&coordVec);
2136:     VecGetArray(coordVec,&coords);

2138:     PetscSectionSetChart(parentSection,pStart,pEnd);
2139:     PetscSectionSetUp(parentSection);
2140:     DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);
2141:     DMPlexSetReferenceTree(*ncdm,K);
2142:     DMPlexSetTree(*ncdm,parentSection,NULL,NULL);
2143:     VecRestoreArray(coordVec,&coords);
2144:   }
2145:   PetscSectionDestroy(&parentSection);

2147:   return 0;
2148: }

2150: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2151: {
2152:   PetscSF           coarseToFineEmbedded;
2153:   PetscSection      globalCoarse, globalFine;
2154:   PetscSection      localCoarse, localFine;
2155:   PetscSection      aSec, cSec;
2156:   PetscSection      rootIndicesSec, rootMatricesSec;
2157:   PetscSection      leafIndicesSec, leafMatricesSec;
2158:   PetscInt          *rootIndices, *leafIndices;
2159:   PetscScalar       *rootMatrices, *leafMatrices;
2160:   IS                aIS;
2161:   const PetscInt    *anchors;
2162:   Mat               cMat;
2163:   PetscInt          numFields, maxFields;
2164:   PetscInt          pStartC, pEndC, pStartF, pEndF, p;
2165:   PetscInt          aStart, aEnd, cStart, cEnd;
2166:   PetscInt          *maxChildIds;
2167:   PetscInt          *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2168:   const PetscInt    ***perms;
2169:   const PetscScalar ***flips;

2171:   DMPlexGetChart(coarse,&pStartC,&pEndC);
2172:   DMPlexGetChart(fine,&pStartF,&pEndF);
2173:   DMGetGlobalSection(fine,&globalFine);
2174:   { /* winnow fine points that don't have global dofs out of the sf */
2175:     PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2176:     const PetscInt *leaves;

2178:     PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);
2179:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2180:       p = leaves ? leaves[l] : l;
2181:       PetscSectionGetDof(globalFine,p,&dof);
2182:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
2183:       if ((dof - cdof) > 0) {
2184:         numPointsWithDofs++;
2185:       }
2186:     }
2187:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
2188:     for (l = 0, offset = 0; l < nleaves; l++) {
2189:       p = leaves ? leaves[l] : l;
2190:       PetscSectionGetDof(globalFine,p,&dof);
2191:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
2192:       if ((dof - cdof) > 0) {
2193:         pointsWithDofs[offset++] = l;
2194:       }
2195:     }
2196:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
2197:     PetscFree(pointsWithDofs);
2198:   }
2199:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2200:   PetscMalloc1(pEndC-pStartC,&maxChildIds);
2201:   for (p = pStartC; p < pEndC; p++) {
2202:     maxChildIds[p - pStartC] = -2;
2203:   }
2204:   PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);
2205:   PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);

2207:   DMGetLocalSection(coarse,&localCoarse);
2208:   DMGetGlobalSection(coarse,&globalCoarse);

2210:   DMPlexGetAnchors(coarse,&aSec,&aIS);
2211:   ISGetIndices(aIS,&anchors);
2212:   PetscSectionGetChart(aSec,&aStart,&aEnd);

2214:   DMGetDefaultConstraints(coarse,&cSec,&cMat,NULL);
2215:   PetscSectionGetChart(cSec,&cStart,&cEnd);

2217:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2218:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
2219:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);
2220:   PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);
2221:   PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);
2222:   PetscSectionGetNumFields(localCoarse,&numFields);
2223:   maxFields = PetscMax(1,numFields);
2224:   PetscMalloc7(maxFields+1,&offsets,maxFields+1,&offsetsCopy,maxFields+1,&newOffsets,maxFields+1,&newOffsetsCopy,maxFields+1,&rowOffsets,maxFields+1,&numD,maxFields+1,&numO);
2225:   PetscMalloc2(maxFields+1,(PetscInt****)&perms,maxFields+1,(PetscScalar****)&flips);
2226:   PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **));
2227:   PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **));

2229:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2230:     PetscInt dof, matSize   = 0;
2231:     PetscInt aDof           = 0;
2232:     PetscInt cDof           = 0;
2233:     PetscInt maxChildId     = maxChildIds[p - pStartC];
2234:     PetscInt numRowIndices  = 0;
2235:     PetscInt numColIndices  = 0;
2236:     PetscInt f;

2238:     PetscSectionGetDof(globalCoarse,p,&dof);
2239:     if (dof < 0) {
2240:       dof = -(dof + 1);
2241:     }
2242:     if (p >= aStart && p < aEnd) {
2243:       PetscSectionGetDof(aSec,p,&aDof);
2244:     }
2245:     if (p >= cStart && p < cEnd) {
2246:       PetscSectionGetDof(cSec,p,&cDof);
2247:     }
2248:     for (f = 0; f <= numFields; f++) offsets[f] = 0;
2249:     for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2250:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2251:       PetscInt *closure = NULL, closureSize, cl;

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

2257:         PetscSectionGetDof(localCoarse,c,&clDof);
2258:         numRowIndices += clDof;
2259:         for (f = 0; f < numFields; f++) {
2260:           PetscSectionGetFieldDof(localCoarse,c,f,&clDof);
2261:           offsets[f + 1] += clDof;
2262:         }
2263:       }
2264:       for (f = 0; f < numFields; f++) {
2265:         offsets[f + 1]   += offsets[f];
2266:         newOffsets[f + 1] = offsets[f + 1];
2267:       }
2268:       /* get the number of indices needed and their field offsets */
2269:       DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);
2270:       DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
2271:       if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2272:         numColIndices = numRowIndices;
2273:         matSize = 0;
2274:       }
2275:       else if (numFields) { /* we send one submat for each field: sum their sizes */
2276:         matSize = 0;
2277:         for (f = 0; f < numFields; f++) {
2278:           PetscInt numRow, numCol;

2280:           numRow = offsets[f + 1] - offsets[f];
2281:           numCol = newOffsets[f + 1] - newOffsets[f];
2282:           matSize += numRow * numCol;
2283:         }
2284:       }
2285:       else {
2286:         matSize = numRowIndices * numColIndices;
2287:       }
2288:     } else if (maxChildId == -1) {
2289:       if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2290:         PetscInt aOff, a;

2292:         PetscSectionGetOffset(aSec,p,&aOff);
2293:         for (f = 0; f < numFields; f++) {
2294:           PetscInt fDof;

2296:           PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2297:           offsets[f+1] = fDof;
2298:         }
2299:         for (a = 0; a < aDof; a++) {
2300:           PetscInt anchor = anchors[a + aOff], aLocalDof;

2302:           PetscSectionGetDof(localCoarse,anchor,&aLocalDof);
2303:           numColIndices += aLocalDof;
2304:           for (f = 0; f < numFields; f++) {
2305:             PetscInt fDof;

2307:             PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2308:             newOffsets[f+1] += fDof;
2309:           }
2310:         }
2311:         if (numFields) {
2312:           matSize = 0;
2313:           for (f = 0; f < numFields; f++) {
2314:             matSize += offsets[f+1] * newOffsets[f+1];
2315:           }
2316:         }
2317:         else {
2318:           matSize = numColIndices * dof;
2319:         }
2320:       }
2321:       else { /* no children, and no constraints on dofs: just get the global indices */
2322:         numColIndices = dof;
2323:         matSize       = 0;
2324:       }
2325:     }
2326:     /* we will pack the column indices with the field offsets */
2327:     PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);
2328:     PetscSectionSetDof(rootMatricesSec,p,matSize);
2329:   }
2330:   PetscSectionSetUp(rootIndicesSec);
2331:   PetscSectionSetUp(rootMatricesSec);
2332:   {
2333:     PetscInt numRootIndices, numRootMatrices;

2335:     PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
2336:     PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);
2337:     PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);
2338:     for (p = pStartC; p < pEndC; p++) {
2339:       PetscInt    numRowIndices, numColIndices, matSize, dof;
2340:       PetscInt    pIndOff, pMatOff, f;
2341:       PetscInt    *pInd;
2342:       PetscInt    maxChildId = maxChildIds[p - pStartC];
2343:       PetscScalar *pMat = NULL;

2345:       PetscSectionGetDof(rootIndicesSec,p,&numColIndices);
2346:       if (!numColIndices) {
2347:         continue;
2348:       }
2349:       for (f = 0; f <= numFields; f++) {
2350:         offsets[f]        = 0;
2351:         newOffsets[f]     = 0;
2352:         offsetsCopy[f]    = 0;
2353:         newOffsetsCopy[f] = 0;
2354:       }
2355:       numColIndices -= 2 * numFields;
2356:       PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);
2357:       pInd = &(rootIndices[pIndOff]);
2358:       PetscSectionGetDof(rootMatricesSec,p,&matSize);
2359:       if (matSize) {
2360:         PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);
2361:         pMat = &rootMatrices[pMatOff];
2362:       }
2363:       PetscSectionGetDof(globalCoarse,p,&dof);
2364:       if (dof < 0) {
2365:         dof = -(dof + 1);
2366:       }
2367:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2368:         PetscInt i, j;
2369:         PetscInt numRowIndices = matSize / numColIndices;

2371:         if (!numRowIndices) { /* don't need to calculate the mat, just the indices */
2372:           PetscInt numIndices, *indices;
2373:           DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,PETSC_TRUE,&numIndices,&indices,offsets,NULL);
2375:           for (i = 0; i < numColIndices; i++) {
2376:             pInd[i] = indices[i];
2377:           }
2378:           for (i = 0; i < numFields; i++) {
2379:             pInd[numColIndices + i]             = offsets[i+1];
2380:             pInd[numColIndices + numFields + i] = offsets[i+1];
2381:           }
2382:           DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,PETSC_TRUE,&numIndices,&indices,offsets,NULL);
2383:         }
2384:         else {
2385:           PetscInt closureSize, *closure = NULL, cl;
2386:           PetscScalar *pMatIn, *pMatModified;
2387:           PetscInt numPoints,*points;

2389:           DMGetWorkArray(coarse,numRowIndices * numRowIndices,MPIU_SCALAR,&pMatIn);
2390:           for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2391:             for (j = 0; j < numRowIndices; j++) {
2392:               pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2393:             }
2394:           }
2395:           DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2396:           for (f = 0; f < maxFields; f++) {
2397:             if (numFields) PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);
2398:             else           PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);
2399:           }
2400:           if (numFields) {
2401:             for (cl = 0; cl < closureSize; cl++) {
2402:               PetscInt c = closure[2 * cl];

2404:               for (f = 0; f < numFields; f++) {
2405:                 PetscInt fDof;

2407:                 PetscSectionGetFieldDof(localCoarse,c,f,&fDof);
2408:                 offsets[f + 1] += fDof;
2409:               }
2410:             }
2411:             for (f = 0; f < numFields; f++) {
2412:               offsets[f + 1]   += offsets[f];
2413:               newOffsets[f + 1] = offsets[f + 1];
2414:             }
2415:           }
2416:           /* TODO : flips here ? */
2417:           /* apply hanging node constraints on the right, get the new points and the new offsets */
2418:           DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);
2419:           for (f = 0; f < maxFields; f++) {
2420:             if (numFields) PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);
2421:             else           PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);
2422:           }
2423:           for (f = 0; f < maxFields; f++) {
2424:             if (numFields) PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);
2425:             else           PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);
2426:           }
2427:           if (!numFields) {
2428:             for (i = 0; i < numRowIndices * numColIndices; i++) {
2429:               pMat[i] = pMatModified[i];
2430:             }
2431:           }
2432:           else {
2433:             PetscInt i, j, count;
2434:             for (f = 0, count = 0; f < numFields; f++) {
2435:               for (i = offsets[f]; i < offsets[f+1]; i++) {
2436:                 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) {
2437:                   pMat[count] = pMatModified[i * numColIndices + j];
2438:                 }
2439:               }
2440:             }
2441:           }
2442:           DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatModified);
2443:           DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2444:           DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatIn);
2445:           if (numFields) {
2446:             for (f = 0; f < numFields; f++) {
2447:               pInd[numColIndices + f]             = offsets[f+1];
2448:               pInd[numColIndices + numFields + f] = newOffsets[f+1];
2449:             }
2450:             for (cl = 0; cl < numPoints; cl++) {
2451:               PetscInt globalOff, c = points[2*cl];
2452:               PetscSectionGetOffset(globalCoarse, c, &globalOff);
2453:               DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd);
2454:             }
2455:           } else {
2456:             for (cl = 0; cl < numPoints; cl++) {
2457:               PetscInt c = points[2*cl], globalOff;
2458:               const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;

2460:               PetscSectionGetOffset(globalCoarse, c, &globalOff);
2461:               DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd);
2462:             }
2463:           }
2464:           for (f = 0; f < maxFields; f++) {
2465:             if (numFields) PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);
2466:             else           PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);
2467:           }
2468:           DMRestoreWorkArray(coarse,numPoints,MPIU_SCALAR,&points);
2469:         }
2470:       }
2471:       else if (matSize) {
2472:         PetscInt cOff;
2473:         PetscInt *rowIndices, *colIndices, a, aDof, aOff;

2475:         numRowIndices = matSize / numColIndices;
2477:         DMGetWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);
2478:         DMGetWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);
2479:         PetscSectionGetOffset(cSec,p,&cOff);
2480:         PetscSectionGetDof(aSec,p,&aDof);
2481:         PetscSectionGetOffset(aSec,p,&aOff);
2482:         if (numFields) {
2483:           for (f = 0; f < numFields; f++) {
2484:             PetscInt fDof;

2486:             PetscSectionGetFieldDof(cSec,p,f,&fDof);
2487:             offsets[f + 1] = fDof;
2488:             for (a = 0; a < aDof; a++) {
2489:               PetscInt anchor = anchors[a + aOff];
2490:               PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2491:               newOffsets[f + 1] += fDof;
2492:             }
2493:           }
2494:           for (f = 0; f < numFields; f++) {
2495:             offsets[f + 1]       += offsets[f];
2496:             offsetsCopy[f + 1]    = offsets[f + 1];
2497:             newOffsets[f + 1]    += newOffsets[f];
2498:             newOffsetsCopy[f + 1] = newOffsets[f + 1];
2499:           }
2500:           DMPlexGetIndicesPointFields_Internal(cSec,PETSC_TRUE,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1, NULL,rowIndices);
2501:           for (a = 0; a < aDof; a++) {
2502:             PetscInt anchor = anchors[a + aOff], lOff;
2503:             PetscSectionGetOffset(localCoarse,anchor,&lOff);
2504:             DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_TRUE,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1, NULL,colIndices);
2505:           }
2506:         }
2507:         else {
2508:           DMPlexGetIndicesPoint_Internal(cSec,PETSC_TRUE,p,cOff,offsetsCopy,PETSC_TRUE,NULL, NULL,rowIndices);
2509:           for (a = 0; a < aDof; a++) {
2510:             PetscInt anchor = anchors[a + aOff], lOff;
2511:             PetscSectionGetOffset(localCoarse,anchor,&lOff);
2512:             DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_TRUE,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL, NULL,colIndices);
2513:           }
2514:         }
2515:         if (numFields) {
2516:           PetscInt count, a;

2518:           for (f = 0, count = 0; f < numFields; f++) {
2519:             PetscInt iSize = offsets[f + 1] - offsets[f];
2520:             PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2521:             MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);
2522:             count += iSize * jSize;
2523:             pInd[numColIndices + f]             = offsets[f+1];
2524:             pInd[numColIndices + numFields + f] = newOffsets[f+1];
2525:           }
2526:           for (a = 0; a < aDof; a++) {
2527:             PetscInt anchor = anchors[a + aOff];
2528:             PetscInt gOff;
2529:             PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2530:             DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1, NULL,pInd);
2531:           }
2532:         }
2533:         else {
2534:           PetscInt a;
2535:           MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);
2536:           for (a = 0; a < aDof; a++) {
2537:             PetscInt anchor = anchors[a + aOff];
2538:             PetscInt gOff;
2539:             PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2540:             DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL, NULL,pInd);
2541:           }
2542:         }
2543:         DMRestoreWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);
2544:         DMRestoreWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);
2545:       }
2546:       else {
2547:         PetscInt gOff;

2549:         PetscSectionGetOffset(globalCoarse,p,&gOff);
2550:         if (numFields) {
2551:           for (f = 0; f < numFields; f++) {
2552:             PetscInt fDof;
2553:             PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2554:             offsets[f + 1] = fDof + offsets[f];
2555:           }
2556:           for (f = 0; f < numFields; f++) {
2557:             pInd[numColIndices + f]             = offsets[f+1];
2558:             pInd[numColIndices + numFields + f] = offsets[f+1];
2559:           }
2560:           DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1, NULL,pInd);
2561:         } else {
2562:           DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL, NULL,pInd);
2563:         }
2564:       }
2565:     }
2566:     PetscFree(maxChildIds);
2567:   }
2568:   {
2569:     PetscSF  indicesSF, matricesSF;
2570:     PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;

2572:     PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
2573:     PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);
2574:     PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);
2575:     PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);
2576:     PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);
2577:     PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);
2578:     PetscSFDestroy(&coarseToFineEmbedded);
2579:     PetscFree(remoteOffsetsIndices);
2580:     PetscFree(remoteOffsetsMatrices);
2581:     PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);
2582:     PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);
2583:     PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);
2584:     PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices,MPI_REPLACE);
2585:     PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices,MPI_REPLACE);
2586:     PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices,MPI_REPLACE);
2587:     PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices,MPI_REPLACE);
2588:     PetscSFDestroy(&matricesSF);
2589:     PetscSFDestroy(&indicesSF);
2590:     PetscFree2(rootIndices,rootMatrices);
2591:     PetscSectionDestroy(&rootIndicesSec);
2592:     PetscSectionDestroy(&rootMatricesSec);
2593:   }
2594:   /* count to preallocate */
2595:   DMGetLocalSection(fine,&localFine);
2596:   {
2597:     PetscInt    nGlobal;
2598:     PetscInt    *dnnz, *onnz;
2599:     PetscLayout rowMap, colMap;
2600:     PetscInt    rowStart, rowEnd, colStart, colEnd;
2601:     PetscInt    maxDof;
2602:     PetscInt    *rowIndices;
2603:     DM           refTree;
2604:     PetscInt     **refPointFieldN;
2605:     PetscScalar  ***refPointFieldMats;
2606:     PetscSection refConSec, refAnSec;
2607:     PetscInt     pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd;
2608:     PetscScalar  *pointWork;

2610:     PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);
2611:     PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);
2612:     MatGetLayouts(mat,&rowMap,&colMap);
2613:     PetscLayoutSetUp(rowMap);
2614:     PetscLayoutSetUp(colMap);
2615:     PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
2616:     PetscLayoutGetRange(colMap,&colStart,&colEnd);
2617:     PetscSectionGetMaxDof(localFine,&maxDof);
2618:     PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd);
2619:     DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
2620:     for (p = leafStart; p < leafEnd; p++) {
2621:       PetscInt    gDof, gcDof, gOff;
2622:       PetscInt    numColIndices, pIndOff, *pInd;
2623:       PetscInt    matSize;
2624:       PetscInt    i;

2626:       PetscSectionGetDof(globalFine,p,&gDof);
2627:       PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2628:       if ((gDof - gcDof) <= 0) {
2629:         continue;
2630:       }
2631:       PetscSectionGetOffset(globalFine,p,&gOff);
2634:       PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2635:       PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2636:       numColIndices -= 2 * numFields;
2638:       pInd = &leafIndices[pIndOff];
2639:       offsets[0]        = 0;
2640:       offsetsCopy[0]    = 0;
2641:       newOffsets[0]     = 0;
2642:       newOffsetsCopy[0] = 0;
2643:       if (numFields) {
2644:         PetscInt f;
2645:         for (f = 0; f < numFields; f++) {
2646:           PetscInt rowDof;

2648:           PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2649:           offsets[f + 1]        = offsets[f] + rowDof;
2650:           offsetsCopy[f + 1]    = offsets[f + 1];
2651:           newOffsets[f + 1]     = pInd[numColIndices + numFields + f];
2652:           numD[f] = 0;
2653:           numO[f] = 0;
2654:         }
2655:         DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,rowIndices);
2656:         for (f = 0; f < numFields; f++) {
2657:           PetscInt colOffset    = newOffsets[f];
2658:           PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];

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

2663:             if (gInd >= colStart && gInd < colEnd) {
2664:               numD[f]++;
2665:             }
2666:             else if (gInd >= 0) { /* negative means non-entry */
2667:               numO[f]++;
2668:             }
2669:           }
2670:         }
2671:       }
2672:       else {
2673:         DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,rowIndices);
2674:         numD[0] = 0;
2675:         numO[0] = 0;
2676:         for (i = 0; i < numColIndices; i++) {
2677:           PetscInt gInd = pInd[i];

2679:           if (gInd >= colStart && gInd < colEnd) {
2680:             numD[0]++;
2681:           }
2682:           else if (gInd >= 0) { /* negative means non-entry */
2683:             numO[0]++;
2684:           }
2685:         }
2686:       }
2687:       PetscSectionGetDof(leafMatricesSec,p,&matSize);
2688:       if (!matSize) { /* incoming matrix is identity */
2689:         PetscInt childId;

2691:         childId = childIds[p-pStartF];
2692:         if (childId < 0) { /* no child interpolation: one nnz per */
2693:           if (numFields) {
2694:             PetscInt f;
2695:             for (f = 0; f < numFields; f++) {
2696:               PetscInt numRows = offsets[f+1] - offsets[f], row;
2697:               for (row = 0; row < numRows; row++) {
2698:                 PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2699:                 PetscInt gIndFine   = rowIndices[offsets[f] + row];
2700:                 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2702:                   dnnz[gIndFine - rowStart] = 1;
2703:                 }
2704:                 else if (gIndCoarse >= 0) { /* remote */
2706:                   onnz[gIndFine - rowStart] = 1;
2707:                 }
2708:                 else { /* constrained */
2710:                 }
2711:               }
2712:             }
2713:           }
2714:           else {
2715:             PetscInt i;
2716:             for (i = 0; i < gDof; i++) {
2717:               PetscInt gIndCoarse = pInd[i];
2718:               PetscInt gIndFine   = rowIndices[i];
2719:               if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2721:                 dnnz[gIndFine - rowStart] = 1;
2722:               }
2723:               else if (gIndCoarse >= 0) { /* remote */
2725:                 onnz[gIndFine - rowStart] = 1;
2726:               }
2727:               else { /* constrained */
2729:               }
2730:             }
2731:           }
2732:         }
2733:         else { /* interpolate from all */
2734:           if (numFields) {
2735:             PetscInt f;
2736:             for (f = 0; f < numFields; f++) {
2737:               PetscInt numRows = offsets[f+1] - offsets[f], row;
2738:               for (row = 0; row < numRows; row++) {
2739:                 PetscInt gIndFine = rowIndices[offsets[f] + row];
2740:                 if (gIndFine >= 0) {
2742:                   dnnz[gIndFine - rowStart] = numD[f];
2743:                   onnz[gIndFine - rowStart] = numO[f];
2744:                 }
2745:               }
2746:             }
2747:           }
2748:           else {
2749:             PetscInt i;
2750:             for (i = 0; i < gDof; i++) {
2751:               PetscInt gIndFine = rowIndices[i];
2752:               if (gIndFine >= 0) {
2754:                 dnnz[gIndFine - rowStart] = numD[0];
2755:                 onnz[gIndFine - rowStart] = numO[0];
2756:               }
2757:             }
2758:           }
2759:         }
2760:       }
2761:       else { /* interpolate from all */
2762:         if (numFields) {
2763:           PetscInt f;
2764:           for (f = 0; f < numFields; f++) {
2765:             PetscInt numRows = offsets[f+1] - offsets[f], row;
2766:             for (row = 0; row < numRows; row++) {
2767:               PetscInt gIndFine = rowIndices[offsets[f] + row];
2768:               if (gIndFine >= 0) {
2770:                 dnnz[gIndFine - rowStart] = numD[f];
2771:                 onnz[gIndFine - rowStart] = numO[f];
2772:               }
2773:             }
2774:           }
2775:         }
2776:         else { /* every dof get a full row */
2777:           PetscInt i;
2778:           for (i = 0; i < gDof; i++) {
2779:             PetscInt gIndFine = rowIndices[i];
2780:             if (gIndFine >= 0) {
2782:               dnnz[gIndFine - rowStart] = numD[0];
2783:               onnz[gIndFine - rowStart] = numO[0];
2784:             }
2785:           }
2786:         }
2787:       }
2788:     }
2789:     MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);
2790:     PetscFree2(dnnz,onnz);

2792:     DMPlexGetReferenceTree(fine,&refTree);
2793:     DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
2794:     DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL);
2795:     DMPlexGetAnchors(refTree,&refAnSec,NULL);
2796:     PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
2797:     PetscSectionGetMaxDof(refConSec,&maxConDof);
2798:     PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);
2799:     PetscMalloc1(maxConDof*maxColumns,&pointWork);
2800:     for (p = leafStart; p < leafEnd; p++) {
2801:       PetscInt gDof, gcDof, gOff;
2802:       PetscInt numColIndices, pIndOff, *pInd;
2803:       PetscInt matSize;
2804:       PetscInt childId;

2806:       PetscSectionGetDof(globalFine,p,&gDof);
2807:       PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2808:       if ((gDof - gcDof) <= 0) {
2809:         continue;
2810:       }
2811:       childId = childIds[p-pStartF];
2812:       PetscSectionGetOffset(globalFine,p,&gOff);
2813:       PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2814:       PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2815:       numColIndices -= 2 * numFields;
2816:       pInd = &leafIndices[pIndOff];
2817:       offsets[0]        = 0;
2818:       offsetsCopy[0]    = 0;
2819:       newOffsets[0]     = 0;
2820:       newOffsetsCopy[0] = 0;
2821:       rowOffsets[0]     = 0;
2822:       if (numFields) {
2823:         PetscInt f;
2824:         for (f = 0; f < numFields; f++) {
2825:           PetscInt rowDof;

2827:           PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2828:           offsets[f + 1]     = offsets[f] + rowDof;
2829:           offsetsCopy[f + 1] = offsets[f + 1];
2830:           rowOffsets[f + 1]  = pInd[numColIndices + f];
2831:           newOffsets[f + 1]  = pInd[numColIndices + numFields + f];
2832:         }
2833:         DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,rowIndices);
2834:       }
2835:       else {
2836:         DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,rowIndices);
2837:       }
2838:       PetscSectionGetDof(leafMatricesSec,p,&matSize);
2839:       if (!matSize) { /* incoming matrix is identity */
2840:         if (childId < 0) { /* no child interpolation: scatter */
2841:           if (numFields) {
2842:             PetscInt f;
2843:             for (f = 0; f < numFields; f++) {
2844:               PetscInt numRows = offsets[f+1] - offsets[f], row;
2845:               for (row = 0; row < numRows; row++) {
2846:                 MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);
2847:               }
2848:             }
2849:           }
2850:           else {
2851:             PetscInt numRows = gDof, row;
2852:             for (row = 0; row < numRows; row++) {
2853:               MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);
2854:             }
2855:           }
2856:         }
2857:         else { /* interpolate from all */
2858:           if (numFields) {
2859:             PetscInt f;
2860:             for (f = 0; f < numFields; f++) {
2861:               PetscInt numRows = offsets[f+1] - offsets[f];
2862:               PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2863:               MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);
2864:             }
2865:           }
2866:           else {
2867:             MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);
2868:           }
2869:         }
2870:       }
2871:       else { /* interpolate from all */
2872:         PetscInt    pMatOff;
2873:         PetscScalar *pMat;

2875:         PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);
2876:         pMat = &leafMatrices[pMatOff];
2877:         if (childId < 0) { /* copy the incoming matrix */
2878:           if (numFields) {
2879:             PetscInt f, count;
2880:             for (f = 0, count = 0; f < numFields; f++) {
2881:               PetscInt numRows = offsets[f+1]-offsets[f];
2882:               PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2883:               PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2884:               PetscScalar *inMat = &pMat[count];

2886:               MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);
2887:               count += numCols * numInRows;
2888:             }
2889:           }
2890:           else {
2891:             MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);
2892:           }
2893:         }
2894:         else { /* multiply the incoming matrix by the child interpolation */
2895:           if (numFields) {
2896:             PetscInt f, count;
2897:             for (f = 0, count = 0; f < numFields; f++) {
2898:               PetscInt numRows = offsets[f+1]-offsets[f];
2899:               PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2900:               PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2901:               PetscScalar *inMat = &pMat[count];
2902:               PetscInt i, j, k;
2904:               for (i = 0; i < numRows; i++) {
2905:                 for (j = 0; j < numCols; j++) {
2906:                   PetscScalar val = 0.;
2907:                   for (k = 0; k < numInRows; k++) {
2908:                     val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2909:                   }
2910:                   pointWork[i * numCols + j] = val;
2911:                 }
2912:               }
2913:               MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);
2914:               count += numCols * numInRows;
2915:             }
2916:           }
2917:           else { /* every dof gets a full row */
2918:             PetscInt numRows   = gDof;
2919:             PetscInt numCols   = numColIndices;
2920:             PetscInt numInRows = matSize / numColIndices;
2921:             PetscInt i, j, k;
2923:             for (i = 0; i < numRows; i++) {
2924:               for (j = 0; j < numCols; j++) {
2925:                 PetscScalar val = 0.;
2926:                 for (k = 0; k < numInRows; k++) {
2927:                   val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2928:                 }
2929:                 pointWork[i * numCols + j] = val;
2930:               }
2931:             }
2932:             MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);
2933:           }
2934:         }
2935:       }
2936:     }
2937:     DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
2938:     DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
2939:     PetscFree(pointWork);
2940:   }
2941:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
2942:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
2943:   PetscSectionDestroy(&leafIndicesSec);
2944:   PetscSectionDestroy(&leafMatricesSec);
2945:   PetscFree2(leafIndices,leafMatrices);
2946:   PetscFree2(*(PetscInt****)&perms,*(PetscScalar****)&flips);
2947:   PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);
2948:   ISRestoreIndices(aIS,&anchors);
2949:   return 0;
2950: }

2952: /*
2953:  * Assuming a nodal basis (w.r.t. the dual basis) basis:
2954:  *
2955:  * for each coarse dof \phi^c_i:
2956:  *   for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
2957:  *     for each fine dof \phi^f_j;
2958:  *       a_{i,j} = 0;
2959:  *       for each fine dof \phi^f_k:
2960:  *         a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
2961:  *                    [^^^ this is = \phi^c_i ^^^]
2962:  */
2963: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
2964: {
2965:   PetscDS        ds;
2966:   PetscSection   section, cSection;
2967:   DMLabel        canonical, depth;
2968:   Mat            cMat, mat;
2969:   PetscInt       *nnz;
2970:   PetscInt       f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
2971:   PetscInt       m, n;
2972:   PetscScalar    *pointScalar;
2973:   PetscReal      *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;

2975:   DMGetLocalSection(refTree,&section);
2976:   DMGetDimension(refTree, &dim);
2977:   PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);
2978:   PetscMalloc2(dim,&pointScalar,dim,&pointRef);
2979:   DMGetDS(refTree,&ds);
2980:   PetscDSGetNumFields(ds,&numFields);
2981:   PetscSectionGetNumFields(section,&numSecFields);
2982:   DMGetLabel(refTree,"canonical",&canonical);
2983:   DMGetLabel(refTree,"depth",&depth);
2984:   DMGetDefaultConstraints(refTree,&cSection,&cMat,NULL);
2985:   DMPlexGetChart(refTree, &pStart, &pEnd);
2986:   DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);
2987:   MatGetSize(cMat,&n,&m); /* the injector has transpose sizes from the constraint matrix */
2988:   /* Step 1: compute non-zero pattern.  A proper subset of constraint matrix non-zero */
2989:   PetscCalloc1(m,&nnz);
2990:   for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
2991:     const PetscInt *children;
2992:     PetscInt numChildren;
2993:     PetscInt i, numChildDof, numSelfDof;

2995:     if (canonical) {
2996:       PetscInt pCanonical;
2997:       DMLabelGetValue(canonical,p,&pCanonical);
2998:       if (p != pCanonical) continue;
2999:     }
3000:     DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3001:     if (!numChildren) continue;
3002:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
3003:       PetscInt child = children[i];
3004:       PetscInt dof;

3006:       PetscSectionGetDof(section,child,&dof);
3007:       numChildDof += dof;
3008:     }
3009:     PetscSectionGetDof(section,p,&numSelfDof);
3010:     if (!numChildDof || !numSelfDof) continue;
3011:     for (f = 0; f < numFields; f++) {
3012:       PetscInt selfOff;

3014:       if (numSecFields) { /* count the dofs for just this field */
3015:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
3016:           PetscInt child = children[i];
3017:           PetscInt dof;

3019:           PetscSectionGetFieldDof(section,child,f,&dof);
3020:           numChildDof += dof;
3021:         }
3022:         PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3023:         PetscSectionGetFieldOffset(section,p,f,&selfOff);
3024:       }
3025:       else {
3026:         PetscSectionGetOffset(section,p,&selfOff);
3027:       }
3028:       for (i = 0; i < numSelfDof; i++) {
3029:         nnz[selfOff + i] = numChildDof;
3030:       }
3031:     }
3032:   }
3033:   MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);
3034:   PetscFree(nnz);
3035:   /* Setp 2: compute entries */
3036:   for (p = pStart; p < pEnd; p++) {
3037:     const PetscInt *children;
3038:     PetscInt numChildren;
3039:     PetscInt i, numChildDof, numSelfDof;

3041:     /* same conditions about when entries occur */
3042:     if (canonical) {
3043:       PetscInt pCanonical;
3044:       DMLabelGetValue(canonical,p,&pCanonical);
3045:       if (p != pCanonical) continue;
3046:     }
3047:     DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3048:     if (!numChildren) continue;
3049:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
3050:       PetscInt child = children[i];
3051:       PetscInt dof;

3053:       PetscSectionGetDof(section,child,&dof);
3054:       numChildDof += dof;
3055:     }
3056:     PetscSectionGetDof(section,p,&numSelfDof);
3057:     if (!numChildDof || !numSelfDof) continue;

3059:     for (f = 0; f < numFields; f++) {
3060:       PetscInt       pI = -1, cI = -1;
3061:       PetscInt       selfOff, Nc, parentCell;
3062:       PetscInt       cellShapeOff;
3063:       PetscObject    disc;
3064:       PetscDualSpace dsp;
3065:       PetscClassId   classId;
3066:       PetscScalar    *pointMat;
3067:       PetscInt       *matRows, *matCols;
3068:       PetscInt       pO = PETSC_MIN_INT;
3069:       const PetscInt *depthNumDof;

3071:       if (numSecFields) {
3072:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
3073:           PetscInt child = children[i];
3074:           PetscInt dof;

3076:           PetscSectionGetFieldDof(section,child,f,&dof);
3077:           numChildDof += dof;
3078:         }
3079:         PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3080:         PetscSectionGetFieldOffset(section,p,f,&selfOff);
3081:       }
3082:       else {
3083:         PetscSectionGetOffset(section,p,&selfOff);
3084:       }

3086:       /* find a cell whose closure contains p */
3087:       if (p >= cStart && p < cEnd) {
3088:         parentCell = p;
3089:       }
3090:       else {
3091:         PetscInt *star = NULL;
3092:         PetscInt numStar;

3094:         parentCell = -1;
3095:         DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3096:         for (i = numStar - 1; i >= 0; i--) {
3097:           PetscInt c = star[2 * i];

3099:           if (c >= cStart && c < cEnd) {
3100:             parentCell = c;
3101:             break;
3102:           }
3103:         }
3104:         DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3105:       }
3106:       /* determine the offset of p's shape functions within parentCell's shape functions */
3107:       PetscDSGetDiscretization(ds,f,&disc);
3108:       PetscObjectGetClassId(disc,&classId);
3109:       if (classId == PETSCFE_CLASSID) {
3110:         PetscFEGetDualSpace((PetscFE)disc,&dsp);
3111:       }
3112:       else if (classId == PETSCFV_CLASSID) {
3113:         PetscFVGetDualSpace((PetscFV)disc,&dsp);
3114:       }
3115:       else {
3116:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object");
3117:       }
3118:       PetscDualSpaceGetNumDof(dsp,&depthNumDof);
3119:       PetscDualSpaceGetNumComponents(dsp,&Nc);
3120:       {
3121:         PetscInt *closure = NULL;
3122:         PetscInt numClosure;

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

3128:           pO = closure[2 * i + 1];
3129:           if (point == p) {
3130:             pI = i;
3131:             break;
3132:           }
3133:           DMLabelGetValue(depth,point,&pointDepth);
3134:           cellShapeOff += depthNumDof[pointDepth];
3135:         }
3136:         DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);
3137:       }

3139:       DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);
3140:       DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);
3141:       matCols = matRows + numSelfDof;
3142:       for (i = 0; i < numSelfDof; i++) {
3143:         matRows[i] = selfOff + i;
3144:       }
3145:       for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
3146:       {
3147:         PetscInt colOff = 0;

3149:         for (i = 0; i < numChildren; i++) {
3150:           PetscInt child = children[i];
3151:           PetscInt dof, off, j;

3153:           if (numSecFields) {
3154:             PetscSectionGetFieldDof(cSection,child,f,&dof);
3155:             PetscSectionGetFieldOffset(cSection,child,f,&off);
3156:           }
3157:           else {
3158:             PetscSectionGetDof(cSection,child,&dof);
3159:             PetscSectionGetOffset(cSection,child,&off);
3160:           }

3162:           for (j = 0; j < dof; j++) {
3163:             matCols[colOff++] = off + j;
3164:           }
3165:         }
3166:       }
3167:       if (classId == PETSCFE_CLASSID) {
3168:         PetscFE        fe = (PetscFE) disc;
3169:         PetscInt       fSize;
3170:         const PetscInt ***perms;
3171:         const PetscScalar ***flips;
3172:         const PetscInt *pperms;

3174:         PetscFEGetDualSpace(fe,&dsp);
3175:         PetscDualSpaceGetDimension(dsp,&fSize);
3176:         PetscDualSpaceGetSymmetries(dsp, &perms, &flips);
3177:         pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL;
3178:         for (i = 0; i < numSelfDof; i++) { /* for every shape function */
3179:           PetscQuadrature q;
3180:           PetscInt        dim, thisNc, numPoints, j, k;
3181:           const PetscReal *points;
3182:           const PetscReal *weights;
3183:           PetscInt        *closure = NULL;
3184:           PetscInt        numClosure;
3185:           PetscInt        iCell = pperms ? pperms[i] : i;
3186:           PetscInt        parentCellShapeDof = cellShapeOff + iCell;
3187:           PetscTabulation Tparent;

3189:           PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);
3190:           PetscQuadratureGetData(q,&dim,&thisNc,&numPoints,&points,&weights);
3192:           PetscFECreateTabulation(fe,1,numPoints,points,0,&Tparent); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3193:           for (j = 0; j < numPoints; j++) {
3194:             PetscInt          childCell = -1;
3195:             PetscReal         *parentValAtPoint;
3196:             const PetscReal   xi0[3] = {-1.,-1.,-1.};
3197:             const PetscReal   *pointReal = &points[dim * j];
3198:             const PetscScalar *point;
3199:             PetscTabulation Tchild;
3200:             PetscInt          childCellShapeOff, pointMatOff;
3201: #if defined(PETSC_USE_COMPLEX)
3202:             PetscInt          d;

3204:             for (d = 0; d < dim; d++) {
3205:               pointScalar[d] = points[dim * j + d];
3206:             }
3207:             point = pointScalar;
3208: #else
3209:             point = pointReal;
3210: #endif

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

3214:             for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3215:               PetscInt child = children[k];
3216:               PetscInt *star = NULL;
3217:               PetscInt numStar, s;

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

3223:                 if (c < cStart || c >= cEnd) continue;
3224:                 DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);
3225:                 if (childCell >= 0) break;
3226:               }
3227:               DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);
3228:               if (childCell >= 0) break;
3229:             }
3231:             DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);
3232:             DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);
3233:             CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
3234:             CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);

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

3243:               DMLabelGetValue(depth,child,&childDepth);
3244:               childDof = depthNumDof[childDepth];
3245:               for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) {
3246:                 PetscInt point = closure[2 * l];
3247:                 PetscInt pointDepth;

3249:                 childO = closure[2 * l + 1];
3250:                 if (point == child) {
3251:                   cI = l;
3252:                   break;
3253:                 }
3254:                 DMLabelGetValue(depth,point,&pointDepth);
3255:                 childCellShapeOff += depthNumDof[pointDepth];
3256:               }
3257:               if (l == numClosure) {
3258:                 pointMatOff += childDof;
3259:                 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3260:               }
3261:               cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL;
3262:               for (l = 0; l < childDof; l++) {
3263:                 PetscInt    lCell = cperms ? cperms[l] : l;
3264:                 PetscInt    childCellDof = childCellShapeOff + lCell;
3265:                 PetscReal   *childValAtPoint;
3266:                 PetscReal   val = 0.;

3268:                 childValAtPoint = &Tchild->T[0][childCellDof * Nc];
3269:                 for (m = 0; m < Nc; m++) {
3270:                   val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];
3271:                 }

3273:                 pointMat[i * numChildDof + pointMatOff + l] += val;
3274:               }
3275:               pointMatOff += childDof;
3276:             }
3277:             DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3278:             PetscTabulationDestroy(&Tchild);
3279:           }
3280:           PetscTabulationDestroy(&Tparent);
3281:         }
3282:       }
3283:       else { /* just the volume-weighted averages of the children */
3284:         PetscReal parentVol;
3285:         PetscInt  childCell;

3287:         DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);
3288:         for (i = 0, childCell = 0; i < numChildren; i++) {
3289:           PetscInt  child = children[i], j;
3290:           PetscReal childVol;

3292:           if (child < cStart || child >= cEnd) continue;
3293:           DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);
3294:           for (j = 0; j < Nc; j++) {
3295:             pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3296:           }
3297:           childCell++;
3298:         }
3299:       }
3300:       /* Insert pointMat into mat */
3301:       MatSetValues(mat,numSelfDof,matRows,numChildDof,matCols,pointMat,INSERT_VALUES);
3302:       DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);
3303:       DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);
3304:     }
3305:   }
3306:   PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);
3307:   PetscFree2(pointScalar,pointRef);
3308:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3309:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3310:   *inj = mat;
3311:   return 0;
3312: }

3314: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3315: {
3316:   PetscDS        ds;
3317:   PetscInt       numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3318:   PetscScalar    ***refPointFieldMats;
3319:   PetscSection   refConSec, refSection;

3321:   DMGetDS(refTree,&ds);
3322:   PetscDSGetNumFields(ds,&numFields);
3323:   DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL);
3324:   DMGetLocalSection(refTree,&refSection);
3325:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3326:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
3327:   PetscSectionGetMaxDof(refConSec,&maxDof);
3328:   PetscMalloc1(maxDof,&rows);
3329:   PetscMalloc1(maxDof*maxDof,&cols);
3330:   for (p = pRefStart; p < pRefEnd; p++) {
3331:     PetscInt parent, pDof, parentDof;

3333:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
3334:     PetscSectionGetDof(refConSec,p,&pDof);
3335:     PetscSectionGetDof(refSection,parent,&parentDof);
3336:     if (!pDof || !parentDof || parent == p) continue;

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

3342:       if (numFields > 1) {
3343:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3344:         PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
3345:       }
3346:       else {
3347:         PetscSectionGetDof(refConSec,p,&cDof);
3348:         PetscSectionGetOffset(refConSec,p,&cOff);
3349:       }

3351:       for (r = 0; r < cDof; r++) {
3352:         rows[r] = cOff + r;
3353:       }
3354:       numCols = 0;
3355:       {
3356:         PetscInt aDof, aOff, j;

3358:         if (numFields > 1) {
3359:           PetscSectionGetFieldDof(refSection,parent,f,&aDof);
3360:           PetscSectionGetFieldOffset(refSection,parent,f,&aOff);
3361:         }
3362:         else {
3363:           PetscSectionGetDof(refSection,parent,&aDof);
3364:           PetscSectionGetOffset(refSection,parent,&aOff);
3365:         }

3367:         for (j = 0; j < aDof; j++) {
3368:           cols[numCols++] = aOff + j;
3369:         }
3370:       }
3371:       PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
3372:       /* transpose of constraint matrix */
3373:       MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);
3374:     }
3375:   }
3376:   *childrenMats = refPointFieldMats;
3377:   PetscFree(rows);
3378:   PetscFree(cols);
3379:   return 0;
3380: }

3382: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3383: {
3384:   PetscDS        ds;
3385:   PetscScalar    ***refPointFieldMats;
3386:   PetscInt       numFields, pRefStart, pRefEnd, p, f;
3387:   PetscSection   refConSec, refSection;

3389:   refPointFieldMats = *childrenMats;
3390:   *childrenMats = NULL;
3391:   DMGetDS(refTree,&ds);
3392:   DMGetLocalSection(refTree,&refSection);
3393:   PetscDSGetNumFields(ds,&numFields);
3394:   DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL);
3395:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3396:   for (p = pRefStart; p < pRefEnd; p++) {
3397:     PetscInt parent, pDof, parentDof;

3399:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
3400:     PetscSectionGetDof(refConSec,p,&pDof);
3401:     PetscSectionGetDof(refSection,parent,&parentDof);
3402:     if (!pDof || !parentDof || parent == p) continue;

3404:     for (f = 0; f < numFields; f++) {
3405:       PetscInt cDof;

3407:       if (numFields > 1) {
3408:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3409:       }
3410:       else {
3411:         PetscSectionGetDof(refConSec,p,&cDof);
3412:       }

3414:       PetscFree(refPointFieldMats[p - pRefStart][f]);
3415:     }
3416:     PetscFree(refPointFieldMats[p - pRefStart]);
3417:   }
3418:   PetscFree(refPointFieldMats);
3419:   return 0;
3420: }

3422: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef)
3423: {
3424:   Mat            cMatRef;
3425:   PetscObject    injRefObj;

3427:   DMGetDefaultConstraints(refTree,NULL,&cMatRef,NULL);
3428:   PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);
3429:   *injRef = (Mat) injRefObj;
3430:   if (!*injRef) {
3431:     DMPlexComputeInjectorReferenceTree(refTree,injRef);
3432:     PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);
3433:     /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3434:     PetscObjectDereference((PetscObject)*injRef);
3435:   }
3436:   return 0;
3437: }

3439: 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)
3440: {
3441:   PetscInt       pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3442:   PetscSection   globalCoarse, globalFine;
3443:   PetscSection   localCoarse, localFine, leafIndicesSec;
3444:   PetscSection   multiRootSec, rootIndicesSec;
3445:   PetscInt       *leafInds, *rootInds = NULL;
3446:   const PetscInt *rootDegrees;
3447:   PetscScalar    *leafVals = NULL, *rootVals = NULL;
3448:   PetscSF        coarseToFineEmbedded;

3450:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3451:   DMPlexGetChart(fine,&pStartF,&pEndF);
3452:   DMGetLocalSection(fine,&localFine);
3453:   DMGetGlobalSection(fine,&globalFine);
3454:   PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
3455:   PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);
3456:   PetscSectionGetMaxDof(localFine,&maxDof);
3457:   { /* winnow fine points that don't have global dofs out of the sf */
3458:     PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3459:     const PetscInt *leaves;

3461:     PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);
3462:     for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3463:       p    = leaves ? leaves[l] : l;
3464:       PetscSectionGetDof(globalFine,p,&dof);
3465:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3466:       if ((dof - cdof) > 0) {
3467:         numPointsWithDofs++;

3469:         PetscSectionGetDof(localFine,p,&dof);
3470:         PetscSectionSetDof(leafIndicesSec,p,dof + 1);
3471:       }
3472:     }
3473:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3474:     PetscSectionSetUp(leafIndicesSec);
3475:     PetscSectionGetStorageSize(leafIndicesSec,&numIndices);
3476:     PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);
3477:     if (gatheredValues)  PetscMalloc1(numIndices,&leafVals);
3478:     for (l = 0, offset = 0; l < nleaves; l++) {
3479:       p    = leaves ? leaves[l] : l;
3480:       PetscSectionGetDof(globalFine,p,&dof);
3481:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3482:       if ((dof - cdof) > 0) {
3483:         PetscInt    off, gOff;
3484:         PetscInt    *pInd;
3485:         PetscScalar *pVal = NULL;

3487:         pointsWithDofs[offset++] = l;

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

3491:         pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3492:         if (gatheredValues) {
3493:           PetscInt i;

3495:           pVal = &leafVals[off + 1];
3496:           for (i = 0; i < dof; i++) pVal[i] = 0.;
3497:         }
3498:         PetscSectionGetOffset(globalFine,p,&gOff);

3500:         offsets[0] = 0;
3501:         if (numFields) {
3502:           PetscInt f;

3504:           for (f = 0; f < numFields; f++) {
3505:             PetscInt fDof;
3506:             PetscSectionGetFieldDof(localFine,p,f,&fDof);
3507:             offsets[f + 1] = fDof + offsets[f];
3508:           }
3509:           DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1, NULL,pInd);
3510:         } else {
3511:           DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL, NULL,pInd);
3512:         }
3513:         if (gatheredValues) VecGetValues(fineVec,dof,pInd,pVal);
3514:       }
3515:     }
3516:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
3517:     PetscFree(pointsWithDofs);
3518:   }

3520:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3521:   DMGetLocalSection(coarse,&localCoarse);
3522:   DMGetGlobalSection(coarse,&globalCoarse);

3524:   { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3525:     MPI_Datatype threeInt;
3526:     PetscMPIInt  rank;
3527:     PetscInt     (*parentNodeAndIdCoarse)[3];
3528:     PetscInt     (*parentNodeAndIdFine)[3];
3529:     PetscInt     p, nleaves, nleavesToParents;
3530:     PetscSF      pointSF, sfToParents;
3531:     const PetscInt *ilocal;
3532:     const PetscSFNode *iremote;
3533:     PetscSFNode  *iremoteToParents;
3534:     PetscInt     *ilocalToParents;

3536:     MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
3537:     MPI_Type_contiguous(3,MPIU_INT,&threeInt);
3538:     MPI_Type_commit(&threeInt);
3539:     PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);
3540:     DMGetPointSF(coarse,&pointSF);
3541:     PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);
3542:     for (p = pStartC; p < pEndC; p++) {
3543:       PetscInt parent, childId;
3544:       DMPlexGetTreeParent(coarse,p,&parent,&childId);
3545:       parentNodeAndIdCoarse[p - pStartC][0] = rank;
3546:       parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3547:       parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3548:       if (nleaves > 0) {
3549:         PetscInt leaf = -1;

3551:         if (ilocal) {
3552:           PetscFindInt(parent,nleaves,ilocal,&leaf);
3553:         }
3554:         else {
3555:           leaf = p - pStartC;
3556:         }
3557:         if (leaf >= 0) {
3558:           parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3559:           parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3560:         }
3561:       }
3562:     }
3563:     for (p = pStartF; p < pEndF; p++) {
3564:       parentNodeAndIdFine[p - pStartF][0] = -1;
3565:       parentNodeAndIdFine[p - pStartF][1] = -1;
3566:       parentNodeAndIdFine[p - pStartF][2] = -1;
3567:     }
3568:     PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine,MPI_REPLACE);
3569:     PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine,MPI_REPLACE);
3570:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3571:       PetscInt dof;

3573:       PetscSectionGetDof(leafIndicesSec,p,&dof);
3574:       if (dof) {
3575:         PetscInt off;

3577:         PetscSectionGetOffset(leafIndicesSec,p,&off);
3578:         if (gatheredIndices) {
3579:           leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3580:         } else if (gatheredValues) {
3581:           leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3582:         }
3583:       }
3584:       if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3585:         nleavesToParents++;
3586:       }
3587:     }
3588:     PetscMalloc1(nleavesToParents,&ilocalToParents);
3589:     PetscMalloc1(nleavesToParents,&iremoteToParents);
3590:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3591:       if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3592:         ilocalToParents[nleavesToParents] = p - pStartF;
3593:         iremoteToParents[nleavesToParents].rank  = parentNodeAndIdFine[p-pStartF][0];
3594:         iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1];
3595:         nleavesToParents++;
3596:       }
3597:     }
3598:     PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);
3599:     PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);
3600:     PetscSFDestroy(&coarseToFineEmbedded);

3602:     coarseToFineEmbedded = sfToParents;

3604:     PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);
3605:     MPI_Type_free(&threeInt);
3606:   }

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

3612:     for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3613:       PetscSectionGetDof(globalCoarse,p,&dof);
3614:       PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3615:       if ((dof - cdof) > 0) {
3616:         numPointsWithDofs++;
3617:       }
3618:     }
3619:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3620:     for (p = pStartC, offset = 0; p < pEndC; p++) {
3621:       PetscSectionGetDof(globalCoarse,p,&dof);
3622:       PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3623:       if ((dof - cdof) > 0) {
3624:         pointsWithDofs[offset++] = p - pStartC;
3625:       }
3626:     }
3627:     PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);
3628:     PetscSFDestroy(&coarseToFineEmbedded);
3629:     PetscFree(pointsWithDofs);
3630:     coarseToFineEmbedded = sfDofsOnly;
3631:   }

3633:   /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3634:   PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);
3635:   PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);
3636:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);
3637:   PetscSectionSetChart(multiRootSec,pStartC,pEndC);
3638:   for (p = pStartC; p < pEndC; p++) {
3639:     PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);
3640:   }
3641:   PetscSectionSetUp(multiRootSec);
3642:   PetscSectionGetStorageSize(multiRootSec,&numMulti);
3643:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
3644:   { /* distribute the leaf section */
3645:     PetscSF multi, multiInv, indicesSF;
3646:     PetscInt *remoteOffsets, numRootIndices;

3648:     PetscSFGetMultiSF(coarseToFineEmbedded,&multi);
3649:     PetscSFCreateInverseSF(multi,&multiInv);
3650:     PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);
3651:     PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);
3652:     PetscFree(remoteOffsets);
3653:     PetscSFDestroy(&multiInv);
3654:     PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
3655:     if (gatheredIndices) {
3656:       PetscMalloc1(numRootIndices,&rootInds);
3657:       PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds,MPI_REPLACE);
3658:       PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds,MPI_REPLACE);
3659:     }
3660:     if (gatheredValues) {
3661:       PetscMalloc1(numRootIndices,&rootVals);
3662:       PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals,MPI_REPLACE);
3663:       PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals,MPI_REPLACE);
3664:     }
3665:     PetscSFDestroy(&indicesSF);
3666:   }
3667:   PetscSectionDestroy(&leafIndicesSec);
3668:   PetscFree(leafInds);
3669:   PetscFree(leafVals);
3670:   PetscSFDestroy(&coarseToFineEmbedded);
3671:   *rootMultiSec = multiRootSec;
3672:   *multiLeafSec = rootIndicesSec;
3673:   if (gatheredIndices) *gatheredIndices = rootInds;
3674:   if (gatheredValues)  *gatheredValues  = rootVals;
3675:   return 0;
3676: }

3678: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3679: {
3680:   DM             refTree;
3681:   PetscSection   multiRootSec, rootIndicesSec;
3682:   PetscSection   globalCoarse, globalFine;
3683:   PetscSection   localCoarse, localFine;
3684:   PetscSection   cSecRef;
3685:   PetscInt       *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd;
3686:   Mat            injRef;
3687:   PetscInt       numFields, maxDof;
3688:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
3689:   PetscInt       *offsets, *offsetsCopy, *rowOffsets;
3690:   PetscLayout    rowMap, colMap;
3691:   PetscInt       rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3692:   PetscScalar    ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */


3695:   /* get the templates for the fine-to-coarse injection from the reference tree */
3696:   DMPlexGetReferenceTree(coarse,&refTree);
3697:   DMGetDefaultConstraints(refTree,&cSecRef,NULL,NULL);
3698:   PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
3699:   DMPlexReferenceTreeGetInjector(refTree,&injRef);

3701:   DMPlexGetChart(fine,&pStartF,&pEndF);
3702:   DMGetLocalSection(fine,&localFine);
3703:   DMGetGlobalSection(fine,&globalFine);
3704:   PetscSectionGetNumFields(localFine,&numFields);
3705:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3706:   DMGetLocalSection(coarse,&localCoarse);
3707:   DMGetGlobalSection(coarse,&globalCoarse);
3708:   PetscSectionGetMaxDof(localCoarse,&maxDof);
3709:   {
3710:     PetscInt maxFields = PetscMax(1,numFields) + 1;
3711:     PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
3712:   }

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

3716:   PetscMalloc1(maxDof,&parentIndices);

3718:   /* count indices */
3719:   MatGetLayouts(mat,&rowMap,&colMap);
3720:   PetscLayoutSetUp(rowMap);
3721:   PetscLayoutSetUp(colMap);
3722:   PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
3723:   PetscLayoutGetRange(colMap,&colStart,&colEnd);
3724:   PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);
3725:   for (p = pStartC; p < pEndC; p++) {
3726:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3728:     PetscSectionGetDof(globalCoarse,p,&dof);
3729:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3730:     if ((dof - cdof) <= 0) continue;
3731:     PetscSectionGetOffset(globalCoarse,p,&gOff);

3733:     rowOffsets[0] = 0;
3734:     offsetsCopy[0] = 0;
3735:     if (numFields) {
3736:       PetscInt f;

3738:       for (f = 0; f < numFields; f++) {
3739:         PetscInt fDof;
3740:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3741:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3742:       }
3743:       DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,parentIndices);
3744:     } else {
3745:       DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,parentIndices);
3746:       rowOffsets[1] = offsetsCopy[0];
3747:     }

3749:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
3750:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
3751:     leafEnd = leafStart + numLeaves;
3752:     for (l = leafStart; l < leafEnd; l++) {
3753:       PetscInt numIndices, childId, offset;
3754:       const PetscInt *childIndices;

3756:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3757:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
3758:       childId = rootIndices[offset++];
3759:       childIndices = &rootIndices[offset];
3760:       numIndices--;

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

3765:         for (i = 0; i < numIndices; i++) {
3766:           PetscInt colIndex = childIndices[i];
3767:           PetscInt rowIndex = parentIndices[i];
3768:           if (rowIndex < 0) continue;
3770:           if (colIndex >= colStart && colIndex < colEnd) {
3771:             nnzD[rowIndex - rowStart] = 1;
3772:           }
3773:           else {
3774:             nnzO[rowIndex - rowStart] = 1;
3775:           }
3776:         }
3777:       }
3778:       else {
3779:         PetscInt parentId, f, lim;

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

3783:         lim = PetscMax(1,numFields);
3784:         offsets[0] = 0;
3785:         if (numFields) {
3786:           PetscInt f;

3788:           for (f = 0; f < numFields; f++) {
3789:             PetscInt fDof;
3790:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

3792:             offsets[f + 1] = fDof + offsets[f];
3793:           }
3794:         }
3795:         else {
3796:           PetscInt cDof;

3798:           PetscSectionGetDof(cSecRef,childId,&cDof);
3799:           offsets[1] = cDof;
3800:         }
3801:         for (f = 0; f < lim; f++) {
3802:           PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3803:           PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3804:           PetscInt i, numD = 0, numO = 0;

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

3809:             if (colIndex < 0) continue;
3810:             if (colIndex >= colStart && colIndex < colEnd) {
3811:               numD++;
3812:             }
3813:             else {
3814:               numO++;
3815:             }
3816:           }
3817:           for (i = parentStart; i < parentEnd; i++) {
3818:             PetscInt rowIndex = parentIndices[i];

3820:             if (rowIndex < 0) continue;
3821:             nnzD[rowIndex - rowStart] += numD;
3822:             nnzO[rowIndex - rowStart] += numO;
3823:           }
3824:         }
3825:       }
3826:     }
3827:   }
3828:   /* preallocate */
3829:   MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);
3830:   PetscFree2(nnzD,nnzO);
3831:   /* insert values */
3832:   DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3833:   for (p = pStartC; p < pEndC; p++) {
3834:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3836:     PetscSectionGetDof(globalCoarse,p,&dof);
3837:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3838:     if ((dof - cdof) <= 0) continue;
3839:     PetscSectionGetOffset(globalCoarse,p,&gOff);

3841:     rowOffsets[0] = 0;
3842:     offsetsCopy[0] = 0;
3843:     if (numFields) {
3844:       PetscInt f;

3846:       for (f = 0; f < numFields; f++) {
3847:         PetscInt fDof;
3848:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3849:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3850:       }
3851:       DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,parentIndices);
3852:     } else {
3853:       DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,parentIndices);
3854:       rowOffsets[1] = offsetsCopy[0];
3855:     }

3857:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
3858:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
3859:     leafEnd = leafStart + numLeaves;
3860:     for (l = leafStart; l < leafEnd; l++) {
3861:       PetscInt numIndices, childId, offset;
3862:       const PetscInt *childIndices;

3864:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3865:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
3866:       childId = rootIndices[offset++];
3867:       childIndices = &rootIndices[offset];
3868:       numIndices--;

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

3873:         for (i = 0; i < numIndices; i++) {
3874:           MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);
3875:         }
3876:       }
3877:       else {
3878:         PetscInt parentId, f, lim;

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

3882:         lim = PetscMax(1,numFields);
3883:         offsets[0] = 0;
3884:         if (numFields) {
3885:           PetscInt f;

3887:           for (f = 0; f < numFields; f++) {
3888:             PetscInt fDof;
3889:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

3891:             offsets[f + 1] = fDof + offsets[f];
3892:           }
3893:         }
3894:         else {
3895:           PetscInt cDof;

3897:           PetscSectionGetDof(cSecRef,childId,&cDof);
3898:           offsets[1] = cDof;
3899:         }
3900:         for (f = 0; f < lim; f++) {
3901:           PetscScalar    *childMat   = &childrenMats[childId - pRefStart][f][0];
3902:           PetscInt       *rowIndices = &parentIndices[rowOffsets[f]];
3903:           const PetscInt *colIndices = &childIndices[offsets[f]];

3905:           MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);
3906:         }
3907:       }
3908:     }
3909:   }
3910:   PetscSectionDestroy(&multiRootSec);
3911:   PetscSectionDestroy(&rootIndicesSec);
3912:   PetscFree(parentIndices);
3913:   DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3914:   PetscFree(rootIndices);
3915:   PetscFree3(offsets,offsetsCopy,rowOffsets);

3917:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3918:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3919:   return 0;
3920: }

3922: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3923: {
3924:   PetscSF           coarseToFineEmbedded;
3925:   PetscSection      globalCoarse, globalFine;
3926:   PetscSection      localCoarse, localFine;
3927:   PetscSection      aSec, cSec;
3928:   PetscSection      rootValuesSec;
3929:   PetscSection      leafValuesSec;
3930:   PetscScalar       *rootValues, *leafValues;
3931:   IS                aIS;
3932:   const PetscInt    *anchors;
3933:   Mat               cMat;
3934:   PetscInt          numFields;
3935:   PetscInt          pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd;
3936:   PetscInt          aStart, aEnd, cStart, cEnd;
3937:   PetscInt          *maxChildIds;
3938:   PetscInt          *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3939:   PetscFV           fv = NULL;
3940:   PetscInt          dim, numFVcomps = -1, fvField = -1;
3941:   DM                cellDM = NULL, gradDM = NULL;
3942:   const PetscScalar *cellGeomArray = NULL;
3943:   const PetscScalar *gradArray = NULL;

3945:   VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
3946:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3947:   DMPlexGetSimplexOrBoxCells(coarse,0,&cellStart,&cellEnd);
3948:   DMPlexGetChart(fine,&pStartF,&pEndF);
3949:   DMGetGlobalSection(fine,&globalFine);
3950:   DMGetCoordinateDim(coarse,&dim);
3951:   { /* winnow fine points that don't have global dofs out of the sf */
3952:     PetscInt       nleaves, l;
3953:     const PetscInt *leaves;
3954:     PetscInt       dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;

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

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

3961:       PetscSectionGetDof(globalFine,p,&dof);
3962:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3963:       if ((dof - cdof) > 0) {
3964:         numPointsWithDofs++;
3965:       }
3966:     }
3967:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3968:     for (l = 0, offset = 0; l < nleaves; l++) {
3969:       PetscInt p = leaves ? leaves[l] : l;

3971:       PetscSectionGetDof(globalFine,p,&dof);
3972:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3973:       if ((dof - cdof) > 0) {
3974:         pointsWithDofs[offset++] = l;
3975:       }
3976:     }
3977:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
3978:     PetscFree(pointsWithDofs);
3979:   }
3980:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
3981:   PetscMalloc1(pEndC-pStartC,&maxChildIds);
3982:   for (p = pStartC; p < pEndC; p++) {
3983:     maxChildIds[p - pStartC] = -2;
3984:   }
3985:   PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);
3986:   PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);

3988:   DMGetLocalSection(coarse,&localCoarse);
3989:   DMGetGlobalSection(coarse,&globalCoarse);

3991:   DMPlexGetAnchors(coarse,&aSec,&aIS);
3992:   ISGetIndices(aIS,&anchors);
3993:   PetscSectionGetChart(aSec,&aStart,&aEnd);

3995:   DMGetDefaultConstraints(coarse,&cSec,&cMat,NULL);
3996:   PetscSectionGetChart(cSec,&cStart,&cEnd);

3998:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
3999:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);
4000:   PetscSectionSetChart(rootValuesSec,pStartC,pEndC);
4001:   PetscSectionGetNumFields(localCoarse,&numFields);
4002:   {
4003:     PetscInt maxFields = PetscMax(1,numFields) + 1;
4004:     PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);
4005:   }
4006:   if (grad) {
4007:     PetscInt i;

4009:     VecGetDM(cellGeom,&cellDM);
4010:     VecGetArrayRead(cellGeom,&cellGeomArray);
4011:     VecGetDM(grad,&gradDM);
4012:     VecGetArrayRead(grad,&gradArray);
4013:     for (i = 0; i < PetscMax(1,numFields); i++) {
4014:       PetscObject  obj;
4015:       PetscClassId id;

4017:       DMGetField(coarse, i, NULL, &obj);
4018:       PetscObjectGetClassId(obj,&id);
4019:       if (id == PETSCFV_CLASSID) {
4020:         fv      = (PetscFV) obj;
4021:         PetscFVGetNumComponents(fv,&numFVcomps);
4022:         fvField = i;
4023:         break;
4024:       }
4025:     }
4026:   }

4028:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
4029:     PetscInt dof;
4030:     PetscInt maxChildId     = maxChildIds[p - pStartC];
4031:     PetscInt numValues      = 0;

4033:     PetscSectionGetDof(globalCoarse,p,&dof);
4034:     if (dof < 0) {
4035:       dof = -(dof + 1);
4036:     }
4037:     offsets[0]    = 0;
4038:     newOffsets[0] = 0;
4039:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
4040:       PetscInt *closure = NULL, closureSize, cl;

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

4046:         PetscSectionGetDof(localCoarse,c,&clDof);
4047:         numValues += clDof;
4048:       }
4049:       DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
4050:     }
4051:     else if (maxChildId == -1) {
4052:       PetscSectionGetDof(localCoarse,p,&numValues);
4053:     }
4054:     /* we will pack the column indices with the field offsets */
4055:     if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
4056:       /* also send the centroid, and the gradient */
4057:       numValues += dim * (1 + numFVcomps);
4058:     }
4059:     PetscSectionSetDof(rootValuesSec,p,numValues);
4060:   }
4061:   PetscSectionSetUp(rootValuesSec);
4062:   {
4063:     PetscInt          numRootValues;
4064:     const PetscScalar *coarseArray;

4066:     PetscSectionGetStorageSize(rootValuesSec,&numRootValues);
4067:     PetscMalloc1(numRootValues,&rootValues);
4068:     VecGetArrayRead(vecCoarseLocal,&coarseArray);
4069:     for (p = pStartC; p < pEndC; p++) {
4070:       PetscInt    numValues;
4071:       PetscInt    pValOff;
4072:       PetscScalar *pVal;
4073:       PetscInt    maxChildId = maxChildIds[p - pStartC];

4075:       PetscSectionGetDof(rootValuesSec,p,&numValues);
4076:       if (!numValues) {
4077:         continue;
4078:       }
4079:       PetscSectionGetOffset(rootValuesSec,p,&pValOff);
4080:       pVal = &(rootValues[pValOff]);
4081:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
4082:         PetscInt closureSize = numValues;
4083:         DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);
4084:         if (grad && p >= cellStart && p < cellEnd) {
4085:           PetscFVCellGeom *cg;
4086:           PetscScalar     *gradVals = NULL;
4087:           PetscInt        i;

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

4091:           DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);
4092:           for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
4093:           pVal += dim;
4094:           DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);
4095:           for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
4096:         }
4097:       }
4098:       else if (maxChildId == -1) {
4099:         PetscInt lDof, lOff, i;

4101:         PetscSectionGetDof(localCoarse,p,&lDof);
4102:         PetscSectionGetOffset(localCoarse,p,&lOff);
4103:         for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
4104:       }
4105:     }
4106:     VecRestoreArrayRead(vecCoarseLocal,&coarseArray);
4107:     PetscFree(maxChildIds);
4108:   }
4109:   {
4110:     PetscSF  valuesSF;
4111:     PetscInt *remoteOffsetsValues, numLeafValues;

4113:     PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);
4114:     PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);
4115:     PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);
4116:     PetscSFDestroy(&coarseToFineEmbedded);
4117:     PetscFree(remoteOffsetsValues);
4118:     PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);
4119:     PetscMalloc1(numLeafValues,&leafValues);
4120:     PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues,MPI_REPLACE);
4121:     PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues,MPI_REPLACE);
4122:     PetscSFDestroy(&valuesSF);
4123:     PetscFree(rootValues);
4124:     PetscSectionDestroy(&rootValuesSec);
4125:   }
4126:   DMGetLocalSection(fine,&localFine);
4127:   {
4128:     PetscInt    maxDof;
4129:     PetscInt    *rowIndices;
4130:     DM           refTree;
4131:     PetscInt     **refPointFieldN;
4132:     PetscScalar  ***refPointFieldMats;
4133:     PetscSection refConSec, refAnSec;
4134:     PetscInt     pRefStart,pRefEnd,leafStart,leafEnd;
4135:     PetscScalar  *pointWork;

4137:     PetscSectionGetMaxDof(localFine,&maxDof);
4138:     DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
4139:     DMGetWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);
4140:     DMPlexGetReferenceTree(fine,&refTree);
4141:     DMCopyDisc(fine,refTree);
4142:     DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
4143:     DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL);
4144:     DMPlexGetAnchors(refTree,&refAnSec,NULL);
4145:     PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
4146:     PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);
4147:     DMPlexGetSimplexOrBoxCells(fine,0,&cellStart,&cellEnd);
4148:     for (p = leafStart; p < leafEnd; p++) {
4149:       PetscInt          gDof, gcDof, gOff, lDof;
4150:       PetscInt          numValues, pValOff;
4151:       PetscInt          childId;
4152:       const PetscScalar *pVal;
4153:       const PetscScalar *fvGradData = NULL;

4155:       PetscSectionGetDof(globalFine,p,&gDof);
4156:       PetscSectionGetDof(localFine,p,&lDof);
4157:       PetscSectionGetConstraintDof(globalFine,p,&gcDof);
4158:       if ((gDof - gcDof) <= 0) {
4159:         continue;
4160:       }
4161:       PetscSectionGetOffset(globalFine,p,&gOff);
4162:       PetscSectionGetDof(leafValuesSec,p,&numValues);
4163:       if (!numValues) continue;
4164:       PetscSectionGetOffset(leafValuesSec,p,&pValOff);
4165:       pVal = &leafValues[pValOff];
4166:       offsets[0]        = 0;
4167:       offsetsCopy[0]    = 0;
4168:       newOffsets[0]     = 0;
4169:       newOffsetsCopy[0] = 0;
4170:       childId           = cids[p - pStartF];
4171:       if (numFields) {
4172:         PetscInt f;
4173:         for (f = 0; f < numFields; f++) {
4174:           PetscInt rowDof;

4176:           PetscSectionGetFieldDof(localFine,p,f,&rowDof);
4177:           offsets[f + 1]        = offsets[f] + rowDof;
4178:           offsetsCopy[f + 1]    = offsets[f + 1];
4179:           /* TODO: closure indices */
4180:           newOffsets[f + 1]     = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
4181:         }
4182:         DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,NULL,rowIndices);
4183:       }
4184:       else {
4185:         offsets[0]    = 0;
4186:         offsets[1]    = lDof;
4187:         newOffsets[0] = 0;
4188:         newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
4189:         DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,NULL,rowIndices);
4190:       }
4191:       if (childId == -1) { /* no child interpolation: one nnz per */
4192:         VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);
4193:       } else {
4194:         PetscInt f;

4196:         if (grad && p >= cellStart && p < cellEnd) {
4197:           numValues -= (dim * (1 + numFVcomps));
4198:           fvGradData = &pVal[numValues];
4199:         }
4200:         for (f = 0; f < PetscMax(1,numFields); f++) {
4201:           const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
4202:           PetscInt numRows = offsets[f+1] - offsets[f];
4203:           PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
4204:           const PetscScalar *cVal = &pVal[newOffsets[f]];
4205:           PetscScalar *rVal = &pointWork[offsets[f]];
4206:           PetscInt i, j;

4208: #if 0
4209:           PetscInfo(coarse,"childId %D, numRows %D, numCols %D, refPointFieldN %D maxDof %D\n",childId,numRows,numCols,refPointFieldN[childId - pRefStart][f], maxDof);
4210: #endif
4211:           for (i = 0; i < numRows; i++) {
4212:             PetscScalar val = 0.;
4213:             for (j = 0; j < numCols; j++) {
4214:               val += childMat[i * numCols + j] * cVal[j];
4215:             }
4216:             rVal[i] = val;
4217:           }
4218:           if (f == fvField && p >= cellStart && p < cellEnd) {
4219:             PetscReal   centroid[3];
4220:             PetscScalar diff[3];
4221:             const PetscScalar *parentCentroid = &fvGradData[0];
4222:             const PetscScalar *gradient       = &fvGradData[dim];

4224:             DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);
4225:             for (i = 0; i < dim; i++) {
4226:               diff[i] = centroid[i] - parentCentroid[i];
4227:             }
4228:             for (i = 0; i < numFVcomps; i++) {
4229:               PetscScalar val = 0.;

4231:               for (j = 0; j < dim; j++) {
4232:                 val += gradient[dim * i + j] * diff[j];
4233:               }
4234:               rVal[i] += val;
4235:             }
4236:           }
4237:           VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);
4238:         }
4239:       }
4240:     }
4241:     DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
4242:     DMRestoreWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);
4243:     DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
4244:   }
4245:   PetscFree(leafValues);
4246:   PetscSectionDestroy(&leafValuesSec);
4247:   PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);
4248:   ISRestoreIndices(aIS,&anchors);
4249:   return 0;
4250: }

4252: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4253: {
4254:   DM             refTree;
4255:   PetscSection   multiRootSec, rootIndicesSec;
4256:   PetscSection   globalCoarse, globalFine;
4257:   PetscSection   localCoarse, localFine;
4258:   PetscSection   cSecRef;
4259:   PetscInt       *parentIndices, pRefStart, pRefEnd;
4260:   PetscScalar    *rootValues, *parentValues;
4261:   Mat            injRef;
4262:   PetscInt       numFields, maxDof;
4263:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
4264:   PetscInt       *offsets, *offsetsCopy, *rowOffsets;
4265:   PetscLayout    rowMap, colMap;
4266:   PetscInt       rowStart, rowEnd, colStart, colEnd;
4267:   PetscScalar    ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */


4270:   /* get the templates for the fine-to-coarse injection from the reference tree */
4271:   VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
4272:   VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
4273:   DMPlexGetReferenceTree(coarse,&refTree);
4274:   DMCopyDisc(coarse,refTree);
4275:   DMGetDefaultConstraints(refTree,&cSecRef,NULL,NULL);
4276:   PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
4277:   DMPlexReferenceTreeGetInjector(refTree,&injRef);

4279:   DMPlexGetChart(fine,&pStartF,&pEndF);
4280:   DMGetLocalSection(fine,&localFine);
4281:   DMGetGlobalSection(fine,&globalFine);
4282:   PetscSectionGetNumFields(localFine,&numFields);
4283:   DMPlexGetChart(coarse,&pStartC,&pEndC);
4284:   DMGetLocalSection(coarse,&localCoarse);
4285:   DMGetGlobalSection(coarse,&globalCoarse);
4286:   PetscSectionGetMaxDof(localCoarse,&maxDof);
4287:   {
4288:     PetscInt maxFields = PetscMax(1,numFields) + 1;
4289:     PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
4290:   }

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

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

4296:   /* count indices */
4297:   VecGetLayout(vecFine,&colMap);
4298:   VecGetLayout(vecCoarse,&rowMap);
4299:   PetscLayoutSetUp(rowMap);
4300:   PetscLayoutSetUp(colMap);
4301:   PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
4302:   PetscLayoutGetRange(colMap,&colStart,&colEnd);
4303:   /* insert values */
4304:   DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
4305:   for (p = pStartC; p < pEndC; p++) {
4306:     PetscInt  numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4307:     PetscBool contribute = PETSC_FALSE;

4309:     PetscSectionGetDof(globalCoarse,p,&dof);
4310:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
4311:     if ((dof - cdof) <= 0) continue;
4312:     PetscSectionGetDof(localCoarse,p,&dof);
4313:     PetscSectionGetOffset(globalCoarse,p,&gOff);

4315:     rowOffsets[0] = 0;
4316:     offsetsCopy[0] = 0;
4317:     if (numFields) {
4318:       PetscInt f;

4320:       for (f = 0; f < numFields; f++) {
4321:         PetscInt fDof;
4322:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
4323:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
4324:       }
4325:       DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,NULL,parentIndices);
4326:     } else {
4327:       DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,NULL,parentIndices);
4328:       rowOffsets[1] = offsetsCopy[0];
4329:     }

4331:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
4332:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
4333:     leafEnd = leafStart + numLeaves;
4334:     for (l = 0; l < dof; l++) parentValues[l] = 0.;
4335:     for (l = leafStart; l < leafEnd; l++) {
4336:       PetscInt numIndices, childId, offset;
4337:       const PetscScalar *childValues;

4339:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
4340:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
4341:       childId = (PetscInt) PetscRealPart(rootValues[offset++]);
4342:       childValues = &rootValues[offset];
4343:       numIndices--;

4345:       if (childId == -2) { /* skip */
4346:         continue;
4347:       } else if (childId == -1) { /* equivalent points: scatter */
4348:         PetscInt m;

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

4355:         contribute = PETSC_TRUE;
4356:         DMPlexGetTreeParent(refTree,childId,&parentId,NULL);

4358:         lim = PetscMax(1,numFields);
4359:         offsets[0] = 0;
4360:         if (numFields) {
4361:           PetscInt f;

4363:           for (f = 0; f < numFields; f++) {
4364:             PetscInt fDof;
4365:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

4367:             offsets[f + 1] = fDof + offsets[f];
4368:           }
4369:         }
4370:         else {
4371:           PetscInt cDof;

4373:           PetscSectionGetDof(cSecRef,childId,&cDof);
4374:           offsets[1] = cDof;
4375:         }
4376:         for (f = 0; f < lim; f++) {
4377:           PetscScalar       *childMat   = &childrenMats[childId - pRefStart][f][0];
4378:           PetscInt          n           = offsets[f+1]-offsets[f];
4379:           PetscInt          m           = rowOffsets[f+1]-rowOffsets[f];
4380:           PetscInt          i, j;
4381:           const PetscScalar *colValues  = &childValues[offsets[f]];

4383:           for (i = 0; i < m; i++) {
4384:             PetscScalar val = 0.;
4385:             for (j = 0; j < n; j++) {
4386:               val += childMat[n * i + j] * colValues[j];
4387:             }
4388:             parentValues[rowOffsets[f] + i] += val;
4389:           }
4390:         }
4391:       }
4392:     }
4393:     if (contribute) VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);
4394:   }
4395:   PetscSectionDestroy(&multiRootSec);
4396:   PetscSectionDestroy(&rootIndicesSec);
4397:   PetscFree2(parentIndices,parentValues);
4398:   DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
4399:   PetscFree(rootValues);
4400:   PetscFree3(offsets,offsetsCopy,rowOffsets);
4401:   return 0;
4402: }

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

4409:   collective

4411:   Input Parameters:
4412: + dmIn        - The DMPlex mesh for the input vector
4413: . vecIn       - The input vector
4414: . sfRefine    - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in
4415:                 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn
4416: . sfCoarsen   - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in
4417:                 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn
4418: . cidsRefine  - The childIds of the points in dmOut.  These childIds relate back to the reference tree: childid[j] = k implies
4419:                 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference
4420:                 tree was refined from its parent.  childid[j] = -1 indicates that the point j in dmOut is exactly
4421:                 equivalent to its root in dmIn, so no interpolation is necessary.  childid[j] = -2 indicates that this
4422:                 point j in dmOut is not a leaf of sfRefine.
4423: . cidsCoarsen - The childIds of the points in dmIn.  These childIds relate back to the reference tree: childid[j] = k implies
4424:                 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference
4425:                 tree coarsens to its parent.  childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen.
4426: . useBCs      - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer.
4427: - time        - Used if boundary values are time dependent.

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

4435:   Level: developer

4437: .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients()
4438: @*/
4439: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4440: {
4441:   VecSet(vecOut,0.0);
4442:   if (sfRefine) {
4443:     Vec vecInLocal;
4444:     DM  dmGrad = NULL;
4445:     Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;

4447:     DMGetLocalVector(dmIn,&vecInLocal);
4448:     VecSet(vecInLocal,0.0);
4449:     {
4450:       PetscInt  numFields, i;

4452:       DMGetNumFields(dmIn, &numFields);
4453:       for (i = 0; i < numFields; i++) {
4454:         PetscObject  obj;
4455:         PetscClassId classid;

4457:         DMGetField(dmIn, i, NULL, &obj);
4458:         PetscObjectGetClassId(obj, &classid);
4459:         if (classid == PETSCFV_CLASSID) {
4460:           DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);
4461:           break;
4462:         }
4463:       }
4464:     }
4465:     if (useBCs) {
4466:       DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);
4467:     }
4468:     DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);
4469:     DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);
4470:     if (dmGrad) {
4471:       DMGetGlobalVector(dmGrad,&grad);
4472:       DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);
4473:     }
4474:     DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);
4475:     DMRestoreLocalVector(dmIn,&vecInLocal);
4476:     if (dmGrad) {
4477:       DMRestoreGlobalVector(dmGrad,&grad);
4478:     }
4479:   }
4480:   if (sfCoarsen) {
4481:     DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);
4482:   }
4483:   VecAssemblyBegin(vecOut);
4484:   VecAssemblyEnd(vecOut);
4485:   return 0;
4486: }