Actual source code: plextree.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <../src/sys/utils/hash.h>
  3: #include <petsc/private/isimpl.h>
  4: #include <petsc/private/petscfeimpl.h>
  5: #include <petscsf.h>
  6: #include <petscds.h>

  8: /** hierarchy routines */

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

 15:   Not collective

 17:   Input Parameters:
 18: + dm - The DMPlex object
 19: - ref - The reference tree DMPlex object

 21:   Level: intermediate

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

 33:   PetscObjectReference((PetscObject)ref);
 34:   DMDestroy(&mesh->referenceTree);
 35:   mesh->referenceTree = ref;
 36:   return(0);
 37: }

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

 44:   Not collective

 46:   Input Parameters:
 47: . dm - The DMPlex object

 49:   Output Parameters
 50: . ref - The reference tree DMPlex object

 52:   Level: intermediate

 54: .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree()
 55: @*/
 56: PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref)
 57: {
 58:   DM_Plex        *mesh = (DM_Plex *)dm->data;

 63:   *ref = mesh->referenceTree;
 64:   return(0);
 65: }

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

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

 93:     DMPlexGetSupportSize(dm,childA,&size);
 94:     DMPlexGetSupport(dm,childA,&supp);

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

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

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

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

177:       posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize);
178:       if (childB) *childB = children[posB];
179:     }
180:   }
181:   if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap);
182:   return(0);
183: }

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

190:   Input Parameters:
191: + dm - the reference tree DMPlex object
192: . parent - the parent point
193: . parentOrientA - the reference orientation for describing the parent
194: . childOrientA - the reference orientation for describing the child
195: . childA - the reference childID for describing the child
196: - parentOrientB - the new orientation for describing the parent

198:   Output Parameters:
199: + childOrientB - if not NULL, set to the new oreintation for describing the child
200: . childB - if not NULL, the new childID for describing the child

202:   Level: developer

204: .seealso: DMPlexGetReferenceTree(), DMPlexSetReferenceTree(), DMPlexSetTree()
205: @*/
206: PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
207: {
208:   DM_Plex        *mesh = (DM_Plex *)dm->data;

213:   if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented");
214:   mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);
215:   return(0);
216: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

450:   Collective on comm

452:   Input Parameters:
453: + comm    - the MPI communicator
454: . dim     - the spatial dimension
455: - simplex - Flag for simplex, otherwise use a tensor-product cell

457:   Output Parameters:
458: . ref     - the reference tree DMPlex object

460:   Level: intermediate

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

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

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

500: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
501: {
502:   DM_Plex        *mesh = (DM_Plex *)dm->data;
503:   PetscSection   childSec, pSec;
504:   PetscInt       p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT;
505:   PetscInt       *offsets, *children, pStart, pEnd;

510:   PetscSectionDestroy(&mesh->childSection);
511:   PetscFree(mesh->children);
512:   pSec = mesh->parentSection;
513:   if (!pSec) return(0);
514:   PetscSectionGetStorageSize(pSec,&pSize);
515:   for (p = 0; p < pSize; p++) {
516:     PetscInt par = mesh->parents[p];

518:     parMax = PetscMax(parMax,par+1);
519:     parMin = PetscMin(parMin,par);
520:   }
521:   if (parMin > parMax) {
522:     parMin = -1;
523:     parMax = -1;
524:   }
525:   PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);
526:   PetscSectionSetChart(childSec,parMin,parMax);
527:   for (p = 0; p < pSize; p++) {
528:     PetscInt par = mesh->parents[p];

530:     PetscSectionAddDof(childSec,par,1);
531:   }
532:   PetscSectionSetUp(childSec);
533:   PetscSectionGetStorageSize(childSec,&cSize);
534:   PetscMalloc1(cSize,&children);
535:   PetscCalloc1(parMax-parMin,&offsets);
536:   PetscSectionGetChart(pSec,&pStart,&pEnd);
537:   for (p = pStart; p < pEnd; p++) {
538:     PetscInt dof, off, i;

540:     PetscSectionGetDof(pSec,p,&dof);
541:     PetscSectionGetOffset(pSec,p,&off);
542:     for (i = 0; i < dof; i++) {
543:       PetscInt par = mesh->parents[off + i], cOff;

545:       PetscSectionGetOffset(childSec,par,&cOff);
546:       children[cOff + offsets[par-parMin]++] = p;
547:     }
548:   }
549:   mesh->childSection = childSec;
550:   mesh->children = children;
551:   PetscFree(offsets);
552:   return(0);
553: }

557: static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
558: {
559:   PetscInt       pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
560:   const PetscInt *vals;
561:   PetscSection   secNew;
562:   PetscBool      anyNew, globalAnyNew;
563:   PetscBool      compress;

567:   PetscSectionGetChart(section,&pStart,&pEnd);
568:   ISGetLocalSize(is,&size);
569:   ISGetIndices(is,&vals);
570:   PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);
571:   PetscSectionSetChart(secNew,pStart,pEnd);
572:   for (i = 0; i < size; i++) {
573:     PetscInt dof;

575:     p = vals[i];
576:     if (p < pStart || p >= pEnd) continue;
577:     PetscSectionGetDof(section, p, &dof);
578:     if (dof) break;
579:   }
580:   if (i == size) {
581:     PetscSectionSetUp(secNew);
582:     anyNew   = PETSC_FALSE;
583:     compress = PETSC_FALSE;
584:     sizeNew  = 0;
585:   }
586:   else {
587:     anyNew = PETSC_TRUE;
588:     for (p = pStart; p < pEnd; p++) {
589:       PetscInt dof, off;

591:       PetscSectionGetDof(section, p, &dof);
592:       PetscSectionGetOffset(section, p, &off);
593:       for (i = 0; i < dof; i++) {
594:         PetscInt q = vals[off + i], qDof = 0;

596:         if (q >= pStart && q < pEnd) {
597:           PetscSectionGetDof(section, q, &qDof);
598:         }
599:         if (qDof) {
600:           PetscSectionAddDof(secNew, p, qDof);
601:         }
602:         else {
603:           PetscSectionAddDof(secNew, p, 1);
604:         }
605:       }
606:     }
607:     PetscSectionSetUp(secNew);
608:     PetscSectionGetStorageSize(secNew,&sizeNew);
609:     PetscMalloc1(sizeNew,&valsNew);
610:     compress = PETSC_FALSE;
611:     for (p = pStart; p < pEnd; p++) {
612:       PetscInt dof, off, count, offNew, dofNew;

614:       PetscSectionGetDof(section, p, &dof);
615:       PetscSectionGetOffset(section, p, &off);
616:       PetscSectionGetDof(secNew, p, &dofNew);
617:       PetscSectionGetOffset(secNew, p, &offNew);
618:       count = 0;
619:       for (i = 0; i < dof; i++) {
620:         PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;

622:         if (q >= pStart && q < pEnd) {
623:           PetscSectionGetDof(section, q, &qDof);
624:           PetscSectionGetOffset(section, q, &qOff);
625:         }
626:         if (qDof) {
627:           PetscInt oldCount = count;

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

632:             for (k = 0; k < oldCount; k++) {
633:               if (valsNew[offNew + k] == r) {
634:                 break;
635:               }
636:             }
637:             if (k == oldCount) {
638:               valsNew[offNew + count++] = r;
639:             }
640:           }
641:         }
642:         else {
643:           PetscInt k, oldCount = count;

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

671:     MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
672:     if (compress) {
673:       PetscSection secComp;
674:       PetscInt *valsComp = NULL;

676:       PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);
677:       PetscSectionSetChart(secComp,pStart,pEnd);
678:       for (p = pStart; p < pEnd; p++) {
679:         PetscInt dof;

681:         PetscSectionGetDof(secNew, p, &dof);
682:         PetscSectionSetDof(secComp, p, dof);
683:       }
684:       PetscSectionSetUp(secComp);
685:       PetscSectionGetStorageSize(secComp,&sizeNew);
686:       PetscMalloc1(sizeNew,&valsComp);
687:       for (p = pStart; p < pEnd; p++) {
688:         PetscInt dof, off, offNew, j;

690:         PetscSectionGetDof(secNew, p, &dof);
691:         PetscSectionGetOffset(secNew, p, &off);
692:         PetscSectionGetOffset(secComp, p, &offNew);
693:         for (j = 0; j < dof; j++) {
694:           valsComp[offNew + j] = valsNew[off + j];
695:         }
696:       }
697:       PetscSectionDestroy(&secNew);
698:       secNew  = secComp;
699:       PetscFree(valsNew);
700:       valsNew = valsComp;
701:     }
702:     ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);
703:   }
704:   return(0);
705: }

709: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
710: {
711:   PetscInt       p, pStart, pEnd, *anchors, size;
712:   PetscInt       aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
713:   PetscSection   aSec;
714:   DMLabel        canonLabel;
715:   IS             aIS;

720:   DMPlexGetChart(dm,&pStart,&pEnd);
721:   DMGetLabel(dm,"canonical",&canonLabel);
722:   for (p = pStart; p < pEnd; p++) {
723:     PetscInt parent;

725:     if (canonLabel) {
726:       PetscInt canon;

728:       DMLabelGetValue(canonLabel,p,&canon);
729:       if (p != canon) continue;
730:     }
731:     DMPlexGetTreeParent(dm,p,&parent,NULL);
732:     if (parent != p) {
733:       aMin = PetscMin(aMin,p);
734:       aMax = PetscMax(aMax,p+1);
735:     }
736:   }
737:   if (aMin > aMax) {
738:     aMin = -1;
739:     aMax = -1;
740:   }
741:   PetscSectionCreate(PETSC_COMM_SELF,&aSec);
742:   PetscSectionSetChart(aSec,aMin,aMax);
743:   for (p = aMin; p < aMax; p++) {
744:     PetscInt parent, ancestor = p;

746:     if (canonLabel) {
747:       PetscInt canon;

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

760:       DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
761:       PetscSectionSetDof(aSec,p,closureSize);
762:       DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
763:     }
764:   }
765:   PetscSectionSetUp(aSec);
766:   PetscSectionGetStorageSize(aSec,&size);
767:   PetscMalloc1(size,&anchors);
768:   for (p = aMin; p < aMax; p++) {
769:     PetscInt parent, ancestor = p;

771:     if (canonLabel) {
772:       PetscInt canon;

774:       DMLabelGetValue(canonLabel,p,&canon);
775:       if (p != canon) continue;
776:     }
777:     DMPlexGetTreeParent(dm,p,&parent,NULL);
778:     while (parent != ancestor) {
779:       ancestor = parent;
780:       DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
781:     }
782:     if (ancestor != p) {
783:       PetscInt j, closureSize, *closure = NULL, aOff;

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

787:       DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
788:       for (j = 0; j < closureSize; j++) {
789:         anchors[aOff + j] = closure[2*j];
790:       }
791:       DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
792:     }
793:   }
794:   ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);
795:   {
796:     PetscSection aSecNew = aSec;
797:     IS           aISNew  = aIS;

799:     PetscObjectReference((PetscObject)aSec);
800:     PetscObjectReference((PetscObject)aIS);
801:     while (aSecNew) {
802:       PetscSectionDestroy(&aSec);
803:       ISDestroy(&aIS);
804:       aSec    = aSecNew;
805:       aIS     = aISNew;
806:       aSecNew = NULL;
807:       aISNew  = NULL;
808:       AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);
809:     }
810:   }
811:   DMPlexSetAnchors(dm,aSec,aIS);
812:   PetscSectionDestroy(&aSec);
813:   ISDestroy(&aIS);
814:   return(0);
815: }

819: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp)
820: {

824:   if (numTrueSupp[p] == -1) {
825:     PetscInt i, alldof;
826:     const PetscInt *supp;
827:     PetscInt count = 0;

829:     DMPlexGetSupportSize(dm,p,&alldof);
830:     DMPlexGetSupport(dm,p,&supp);
831:     for (i = 0; i < alldof; i++) {
832:       PetscInt q = supp[i], numCones, j;
833:       const PetscInt *cone;

835:       DMPlexGetConeSize(dm,q,&numCones);
836:       DMPlexGetCone(dm,q,&cone);
837:       for (j = 0; j < numCones; j++) {
838:         if (cone[j] == p) break;
839:       }
840:       if (j < numCones) count++;
841:     }
842:     numTrueSupp[p] = count;
843:   }
844:   *dof = numTrueSupp[p];
845:   return(0);
846: }

850: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
851: {
852:   DM_Plex *mesh = (DM_Plex *)dm->data;
853:   PetscSection newSupportSection;
854:   PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth;
855:   PetscInt *numTrueSupp;
856:   PetscInt *offsets;

861:   /* symmetrize the hierarchy */
862:   DMPlexGetDepth(dm,&depth);
863:   PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);
864:   DMPlexGetChart(dm,&pStart,&pEnd);
865:   PetscSectionSetChart(newSupportSection,pStart,pEnd);
866:   PetscCalloc1(pEnd,&offsets);
867:   PetscMalloc1(pEnd,&numTrueSupp);
868:   for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
869:   /* if a point is in the (true) support of q, it should be in the support of
870:    * parent(q) */
871:   for (d = 0; d <= depth; d++) {
872:     DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
873:     for (p = pStart; p < pEnd; ++p) {
874:       PetscInt dof, q, qdof, parent;

876:       DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);
877:       PetscSectionAddDof(newSupportSection, p, dof);
878:       q    = p;
879:       DMPlexGetTreeParent(dm,q,&parent,NULL);
880:       while (parent != q && parent >= pStart && parent < pEnd) {
881:         q = parent;

883:         DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);
884:         PetscSectionAddDof(newSupportSection,p,qdof);
885:         PetscSectionAddDof(newSupportSection,q,dof);
886:         DMPlexGetTreeParent(dm,q,&parent,NULL);
887:       }
888:     }
889:   }
890:   PetscSectionSetUp(newSupportSection);
891:   PetscSectionGetStorageSize(newSupportSection,&newSize);
892:   PetscMalloc1(newSize,&newSupports);
893:   for (d = 0; d <= depth; d++) {
894:     DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
895:     for (p = pStart; p < pEnd; p++) {
896:       PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;

898:       PetscSectionGetDof(mesh->supportSection, p, &dof);
899:       PetscSectionGetOffset(mesh->supportSection, p, &off);
900:       PetscSectionGetDof(newSupportSection, p, &newDof);
901:       PetscSectionGetOffset(newSupportSection, p, &newOff);
902:       for (i = 0; i < dof; i++) {
903:         PetscInt numCones, j;
904:         const PetscInt *cone;
905:         PetscInt q = mesh->supports[off + i];

907:         DMPlexGetConeSize(dm,q,&numCones);
908:         DMPlexGetCone(dm,q,&cone);
909:         for (j = 0; j < numCones; j++) {
910:           if (cone[j] == p) break;
911:         }
912:         if (j < numCones) newSupports[newOff+offsets[p]++] = q;
913:       }
914:       mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof);

916:       q    = p;
917:       DMPlexGetTreeParent(dm,q,&parent,NULL);
918:       while (parent != q && parent >= pStart && parent < pEnd) {
919:         q = parent;
920:         PetscSectionGetDof(mesh->supportSection, q, &qdof);
921:         PetscSectionGetOffset(mesh->supportSection, q, &qoff);
922:         PetscSectionGetOffset(newSupportSection, q, &newqOff);
923:         for (i = 0; i < qdof; i++) {
924:           PetscInt numCones, j;
925:           const PetscInt *cone;
926:           PetscInt r = mesh->supports[qoff + i];

928:           DMPlexGetConeSize(dm,r,&numCones);
929:           DMPlexGetCone(dm,r,&cone);
930:           for (j = 0; j < numCones; j++) {
931:             if (cone[j] == q) break;
932:           }
933:           if (j < numCones) newSupports[newOff+offsets[p]++] = r;
934:         }
935:         for (i = 0; i < dof; i++) {
936:           PetscInt numCones, j;
937:           const PetscInt *cone;
938:           PetscInt r = mesh->supports[off + i];

940:           DMPlexGetConeSize(dm,r,&numCones);
941:           DMPlexGetCone(dm,r,&cone);
942:           for (j = 0; j < numCones; j++) {
943:             if (cone[j] == p) break;
944:           }
945:           if (j < numCones) newSupports[newqOff+offsets[q]++] = r;
946:         }
947:         DMPlexGetTreeParent(dm,q,&parent,NULL);
948:       }
949:     }
950:   }
951:   PetscSectionDestroy(&mesh->supportSection);
952:   mesh->supportSection = newSupportSection;
953:   PetscFree(mesh->supports);
954:   mesh->supports = newSupports;
955:   PetscFree(offsets);
956:   PetscFree(numTrueSupp);

958:   return(0);
959: }

961: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat);
962: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat);

966: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
967: {
968:   DM_Plex       *mesh = (DM_Plex *)dm->data;
969:   DM             refTree;
970:   PetscInt       size;

976:   PetscObjectReference((PetscObject)parentSection);
977:   PetscSectionDestroy(&mesh->parentSection);
978:   mesh->parentSection = parentSection;
979:   PetscSectionGetStorageSize(parentSection,&size);
980:   if (parents != mesh->parents) {
981:     PetscFree(mesh->parents);
982:     PetscMalloc1(size,&mesh->parents);
983:     PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));
984:   }
985:   if (childIDs != mesh->childIDs) {
986:     PetscFree(mesh->childIDs);
987:     PetscMalloc1(size,&mesh->childIDs);
988:     PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));
989:   }
990:   DMPlexGetReferenceTree(dm,&refTree);
991:   if (refTree) {
992:     DMLabel canonLabel;

994:     DMGetLabel(refTree,"canonical",&canonLabel);
995:     if (canonLabel) {
996:       PetscInt i;

998:       for (i = 0; i < size; i++) {
999:         PetscInt canon;
1000:         DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);
1001:         if (canon >= 0) {
1002:           mesh->childIDs[i] = canon;
1003:         }
1004:       }
1005:     }
1006:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
1007:   }
1008:   else {
1009:     mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
1010:   }
1011:   DMPlexTreeSymmetrize(dm);
1012:   if (computeCanonical) {
1013:     PetscInt d, dim;

1015:     /* add the canonical label */
1016:     DMGetDimension(dm,&dim);
1017:     DMCreateLabel(dm,"canonical");
1018:     for (d = 0; d <= dim; d++) {
1019:       PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
1020:       const PetscInt *cChildren;

1022:       DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
1023:       for (p = dStart; p < dEnd; p++) {
1024:         DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);
1025:         if (cNumChildren) {
1026:           canon = p;
1027:           break;
1028:         }
1029:       }
1030:       if (canon == -1) continue;
1031:       for (p = dStart; p < dEnd; p++) {
1032:         PetscInt numChildren, i;
1033:         const PetscInt *children;

1035:         DMPlexGetTreeChildren(dm,p,&numChildren,&children);
1036:         if (numChildren) {
1037:           if (numChildren != cNumChildren) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"All parent points in a stratum should have the same number of children: %d != %d", numChildren, cNumChildren);
1038:           DMSetLabelValue(dm,"canonical",p,canon);
1039:           for (i = 0; i < numChildren; i++) {
1040:             DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);
1041:           }
1042:         }
1043:       }
1044:     }
1045:   }
1046:   if (exchangeSupports) {
1047:     DMPlexTreeExchangeSupports(dm);
1048:   }
1049:   mesh->createanchors = DMPlexCreateAnchors_Tree;
1050:   /* reset anchors */
1051:   DMPlexSetAnchors(dm,NULL,NULL);
1052:   return(0);
1053: }

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

1062:   Collective on dm

1064:   Input Parameters:
1065: + dm - the DMPlex object
1066: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1067:                   offset indexes the parent and childID list; the reference count of parentSection is incremented
1068: . parents - a list of the point parents; copied, can be destroyed
1069: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1070:              the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed

1072:   Level: intermediate

1074: .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1075: @*/
1076: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
1077: {

1081:   DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);
1082:   return(0);
1083: }

1087: /*@
1088:   DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1089:   Collective on dm

1091:   Input Parameters:
1092: . dm - the DMPlex object

1094:   Output Parameters:
1095: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1096:                   offset indexes the parent and childID list
1097: . parents - a list of the point parents
1098: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1099:              the child corresponds to the point in the reference tree with index childID
1100: . childSection - the inverse of the parent section
1101: - children - a list of the point children

1103:   Level: intermediate

1105: .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1106: @*/
1107: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1108: {
1109:   DM_Plex        *mesh = (DM_Plex *)dm->data;

1113:   if (parentSection) *parentSection = mesh->parentSection;
1114:   if (parents)       *parents       = mesh->parents;
1115:   if (childIDs)      *childIDs      = mesh->childIDs;
1116:   if (childSection)  *childSection  = mesh->childSection;
1117:   if (children)      *children      = mesh->children;
1118:   return(0);
1119: }

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

1126:   Input Parameters:
1127: + dm - the DMPlex object
1128: - point - the query point

1130:   Output Parameters:
1131: + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent
1132: - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point
1133:             does not have a parent

1135:   Level: intermediate

1137: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren()
1138: @*/
1139: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1140: {
1141:   DM_Plex       *mesh = (DM_Plex *)dm->data;
1142:   PetscSection   pSec;

1147:   pSec = mesh->parentSection;
1148:   if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1149:     PetscInt dof;

1151:     PetscSectionGetDof (pSec, point, &dof);
1152:     if (dof) {
1153:       PetscInt off;

1155:       PetscSectionGetOffset (pSec, point, &off);
1156:       if (parent)  *parent = mesh->parents[off];
1157:       if (childID) *childID = mesh->childIDs[off];
1158:       return(0);
1159:     }
1160:   }
1161:   if (parent) {
1162:     *parent = point;
1163:   }
1164:   if (childID) {
1165:     *childID = 0;
1166:   }
1167:   return(0);
1168: }

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

1175:   Input Parameters:
1176: + dm - the DMPlex object
1177: - point - the query point

1179:   Output Parameters:
1180: + numChildren - if not NULL, set to the number of children
1181: - children - if not NULL, set to a list children, or set to NULL if the point has no children

1183:   Level: intermediate

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

1189: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent()
1190: @*/
1191: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1192: {
1193:   DM_Plex       *mesh = (DM_Plex *)dm->data;
1194:   PetscSection   childSec;
1195:   PetscInt       dof = 0;

1200:   childSec = mesh->childSection;
1201:   if (childSec && point >= childSec->pStart && point < childSec->pEnd) {
1202:     PetscSectionGetDof (childSec, point, &dof);
1203:   }
1204:   if (numChildren) *numChildren = dof;
1205:   if (children) {
1206:     if (dof) {
1207:       PetscInt off;

1209:       PetscSectionGetOffset (childSec, point, &off);
1210:       *children = &mesh->children[off];
1211:     }
1212:     else {
1213:       *children = NULL;
1214:     }
1215:   }
1216:   return(0);
1217: }

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

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

1243:   for (f = 0; f < numFields; f++) {
1244:     PetscObject disc;
1245:     PetscFE fe = NULL;
1246:     PetscFV fv = NULL;
1247:     PetscClassId id;
1248:     PetscDualSpace space;
1249:     PetscInt i, j, k, nPoints, offset;
1250:     PetscInt fSize, fComp;
1251:     PetscReal *B = NULL;
1252:     PetscReal *weights, *pointsRef, *pointsReal;
1253:     Mat Amat, Bmat, Xmat;

1255:     PetscDSGetDiscretization(ds,f,&disc);
1256:     PetscObjectGetClassId(disc,&id);
1257:     if (id == PETSCFE_CLASSID) {
1258:       fe = (PetscFE) disc;
1259:       PetscFEGetDualSpace(fe,&space);
1260:       PetscDualSpaceGetDimension(space,&fSize);
1261:       PetscFEGetNumComponents(fe,&fComp);
1262:     }
1263:     else if (id == PETSCFV_CLASSID) {
1264:       fv = (PetscFV) disc;
1265:       PetscFVGetDualSpace(fv,&space);
1266:       PetscDualSpaceGetDimension(space,&fSize);
1267:       PetscFVGetNumComponents(fv,&fComp);
1268:     }
1269:     else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);

1271:     MatCreate(PETSC_COMM_SELF,&Amat);
1272:     MatSetSizes(Amat,fSize,fSize,fSize,fSize);
1273:     MatSetType(Amat,MATSEQDENSE);
1274:     MatSetUp(Amat);
1275:     MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);
1276:     MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);
1277:     nPoints = 0;
1278:     for (i = 0; i < fSize; i++) {
1279:       PetscInt        qPoints;
1280:       PetscQuadrature quad;

1282:       PetscDualSpaceGetFunctional(space,i,&quad);
1283:       PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);
1284:       nPoints += qPoints;
1285:     }
1286:     PetscMalloc3(nPoints,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal);
1287:     offset = 0;
1288:     for (i = 0; i < fSize; i++) {
1289:       PetscInt        qPoints;
1290:       const PetscReal    *p, *w;
1291:       PetscQuadrature quad;

1293:       PetscDualSpaceGetFunctional(space,i,&quad);
1294:       PetscQuadratureGetData(quad,NULL,&qPoints,&p,&w);
1295:       PetscMemcpy(weights+offset,w,qPoints*sizeof(*w));
1296:       PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));
1297:       offset += qPoints;
1298:     }
1299:     if (id == PETSCFE_CLASSID) {
1300:       PetscFEGetTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);
1301:     }
1302:     else {
1303:       PetscFVGetTabulation(fv,nPoints,pointsRef,&B,NULL,NULL);
1304:     }
1305:     offset = 0;
1306:     for (i = 0; i < fSize; i++) {
1307:       PetscInt        qPoints;
1308:       PetscQuadrature quad;

1310:       PetscDualSpaceGetFunctional(space,i,&quad);
1311:       PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);
1312:       for (j = 0; j < fSize; j++) {
1313:         PetscScalar val = 0.;

1315:         for (k = 0; k < qPoints; k++) {
1316:           val += B[((offset + k) * fSize + j) * fComp] * weights[k];
1317:         }
1318:         MatSetValue(Amat,i,j,val,INSERT_VALUES);
1319:       }
1320:       offset += qPoints;
1321:     }
1322:     if (id == PETSCFE_CLASSID) {
1323:       PetscFERestoreTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);
1324:     }
1325:     else {
1326:       PetscFVRestoreTabulation(fv,nPoints,pointsRef,&B,NULL,NULL);
1327:     }
1328:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
1329:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
1330:     MatLUFactor(Amat,NULL,NULL,NULL);
1331:     for (c = cStart; c < cEnd; c++) {
1332:       PetscInt parent;
1333:       PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1334:       PetscInt *childOffsets, *parentOffsets;

1336:       DMPlexGetTreeParent(dm,c,&parent,NULL);
1337:       if (parent == c) continue;
1338:       DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1339:       for (i = 0; i < closureSize; i++) {
1340:         PetscInt p = closure[2*i];
1341:         PetscInt conDof;

1343:         if (p < conStart || p >= conEnd) continue;
1344:         if (numFields > 1) {
1345:           PetscSectionGetFieldDof(cSec,p,f,&conDof);
1346:         }
1347:         else {
1348:           PetscSectionGetDof(cSec,p,&conDof);
1349:         }
1350:         if (conDof) break;
1351:       }
1352:       if (i == closureSize) {
1353:         DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1354:         continue;
1355:       }

1357:       DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ);
1358:       DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent);
1359:       for (i = 0; i < nPoints; i++) {
1360:         CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp);
1361:         CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);
1362:       }
1363:       if (id == PETSCFE_CLASSID) {
1364:         PetscFEGetTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);
1365:       }
1366:       else {
1367:         PetscFVGetTabulation(fv,nPoints,pointsReal,&B,NULL,NULL);
1368:       }
1369:       offset = 0;
1370:       for (i = 0; i < fSize; i++) {
1371:         PetscInt        qPoints;
1372:         PetscQuadrature quad;

1374:         PetscDualSpaceGetFunctional(space,i,&quad);
1375:         PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);
1376:         for (j = 0; j < fSize; j++) {
1377:           PetscScalar val = 0.;

1379:           for (k = 0; k < qPoints; k++) {
1380:             val += B[((offset + k) * fSize + j) * fComp] * weights[k];
1381:           }
1382:           MatSetValue(Bmat,i,j,val,INSERT_VALUES);
1383:         }
1384:         offset += qPoints;
1385:       }
1386:       if (id == PETSCFE_CLASSID) {
1387:         PetscFERestoreTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);
1388:       }
1389:       else {
1390:         PetscFVRestoreTabulation(fv,nPoints,pointsReal,&B,NULL,NULL);
1391:       }
1392:       MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY);
1393:       MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY);
1394:       MatMatSolve(Amat,Bmat,Xmat);
1395:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1396:       PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);
1397:       childOffsets[0] = 0;
1398:       for (i = 0; i < closureSize; i++) {
1399:         PetscInt p = closure[2*i];
1400:         PetscInt dof;

1402:         if (numFields > 1) {
1403:           PetscSectionGetFieldDof(section,p,f,&dof);
1404:         }
1405:         else {
1406:           PetscSectionGetDof(section,p,&dof);
1407:         }
1408:         childOffsets[i+1]=childOffsets[i]+dof / fComp;
1409:       }
1410:       parentOffsets[0] = 0;
1411:       for (i = 0; i < closureSizeP; i++) {
1412:         PetscInt p = closureP[2*i];
1413:         PetscInt dof;

1415:         if (numFields > 1) {
1416:           PetscSectionGetFieldDof(section,p,f,&dof);
1417:         }
1418:         else {
1419:           PetscSectionGetDof(section,p,&dof);
1420:         }
1421:         parentOffsets[i+1]=parentOffsets[i]+dof / fComp;
1422:       }
1423:       for (i = 0; i < closureSize; i++) {
1424:         PetscInt conDof, conOff, aDof, aOff;
1425:         PetscInt p = closure[2*i];
1426:         PetscInt o = closure[2*i+1];

1428:         if (p < conStart || p >= conEnd) continue;
1429:         if (numFields > 1) {
1430:           PetscSectionGetFieldDof(cSec,p,f,&conDof);
1431:           PetscSectionGetFieldOffset(cSec,p,f,&conOff);
1432:         }
1433:         else {
1434:           PetscSectionGetDof(cSec,p,&conDof);
1435:           PetscSectionGetOffset(cSec,p,&conOff);
1436:         }
1437:         if (!conDof) continue;
1438:         PetscSectionGetDof(aSec,p,&aDof);
1439:         PetscSectionGetOffset(aSec,p,&aOff);
1440:         for (k = 0; k < aDof; k++) {
1441:           PetscInt a = anchors[aOff + k];
1442:           PetscInt aSecDof, aSecOff;

1444:           if (numFields > 1) {
1445:             PetscSectionGetFieldDof(section,a,f,&aSecDof);
1446:             PetscSectionGetFieldOffset(section,a,f,&aSecOff);
1447:           }
1448:           else {
1449:             PetscSectionGetDof(section,a,&aSecDof);
1450:             PetscSectionGetOffset(section,a,&aSecOff);
1451:           }
1452:           if (!aSecDof) continue;

1454:           for (j = 0; j < closureSizeP; j++) {
1455:             PetscInt q = closureP[2*j];
1456:             PetscInt oq = closureP[2*j+1];

1458:             if (q == a) {
1459:               PetscInt r, s, t;

1461:               for (r = childOffsets[i]; r < childOffsets[i+1]; r++) {
1462:                 for (s = parentOffsets[j]; s < parentOffsets[j+1]; s++) {
1463:                   PetscScalar val;
1464:                   PetscInt insertCol, insertRow;

1466:                   MatGetValue(Xmat,r,s,&val);
1467:                   if (o >= 0) {
1468:                     insertRow = conOff + fComp * (r - childOffsets[i]);
1469:                   }
1470:                   else {
1471:                     insertRow = conOff + fComp * (childOffsets[i + 1] - 1 - r);
1472:                   }
1473:                   if (oq >= 0) {
1474:                     insertCol = aSecOff + fComp * (s - parentOffsets[j]);
1475:                   }
1476:                   else {
1477:                     insertCol = aSecOff + fComp * (parentOffsets[j + 1] - 1 - s);
1478:                   }
1479:                   for (t = 0; t < fComp; t++) {
1480:                     MatSetValue(cMat,insertRow + t,insertCol + t,val,INSERT_VALUES);
1481:                   }
1482:                 }
1483:               }
1484:             }
1485:           }
1486:         }
1487:       }
1488:       PetscFree2(childOffsets,parentOffsets);
1489:       DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1490:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1491:     }
1492:     MatDestroy(&Amat);
1493:     MatDestroy(&Bmat);
1494:     MatDestroy(&Xmat);
1495:     PetscFree3(weights,pointsRef,pointsReal);
1496:   }
1497:   MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1498:   MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1499:   PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);
1500:   ISRestoreIndices(aIS,&anchors);

1502:   return(0);
1503: }

1507: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1508: {
1509:   Mat            refCmat;
1510:   PetscDS        ds;
1511:   PetscInt       numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1512:   PetscScalar    ***refPointFieldMats;
1513:   PetscSection   refConSec, refAnSec, refSection;
1514:   IS             refAnIS;
1515:   const PetscInt *refAnchors;

1519:   DMGetDS(refTree,&ds);
1520:   PetscDSGetNumFields(ds,&numFields);
1521:   DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1522:   DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1523:   ISGetIndices(refAnIS,&refAnchors);
1524:   DMGetDefaultSection(refTree,&refSection);
1525:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1526:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
1527:   PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);
1528:   PetscSectionGetMaxDof(refConSec,&maxDof);
1529:   PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1530:   PetscMalloc1(maxDof,&rows);
1531:   PetscMalloc1(maxDof*maxAnDof,&cols);
1532:   for (p = pRefStart; p < pRefEnd; p++) {
1533:     PetscInt parent, closureSize, *closure = NULL, pDof;

1535:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
1536:     PetscSectionGetDof(refConSec,p,&pDof);
1537:     if (!pDof || parent == p) continue;

1539:     PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);
1540:     PetscCalloc1(numFields,&refPointFieldN[p-pRefStart]);
1541:     DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1542:     for (f = 0; f < numFields; f++) {
1543:       PetscInt cDof, cOff, numCols, r, i, fComp;
1544:       PetscObject disc;
1545:       PetscClassId id;
1546:       PetscFE fe = NULL;
1547:       PetscFV fv = NULL;

1549:       PetscDSGetDiscretization(ds,f,&disc);
1550:       PetscObjectGetClassId(disc,&id);
1551:       if (id == PETSCFE_CLASSID) {
1552:         fe = (PetscFE) disc;
1553:         PetscFEGetNumComponents(fe,&fComp);
1554:       }
1555:       else if (id == PETSCFV_CLASSID) {
1556:         fv = (PetscFV) disc;
1557:         PetscFVGetNumComponents(fv,&fComp);
1558:       }
1559:       else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);

1561:       if (numFields > 1) {
1562:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1563:         PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
1564:       }
1565:       else {
1566:         PetscSectionGetDof(refConSec,p,&cDof);
1567:         PetscSectionGetOffset(refConSec,p,&cOff);
1568:       }

1570:       if (!cDof) continue;
1571:       for (r = 0; r < cDof; r++) {
1572:         rows[r] = cOff + r;
1573:       }
1574:       numCols = 0;
1575:       for (i = 0; i < closureSize; i++) {
1576:         PetscInt q = closure[2*i];
1577:         PetscInt o = closure[2*i+1];
1578:         PetscInt aDof, aOff, j;

1580:         if (numFields > 1) {
1581:           PetscSectionGetFieldDof(refSection,q,f,&aDof);
1582:           PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1583:         }
1584:         else {
1585:           PetscSectionGetDof(refSection,q,&aDof);
1586:           PetscSectionGetOffset(refSection,q,&aOff);
1587:         }

1589:         for (j = 0; j < aDof; j++) {
1590:           PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp);
1591:           PetscInt comp = (j % fComp);

1593:           cols[numCols++] = aOff + node * fComp + comp;
1594:         }
1595:       }
1596:       refPointFieldN[p-pRefStart][f] = numCols;
1597:       PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
1598:       MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);
1599:     }
1600:     DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1601:   }
1602:   *childrenMats = refPointFieldMats;
1603:   *childrenN = refPointFieldN;
1604:   ISRestoreIndices(refAnIS,&refAnchors);
1605:   PetscFree(rows);
1606:   PetscFree(cols);
1607:   return(0);
1608: }

1612: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1613: {
1614:   PetscDS        ds;
1615:   PetscInt       **refPointFieldN;
1616:   PetscScalar    ***refPointFieldMats;
1617:   PetscInt       numFields, pRefStart, pRefEnd, p, f;
1618:   PetscSection   refConSec;

1622:   refPointFieldN = *childrenN;
1623:   *childrenN = NULL;
1624:   refPointFieldMats = *childrenMats;
1625:   *childrenMats = NULL;
1626:   DMGetDS(refTree,&ds);
1627:   PetscDSGetNumFields(ds,&numFields);
1628:   DMGetDefaultConstraints(refTree,&refConSec,NULL);
1629:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1630:   for (p = pRefStart; p < pRefEnd; p++) {
1631:     PetscInt parent, pDof;

1633:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
1634:     PetscSectionGetDof(refConSec,p,&pDof);
1635:     if (!pDof || parent == p) continue;

1637:     for (f = 0; f < numFields; f++) {
1638:       PetscInt cDof;

1640:       if (numFields > 1) {
1641:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1642:       }
1643:       else {
1644:         PetscSectionGetDof(refConSec,p,&cDof);
1645:       }

1647:       if (!cDof) continue;
1648:       PetscFree(refPointFieldMats[p - pRefStart][f]);
1649:     }
1650:     PetscFree(refPointFieldMats[p - pRefStart]);
1651:     PetscFree(refPointFieldN[p - pRefStart]);
1652:   }
1653:   PetscFree(refPointFieldMats);
1654:   PetscFree(refPointFieldN);
1655:   return(0);
1656: }

1660: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1661: {
1662:   DM             refTree;
1663:   PetscDS        ds;
1664:   Mat            refCmat;
1665:   PetscInt       numFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1666:   PetscScalar ***refPointFieldMats, *pointWork;
1667:   PetscSection   refConSec, refAnSec, anSec;
1668:   IS             refAnIS, anIS;
1669:   const PetscInt *anchors;

1674:   DMGetDS(dm,&ds);
1675:   PetscDSGetNumFields(ds,&numFields);
1676:   DMPlexGetReferenceTree(dm,&refTree);
1677:   DMSetDS(refTree,ds);
1678:   DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1679:   DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1680:   DMPlexGetAnchors(dm,&anSec,&anIS);
1681:   ISGetIndices(anIS,&anchors);
1682:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1683:   PetscSectionGetChart(conSec,&conStart,&conEnd);
1684:   PetscSectionGetMaxDof(refConSec,&maxDof);
1685:   PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1686:   PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);

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

1691:   /* step 2: compute the preorder */
1692:   DMPlexGetChart(dm,&pStart,&pEnd);
1693:   PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);
1694:   for (p = pStart; p < pEnd; p++) {
1695:     perm[p - pStart] = p;
1696:     iperm[p - pStart] = p-pStart;
1697:   }
1698:   for (p = 0; p < pEnd - pStart;) {
1699:     PetscInt point = perm[p];
1700:     PetscInt parent;

1702:     DMPlexGetTreeParent(dm,point,&parent,NULL);
1703:     if (parent == point) {
1704:       p++;
1705:     }
1706:     else {
1707:       PetscInt size, closureSize, *closure = NULL, i;

1709:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1710:       for (i = 0; i < closureSize; i++) {
1711:         PetscInt q = closure[2*i];
1712:         if (iperm[q-pStart] > iperm[point-pStart]) {
1713:           /* swap */
1714:           perm[p]               = q;
1715:           perm[iperm[q-pStart]] = point;
1716:           iperm[point-pStart]   = iperm[q-pStart];
1717:           iperm[q-pStart]       = p;
1718:           break;
1719:         }
1720:       }
1721:       size = closureSize;
1722:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1723:       if (i == size) {
1724:         p++;
1725:       }
1726:     }
1727:   }

1729:   /* step 3: fill the constraint matrix */
1730:   /* we are going to use a preorder progressive fill strategy.  Mat doesn't
1731:    * allow progressive fill without assembly, so we are going to set up the
1732:    * values outside of the Mat first.
1733:    */
1734:   {
1735:     PetscInt nRows, row, nnz;
1736:     PetscBool done;
1737:     const PetscInt *ia, *ja;
1738:     PetscScalar *vals;

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

1750:       DMPlexGetTreeParent(dm,point,&parent,&childid);
1751:       if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1752:       PetscSectionGetDof(conSec,point,&pointDof);
1753:       if (!pointDof) continue;
1754:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1755:       for (f = 0; f < numFields; f++) {
1756:         PetscInt cDof, cOff, numCols, numFillCols, i, r, fComp = -1, matOffset, offset;
1757:         PetscScalar *pointMat;
1758:         PetscObject disc;
1759:         PetscClassId id;
1760:         PetscFE fe = NULL;
1761:         PetscFV fv = NULL;

1763:         PetscDSGetDiscretization(ds,f,&disc);
1764:         PetscObjectGetClassId(disc,&id);
1765:         if (id == PETSCFE_CLASSID) {
1766:           fe = (PetscFE) disc;
1767:           PetscFEGetNumComponents(fe,&fComp);
1768:         }
1769:         else if (id == PETSCFV_CLASSID) {
1770:           fv = (PetscFV) disc;
1771:           PetscFVGetNumComponents(fv,&fComp);
1772:         }
1773:         else {
1774:           SETERRQ(PetscObjectComm((PetscObject)disc),PETSC_ERR_SUP,"Unsupported discretization");
1775:         }

1777:         if (numFields > 1) {
1778:           PetscSectionGetFieldDof(conSec,point,f,&cDof);
1779:           PetscSectionGetFieldOffset(conSec,point,f,&cOff);
1780:         }
1781:         else {
1782:           PetscSectionGetDof(conSec,point,&cDof);
1783:           PetscSectionGetOffset(conSec,point,&cOff);
1784:         }
1785:         if (!cDof) continue;

1787:         /* make sure that every row for this point is the same size */
1788: #if defined(PETSC_USE_DEBUG)
1789:         for (r = 0; r < cDof; r++) {
1790:           if (cDof > 1 && r) {
1791:             if ((ia[cOff+r+1]-ia[cOff+r]) != (ia[cOff+r]-ia[cOff+r-1])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Two point rows have different nnz: %D vs. %D", (ia[cOff+r+1]-ia[cOff+r]), (ia[cOff+r]-ia[cOff+r-1]));
1792:           }
1793:         }
1794: #endif
1795:         /* zero rows */
1796:         for (i = ia[cOff] ; i< ia[cOff+cDof];i++) {
1797:           vals[i] = 0.;
1798:         }
1799:         matOffset = ia[cOff];
1800:         numFillCols = ia[cOff+1] - matOffset;
1801:         pointMat = refPointFieldMats[childid-pRefStart][f];
1802:         numCols = refPointFieldN[childid-pRefStart][f];
1803:         offset = 0;
1804:         for (i = 0; i < closureSize; i++) {
1805:           PetscInt q = closure[2*i];
1806:           PetscInt o = closure[2*i+1];
1807:           PetscInt aDof, aOff, j, k, qConDof, qConOff;

1809:           qConDof = qConOff = 0;
1810:           if (numFields > 1) {
1811:             PetscSectionGetFieldDof(section,q,f,&aDof);
1812:             PetscSectionGetFieldOffset(section,q,f,&aOff);
1813:             if (q >= conStart && q < conEnd) {
1814:               PetscSectionGetFieldDof(conSec,q,f,&qConDof);
1815:               PetscSectionGetFieldOffset(conSec,q,f,&qConOff);
1816:             }
1817:           }
1818:           else {
1819:             PetscSectionGetDof(section,q,&aDof);
1820:             PetscSectionGetOffset(section,q,&aOff);
1821:             if (q >= conStart && q < conEnd) {
1822:               PetscSectionGetDof(conSec,q,&qConDof);
1823:               PetscSectionGetOffset(conSec,q,&qConOff);
1824:             }
1825:           }
1826:           if (!aDof) continue;
1827:           if (qConDof) {
1828:             /* this point has anchors: its rows of the matrix should already
1829:              * be filled, thanks to preordering */
1830:             /* first multiply into pointWork, then set in matrix */
1831:             PetscInt aMatOffset = ia[qConOff];
1832:             PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1833:             for (r = 0; r < cDof; r++) {
1834:               for (j = 0; j < aNumFillCols; j++) {
1835:                 PetscScalar inVal = 0;
1836:                 for (k = 0; k < aDof; k++) {
1837:                   PetscInt node = (o >= 0) ? (k / fComp) : ((aDof - 1 - k) / fComp);
1838:                   PetscInt comp = (k % fComp);
1839:                   PetscInt col  = node * fComp + comp;

1841:                   inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j];
1842:                 }
1843:                 pointWork[r * aNumFillCols + j] = inVal;
1844:               }
1845:             }
1846:             /* assume that the columns are sorted, spend less time searching */
1847:             for (j = 0, k = 0; j < aNumFillCols; j++) {
1848:               PetscInt col = ja[aMatOffset + j];
1849:               for (;k < numFillCols; k++) {
1850:                 if (ja[matOffset + k] == col) {
1851:                   break;
1852:                 }
1853:               }
1854:               if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col);
1855:               for (r = 0; r < cDof; r++) {
1856:                 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1857:               }
1858:             }
1859:           }
1860:           else {
1861:             /* find where to put this portion of pointMat into the matrix */
1862:             for (k = 0; k < numFillCols; k++) {
1863:               if (ja[matOffset + k] == aOff) {
1864:                 break;
1865:               }
1866:             }
1867:             if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff);
1868:             for (j = 0; j < aDof; j++) {
1869:               PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp);
1870:               PetscInt comp = (j % fComp);
1871:               PetscInt col  = node * fComp + comp;
1872:               for (r = 0; r < cDof; r++) {
1873:                 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + col];
1874:               }
1875:             }
1876:           }
1877:           offset += aDof;
1878:         }
1879:       }
1880:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1881:     }
1882:     for (row = 0; row < nRows; row++) {
1883:       MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);
1884:     }
1885:     MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1886:     if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix");
1887:     MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1888:     MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1889:     PetscFree(vals);
1890:   }

1892:   /* clean up */
1893:   ISRestoreIndices(anIS,&anchors);
1894:   PetscFree2(perm,iperm);
1895:   PetscFree(pointWork);
1896:   DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
1897:   return(0);
1898: }

1902: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1903:  * a non-conforming mesh.  Local refinement comes later */
1904: PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm)
1905: {
1906:   DM K;
1907:   PetscMPIInt rank;
1908:   PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1909:   PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1910:   PetscInt *Kembedding;
1911:   PetscInt *cellClosure=NULL, nc;
1912:   PetscScalar *newVertexCoords;
1913:   PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1914:   PetscSection parentSection;

1918:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1919:   DMGetDimension(dm,&dim);
1920:   DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);
1921:   DMSetDimension(*ncdm,dim);

1923:   DMPlexGetChart(dm, &pStart, &pEnd);
1924:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);
1925:   DMPlexGetReferenceTree(dm,&K);
1926:   if (!rank) {
1927:     /* compute the new charts */
1928:     PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);
1929:     offset = 0;
1930:     for (d = 0; d <= dim; d++) {
1931:       PetscInt pOldCount, kStart, kEnd, k;

1933:       pNewStart[d] = offset;
1934:       DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);
1935:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1936:       pOldCount = pOldEnd[d] - pOldStart[d];
1937:       /* adding the new points */
1938:       pNewCount[d] = pOldCount + kEnd - kStart;
1939:       if (!d) {
1940:         /* removing the cell */
1941:         pNewCount[d]--;
1942:       }
1943:       for (k = kStart; k < kEnd; k++) {
1944:         PetscInt parent;
1945:         DMPlexGetTreeParent(K,k,&parent,NULL);
1946:         if (parent == k) {
1947:           /* avoid double counting points that won't actually be new */
1948:           pNewCount[d]--;
1949:         }
1950:       }
1951:       pNewEnd[d] = pNewStart[d] + pNewCount[d];
1952:       offset = pNewEnd[d];

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

1959:     PetscMalloc1(pNewEnd[dim],&newConeSizes);
1960:     {
1961:       PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;

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

1966:       for (k = kStart; k < kEnd; k++) {
1967:         perm[k - kStart] = k;
1968:         iperm [k - kStart] = k - kStart;
1969:         preOrient[k - kStart] = 0;
1970:       }

1972:       DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1973:       for (j = 1; j < closureSizeK; j++) {
1974:         PetscInt parentOrientA = closureK[2*j+1];
1975:         PetscInt parentOrientB = cellClosure[2*j+1];
1976:         PetscInt p, q;

1978:         p = closureK[2*j];
1979:         q = cellClosure[2*j];
1980:         for (d = 0; d <= dim; d++) {
1981:           if (q >= pOldStart[d] && q < pOldEnd[d]) {
1982:             Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1983:           }
1984:         }
1985:         if (parentOrientA != parentOrientB) {
1986:           PetscInt numChildren, i;
1987:           const PetscInt *children;

1989:           DMPlexGetTreeChildren(K,p,&numChildren,&children);
1990:           for (i = 0; i < numChildren; i++) {
1991:             PetscInt kPerm, oPerm;

1993:             k    = children[i];
1994:             DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);
1995:             /* perm = what refTree position I'm in */
1996:             perm[kPerm-kStart]      = k;
1997:             /* iperm = who is at this position */
1998:             iperm[k-kStart]         = kPerm-kStart;
1999:             preOrient[kPerm-kStart] = oPerm;
2000:           }
2001:         }
2002:       }
2003:       DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
2004:     }
2005:     PetscSectionSetChart(parentSection,0,pNewEnd[dim]);
2006:     offset = 0;
2007:     numNewCones = 0;
2008:     for (d = 0; d <= dim; d++) {
2009:       PetscInt kStart, kEnd, k;
2010:       PetscInt p;
2011:       PetscInt size;

2013:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
2014:         /* skip cell 0 */
2015:         if (p == cell) continue;
2016:         /* old cones to new cones */
2017:         DMPlexGetConeSize(dm,p,&size);
2018:         newConeSizes[offset++] = size;
2019:         numNewCones += size;
2020:       }

2022:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
2023:       for (k = kStart; k < kEnd; k++) {
2024:         PetscInt kParent;

2026:         DMPlexGetTreeParent(K,k,&kParent,NULL);
2027:         if (kParent != k) {
2028:           Kembedding[k] = offset;
2029:           DMPlexGetConeSize(K,k,&size);
2030:           newConeSizes[offset++] = size;
2031:           numNewCones += size;
2032:           if (kParent != 0) {
2033:             PetscSectionSetDof(parentSection,Kembedding[k],1);
2034:           }
2035:         }
2036:       }
2037:     }

2039:     PetscSectionSetUp(parentSection);
2040:     PetscSectionGetStorageSize(parentSection,&numPointsWithParents);
2041:     PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);
2042:     PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);

2044:     /* fill new cones */
2045:     offset = 0;
2046:     for (d = 0; d <= dim; d++) {
2047:       PetscInt kStart, kEnd, k, l;
2048:       PetscInt p;
2049:       PetscInt size;
2050:       const PetscInt *cone, *orientation;

2052:       for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
2053:         /* skip cell 0 */
2054:         if (p == cell) continue;
2055:         /* old cones to new cones */
2056:         DMPlexGetConeSize(dm,p,&size);
2057:         DMPlexGetCone(dm,p,&cone);
2058:         DMPlexGetConeOrientation(dm,p,&orientation);
2059:         for (l = 0; l < size; l++) {
2060:           newCones[offset]          = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
2061:           newOrientations[offset++] = orientation[l];
2062:         }
2063:       }

2065:       DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
2066:       for (k = kStart; k < kEnd; k++) {
2067:         PetscInt kPerm = perm[k], kParent;
2068:         PetscInt preO  = preOrient[k];

2070:         DMPlexGetTreeParent(K,k,&kParent,NULL);
2071:         if (kParent != k) {
2072:           /* embed new cones */
2073:           DMPlexGetConeSize(K,k,&size);
2074:           DMPlexGetCone(K,kPerm,&cone);
2075:           DMPlexGetConeOrientation(K,kPerm,&orientation);
2076:           for (l = 0; l < size; l++) {
2077:             PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size);
2078:             PetscInt newO, lSize, oTrue;

2080:             q                         = iperm[cone[m]];
2081:             newCones[offset]          = Kembedding[q];
2082:             DMPlexGetConeSize(K,q,&lSize);
2083:             oTrue                     = orientation[m];
2084:             oTrue                     = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
2085:             newO                      = DihedralCompose(lSize,oTrue,preOrient[q]);
2086:             newOrientations[offset++] = newO;
2087:           }
2088:           if (kParent != 0) {
2089:             PetscInt newPoint = Kembedding[kParent];
2090:             PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);
2091:             parents[pOffset]  = newPoint;
2092:             childIDs[pOffset] = k;
2093:           }
2094:         }
2095:       }
2096:     }

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

2100:     /* fill coordinates */
2101:     offset = 0;
2102:     {
2103:       PetscInt kStart, kEnd, l;
2104:       PetscSection vSection;
2105:       PetscInt v;
2106:       Vec coords;
2107:       PetscScalar *coordvals;
2108:       PetscInt dof, off;
2109:       PetscReal v0[3], J[9], detJ;

2111: #if defined(PETSC_USE_DEBUG)
2112:       {
2113:         PetscInt k;
2114:         DMPlexGetHeightStratum(K,0,&kStart,&kEnd);
2115:         for (k = kStart; k < kEnd; k++) {
2116:           DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);
2117:           if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k);
2118:         }
2119:       }
2120: #endif
2121:       DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);
2122:       DMGetCoordinateSection(dm,&vSection);
2123:       DMGetCoordinatesLocal(dm,&coords);
2124:       VecGetArray(coords,&coordvals);
2125:       for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {

2127:         PetscSectionGetDof(vSection,v,&dof);
2128:         PetscSectionGetOffset(vSection,v,&off);
2129:         for (l = 0; l < dof; l++) {
2130:           newVertexCoords[offset++] = coordvals[off + l];
2131:         }
2132:       }
2133:       VecRestoreArray(coords,&coordvals);

2135:       DMGetCoordinateSection(K,&vSection);
2136:       DMGetCoordinatesLocal(K,&coords);
2137:       VecGetArray(coords,&coordvals);
2138:       DMPlexGetDepthStratum(K,0,&kStart,&kEnd);
2139:       for (v = kStart; v < kEnd; v++) {
2140:         PetscReal coord[3], newCoord[3];
2141:         PetscInt  vPerm = perm[v];
2142:         PetscInt  kParent;

2144:         DMPlexGetTreeParent(K,v,&kParent,NULL);
2145:         if (kParent != v) {
2146:           /* this is a new vertex */
2147:           PetscSectionGetOffset(vSection,vPerm,&off);
2148:           for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]);
2149:           CoordinatesRefToReal(dim, dim, v0, J, coord, newCoord);
2150:           for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l];
2151:           offset += dim;
2152:         }
2153:       }
2154:       VecRestoreArray(coords,&coordvals);
2155:     }

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

2161:       tmp = pNewCount[d];
2162:       pNewCount[d] = pNewCount[dim - d];
2163:       pNewCount[dim - d] = tmp;
2164:     }

2166:     DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);
2167:     DMPlexSetReferenceTree(*ncdm,K);
2168:     DMPlexSetTree(*ncdm,parentSection,parents,childIDs);

2170:     /* clean up */
2171:     DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
2172:     PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);
2173:     PetscFree(newConeSizes);
2174:     PetscFree2(newCones,newOrientations);
2175:     PetscFree(newVertexCoords);
2176:     PetscFree2(parents,childIDs);
2177:     PetscFree4(Kembedding,perm,iperm,preOrient);
2178:   }
2179:   else {
2180:     PetscInt    p, counts[4];
2181:     PetscInt    *coneSizes, *cones, *orientations;
2182:     Vec         coordVec;
2183:     PetscScalar *coords;

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

2188:       DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
2189:       counts[d] = dEnd - dStart;
2190:     }
2191:     PetscMalloc1(pEnd-pStart,&coneSizes);
2192:     for (p = pStart; p < pEnd; p++) {
2193:       DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);
2194:     }
2195:     DMPlexGetCones(dm, &cones);
2196:     DMPlexGetConeOrientations(dm, &orientations);
2197:     DMGetCoordinatesLocal(dm,&coordVec);
2198:     VecGetArray(coordVec,&coords);

2200:     PetscSectionSetChart(parentSection,pStart,pEnd);
2201:     PetscSectionSetUp(parentSection);
2202:     DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);
2203:     DMPlexSetReferenceTree(*ncdm,K);
2204:     DMPlexSetTree(*ncdm,parentSection,NULL,NULL);
2205:     VecRestoreArray(coordVec,&coords);
2206:   }
2207:   PetscSectionDestroy(&parentSection);

2209:   return(0);
2210: }

2212: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt[],PetscBool,PetscInt,PetscInt[]);
2213: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt[]);

2217: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2218: {
2219:   PetscSF        coarseToFineEmbedded;
2220:   PetscSection   globalCoarse, globalFine;
2221:   PetscSection   localCoarse, localFine;
2222:   PetscSection   aSec, cSec;
2223:   PetscSection   rootIndicesSec, rootMatricesSec;
2224:   PetscSection   leafIndicesSec, leafMatricesSec;
2225:   PetscInt       *rootIndices, *leafIndices;
2226:   PetscScalar    *rootMatrices, *leafMatrices;
2227:   IS             aIS;
2228:   const PetscInt *anchors;
2229:   Mat            cMat;
2230:   PetscInt       numFields;
2231:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
2232:   PetscInt       aStart, aEnd, cStart, cEnd;
2233:   PetscInt       *maxChildIds;
2234:   PetscInt       *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;

2238:   DMPlexGetChart(coarse,&pStartC,&pEndC);
2239:   DMPlexGetChart(fine,&pStartF,&pEndF);
2240:   DMGetDefaultGlobalSection(fine,&globalFine);
2241:   { /* winnow fine points that don't have global dofs out of the sf */
2242:     PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;

2244:     for (p = pStartF, numPointsWithDofs = 0; p < pEndF; p++) {
2245:       PetscSectionGetDof(globalFine,p,&dof);
2246:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
2247:       if ((dof - cdof) > 0) {
2248:         numPointsWithDofs++;
2249:       }
2250:     }
2251:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
2252:     for (p = pStartF, offset = 0; p < pEndF; p++) {
2253:       PetscSectionGetDof(globalFine,p,&dof);
2254:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
2255:       if ((dof - cdof) > 0) {
2256:         pointsWithDofs[offset++] = p - pStartF;
2257:       }
2258:     }
2259:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
2260:     PetscFree(pointsWithDofs);
2261:   }
2262:   /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2263:   PetscMalloc1(pEndC-pStartC,&maxChildIds);
2264:   for (p = pStartC; p < pEndC; p++) {
2265:     maxChildIds[p - pStartC] = -2;
2266:   }
2267:   PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);
2268:   PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);

2270:   DMGetDefaultSection(coarse,&localCoarse);
2271:   DMGetDefaultGlobalSection(coarse,&globalCoarse);

2273:   DMPlexGetAnchors(coarse,&aSec,&aIS);
2274:   ISGetIndices(aIS,&anchors);
2275:   PetscSectionGetChart(aSec,&aStart,&aEnd);

2277:   DMGetDefaultConstraints(coarse,&cSec,&cMat);
2278:   PetscSectionGetChart(cSec,&cStart,&cEnd);

2280:   /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2281:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
2282:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);
2283:   PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);
2284:   PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);
2285:   PetscSectionGetNumFields(globalCoarse,&numFields);
2286:   {
2287:     PetscInt maxFields = PetscMax(1,numFields) + 1;
2288:     PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);
2289:   }

2291:   for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2292:     PetscInt dof, matSize   = 0;
2293:     PetscInt aDof           = 0;
2294:     PetscInt cDof           = 0;
2295:     PetscInt maxChildId     = maxChildIds[p - pStartC];
2296:     PetscInt numRowIndices  = 0;
2297:     PetscInt numColIndices  = 0;

2299:     PetscSectionGetDof(globalCoarse,p,&dof);
2300:     if (dof < 0) {
2301:       dof = -(dof + 1);
2302:     }
2303:     if (p >= aStart && p < aEnd) {
2304:       PetscSectionGetDof(aSec,p,&aDof);
2305:     }
2306:     if (p >= cStart && p < cEnd) {
2307:       PetscSectionGetDof(cSec,p,&cDof);
2308:     }
2309:     offsets[0]    = 0;
2310:     newOffsets[0] = 0;
2311:     if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2312:       PetscInt *closure = NULL, closureSize, cl, f;

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

2318:         PetscSectionGetDof(localCoarse,c,&clDof);
2319:         numRowIndices += clDof;
2320:         for (f = 0; f < numFields; f++) {
2321:           PetscSectionGetFieldDof(localCoarse,c,f,&clDof);
2322:           offsets[f + 1] += clDof;
2323:         }
2324:       }
2325:       for (f = 0; f < numFields; f++) {
2326:         offsets[f + 1]   += offsets[f];
2327:         newOffsets[f + 1] = offsets[f + 1];
2328:       }
2329:       /* get the number of indices needed and their field offsets */
2330:       DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);
2331:       DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
2332:       if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2333:         numColIndices = numRowIndices;
2334:         matSize = 0;
2335:       }
2336:       else if (numFields) { /* we send one submat for each field: sum their sizes */
2337:         matSize = 0;
2338:         for (f = 0; f < numFields; f++) {
2339:           PetscInt numRow, numCol;

2341:           numRow = offsets[f + 1] - offsets[f];
2342:           numCol = newOffsets[f + 1] - newOffsets[f + 1];
2343:           matSize += numRow * numCol;
2344:         }
2345:       }
2346:       else {
2347:         matSize = numRowIndices * numColIndices;
2348:       }
2349:     }
2350:     else if (maxChildId == -1) {
2351:       if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2352:         PetscInt aOff, a, f;

2354:         PetscSectionGetOffset(aSec,p,&aOff);
2355:         for (f = 0; f < numFields; f++) {
2356:           PetscInt fDof;

2358:           PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2359:           offsets[f+1] = fDof;
2360:         }
2361:         for (a = 0; a < aDof; a++) {
2362:           PetscInt anchor = anchors[a + aOff], aLocalDof;

2364:           PetscSectionGetDof(localCoarse,anchor,&aLocalDof);
2365:           numColIndices += aLocalDof;
2366:           for (f = 0; f < numFields; f++) {
2367:             PetscInt fDof;

2369:             PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2370:             newOffsets[f+1] += fDof;
2371:           }
2372:         }
2373:         if (numFields) {
2374:           matSize = 0;
2375:           for (f = 0; f < numFields; f++) {
2376:             matSize += offsets[f+1] * newOffsets[f+1];
2377:           }
2378:         }
2379:         else {
2380:           matSize = numColIndices * dof;
2381:         }
2382:       }
2383:       else { /* no children, and no constraints on dofs: just get the global indices */
2384:         numColIndices = dof;
2385:         matSize       = 0;
2386:       }
2387:     }
2388:     /* we will pack the column indices with the field offsets */
2389:     PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);
2390:     PetscSectionSetDof(rootMatricesSec,p,matSize);
2391:   }
2392:   PetscSectionSetUp(rootIndicesSec);
2393:   PetscSectionSetUp(rootMatricesSec);
2394:   {
2395:     PetscInt    numRootIndices, numRootMatrices;

2397:     PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
2398:     PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);
2399:     PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);
2400:     for (p = pStartC; p < pEndC; p++) {
2401:       PetscInt    numRowIndices, numColIndices, matSize, dof;
2402:       PetscInt    pIndOff, pMatOff;
2403:       PetscInt    *pInd;
2404:       PetscInt    maxChildId = maxChildIds[p - pStartC];
2405:       PetscScalar *pMat = NULL;

2407:       PetscSectionGetDof(rootIndicesSec,p,&numColIndices);
2408:       if (!numColIndices) {
2409:         continue;
2410:       }
2411:       offsets[0]        = 0;
2412:       newOffsets[0]     = 0;
2413:       offsetsCopy[0]    = 0;
2414:       newOffsetsCopy[0] = 0;
2415:       numColIndices -= 2 * numFields;
2416:       PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);
2417:       pInd = &(rootIndices[pIndOff]);
2418:       PetscSectionGetDof(rootMatricesSec,p,&matSize);
2419:       if (matSize) {
2420:         PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);
2421:         pMat = &rootMatrices[pMatOff];
2422:       }
2423:       PetscSectionGetDof(globalCoarse,p,&dof);
2424:       if (dof < 0) {
2425:         dof = -(dof + 1);
2426:       }
2427:       if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2428:         PetscInt i, j;
2429:         PetscInt numRowIndices = matSize / numColIndices;

2431:         if (!numRowIndices) { /* don't need to calculate the mat, just the indices */
2432:           PetscInt numIndices, *indices;
2433:           DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2434:           if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations");
2435:           for (i = 0; i < numColIndices; i++) {
2436:             pInd[i] = indices[i];
2437:           }
2438:           for (i = 0; i < numFields; i++) {
2439:             pInd[numColIndices + i]             = offsets[i+1];
2440:             pInd[numColIndices + numFields + i] = offsets[i+1];
2441:           }
2442:           DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2443:         }
2444:         else {
2445:           PetscInt closureSize, *closure = NULL, cl;
2446:           PetscScalar *pMatIn, *pMatModified;
2447:           PetscInt numPoints,*points;

2449:           DMGetWorkArray(coarse,numRowIndices * numRowIndices,PETSC_SCALAR,&pMatIn);
2450:           for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2451:             for (j = 0; j < numRowIndices; j++) {
2452:               pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2453:             }
2454:           }
2455:           DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2456:           if (numFields) {
2457:             PetscInt f;

2459:             for (cl = 0; cl < closureSize; cl++) {
2460:               PetscInt c = closure[2 * cl];

2462:               for (f = 0; f < numFields; f++) {
2463:                 PetscInt fDof;

2465:                 PetscSectionGetFieldDof(localCoarse,c,f,&fDof);
2466:                 offsets[f + 1] += fDof;
2467:               }
2468:             }
2469:             for (f = 0; f < numFields; f++) {
2470:               offsets[f + 1]   += offsets[f];
2471:               newOffsets[f + 1] = offsets[f + 1];
2472:             }
2473:           }
2474:           /* apply hanging node constraints on the right, get the new points and the new offsets */
2475:           DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);
2476:           if (!numFields) {
2477:             for (i = 0; i < numRowIndices * numColIndices; i++) {
2478:               pMat[i] = pMatModified[i];
2479:             }
2480:           }
2481:           else {
2482:             PetscInt f, i, j, count;
2483:             for (f = 0, count = 0; f < numFields; f++) {
2484:               for (i = offsets[f]; i < offsets[f+1]; i++) {
2485:                 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) {
2486:                   pMat[count] = pMatModified[i * numColIndices + j];
2487:                 }
2488:               }
2489:             }
2490:           }
2491:           DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatModified);
2492:           DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2493:           DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatIn);
2494:           if (numFields) {
2495:             PetscInt f;
2496:             for (f = 0; f < numFields; f++) {
2497:               pInd[numColIndices + f]             = offsets[f+1];
2498:               pInd[numColIndices + numFields + f] = newOffsets[f+1];
2499:             }
2500:             for (cl = 0; cl < numPoints*2; cl += 2) {
2501:               PetscInt o = points[cl+1], globalOff, c = points[cl];
2502:               PetscSectionGetOffset(globalCoarse, c, &globalOff);
2503:               indicesPointFields_private(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, o, pInd);
2504:             }
2505:           } else {
2506:             for (cl = 0; cl < numPoints*2; cl += 2) {
2507:               PetscInt o = points[cl+1], c = points[cl], globalOff;
2508:               PetscSectionGetOffset(globalCoarse, c, &globalOff);
2509:               indicesPoint_private(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, o, pInd);
2510:             }
2511:           }
2512:           DMRestoreWorkArray(coarse,numPoints,PETSC_SCALAR,&points);
2513:         }
2514:       }
2515:       else if (matSize) {
2516:         PetscInt cOff;
2517:         PetscInt *rowIndices, *colIndices, a, aDof, aOff;

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

2529:           for (f = 0; f < numFields; f++) {
2530:             PetscInt fDof;
2531:             PetscSectionGetFieldDof(cSec,p,f,&fDof);
2532:             offsets[f + 1] = fDof;
2533:             for (a = 0; a < aDof; a++) {
2534:               PetscInt anchor = anchors[a + aOff];
2535:               PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2536:               newOffsets[f + 1] += fDof;
2537:             }
2538:           }
2539:           for (f = 0; f < numFields; f++) {
2540:             offsets[f + 1]       += offsets[f];
2541:             offsetsCopy[f + 1]    = offsets[f + 1];
2542:             newOffsets[f + 1]    += newOffsets[f];
2543:             newOffsetsCopy[f + 1] = newOffsets[f + 1];
2544:           }
2545:           indicesPointFields_private(cSec,p,cOff,offsetsCopy,PETSC_TRUE,0,rowIndices);
2546:           for (a = 0; a < aDof; a++) {
2547:             PetscInt anchor = anchors[a + aOff], lOff;
2548:             PetscSectionGetOffset(localCoarse,anchor,&lOff);
2549:             indicesPointFields_private(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,0,colIndices);
2550:           }
2551:         }
2552:         else {
2553:           indicesPoint_private(cSec,p,cOff,offsetsCopy,PETSC_TRUE,0,rowIndices);
2554:           for (a = 0; a < aDof; a++) {
2555:             PetscInt anchor = anchors[a + aOff], lOff;
2556:             PetscSectionGetOffset(localCoarse,anchor,&lOff);
2557:             indicesPoint_private(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,0,colIndices);
2558:           }
2559:         }
2560:         if (numFields) {
2561:           PetscInt count, f, a;
2562:           for (f = 0, count = 0; f < numFields; f++) {
2563:             PetscInt iSize = offsets[f + 1] - offsets[f];
2564:             PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2565:             MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);
2566:             count += iSize * jSize;
2567:             pInd[numColIndices + f]             = offsets[f+1];
2568:             pInd[numColIndices + numFields + f] = newOffsets[f+1];
2569:           }
2570:           for (a = 0; a < aDof; a++) {
2571:             PetscInt anchor = anchors[a + aOff];
2572:             PetscInt gOff;
2573:             PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2574:             indicesPointFields_private(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,0,pInd);
2575:           }
2576:         }
2577:         else {
2578:           PetscInt a;
2579:           MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);
2580:           for (a = 0; a < aDof; a++) {
2581:             PetscInt anchor = anchors[a + aOff];
2582:             PetscInt gOff;
2583:             PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2584:             indicesPoint_private(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,0,pInd);
2585:           }
2586:         }
2587:         DMRestoreWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);
2588:         DMRestoreWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);
2589:       }
2590:       else {
2591:         PetscInt gOff;

2593:         PetscSectionGetOffset(globalCoarse,p,&gOff);
2594:         if (numFields) {
2595:           PetscInt f;

2597:           for (f = 0; f < numFields; f++) {
2598:             PetscInt fDof;
2599:             PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2600:             offsets[f + 1] = fDof + offsets[f];
2601:           }
2602:           for (f = 0; f < numFields; f++) {
2603:             pInd[numColIndices + f]             = offsets[f+1];
2604:             pInd[numColIndices + numFields + f] = offsets[f+1];
2605:           }
2606:           indicesPointFields_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);
2607:         }
2608:         else {
2609:           indicesPoint_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);
2610:         }
2611:       }
2612:     }
2613:     PetscFree(maxChildIds);
2614:   }
2615:   {
2616:     PetscSF  indicesSF, matricesSF;
2617:     PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;

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

2657:     PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);
2658:     PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);
2659:     MatGetLayouts(mat,&rowMap,&colMap);
2660:     PetscLayoutSetUp(rowMap);
2661:     PetscLayoutSetUp(colMap);
2662:     PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
2663:     PetscLayoutGetRange(colMap,&colStart,&colEnd);
2664:     PetscSectionGetMaxDof(globalFine,&maxDof);
2665:     DMGetWorkArray(fine,maxDof,PETSC_INT,&rowIndices);
2666:     for (p = pStartF; p < pEndF; p++) {
2667:       PetscInt    gDof, gcDof, gOff;
2668:       PetscInt    numColIndices, pIndOff, *pInd;
2669:       PetscInt    matSize;
2670:       PetscInt    i;

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

2694:           PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2695:           offsets[f + 1]        = offsets[f] + rowDof;
2696:           offsetsCopy[f + 1]    = offsets[f + 1];
2697:           newOffsets[f + 1]     = pInd[numColIndices + numFields + f];
2698:           numD[f] = 0;
2699:           numO[f] = 0;
2700:         }
2701:         indicesPointFields_private(localFine,p,gOff,offsetsCopy,PETSC_FALSE,0,rowIndices);
2702:         for (f = 0; f < numFields; f++) {
2703:           PetscInt colOffset    = newOffsets[f];
2704:           PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];

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

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

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

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

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


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

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

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

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

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

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

3045:     if (canonical) {
3046:       PetscInt pCanonical;
3047:       DMLabelGetValue(canonical,p,&pCanonical);
3048:       if (p != pCanonical) continue;
3049:     }
3050:     DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3051:     if (!numChildren) continue;
3052:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
3053:       PetscInt child = children[i];
3054:       PetscInt dof;

3056:       PetscSectionGetDof(section,child,&dof);
3057:       numChildDof += dof;
3058:     }
3059:     PetscSectionGetDof(section,p,&numSelfDof);
3060:     if (!numChildDof || !numSelfDof) continue;
3061:     for (f = 0; f < numFields; f++) {
3062:       PetscInt selfOff;

3064:       if (numSecFields) { /* count the dofs for just this field */
3065:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
3066:           PetscInt child = children[i];
3067:           PetscInt dof;

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

3091:     /* same conditions about when entries occur */
3092:     if (canonical) {
3093:       PetscInt pCanonical;
3094:       DMLabelGetValue(canonical,p,&pCanonical);
3095:       if (p != pCanonical) continue;
3096:     }
3097:     DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3098:     if (!numChildren) continue;
3099:     for (i = 0, numChildDof = 0; i < numChildren; i++) {
3100:       PetscInt child = children[i];
3101:       PetscInt dof;

3103:       PetscSectionGetDof(section,child,&dof);
3104:       numChildDof += dof;
3105:     }
3106:     PetscSectionGetDof(section,p,&numSelfDof);
3107:     if (!numChildDof || !numSelfDof) continue;

3109:     for (f = 0; f < numFields; f++) {
3110:       PetscInt       selfOff, fComp, numSelfShapes, numChildShapes, parentCell;
3111:       PetscInt       cellShapeOff;
3112:       PetscObject    disc;
3113:       PetscDualSpace dsp;
3114:       PetscClassId   classId;
3115:       PetscScalar    *pointMat;
3116:       PetscInt       *matRows, *matCols;
3117:       PetscInt       pO = PETSC_MIN_INT;
3118:       const PetscInt *depthNumDof;

3120:       if (numSecFields) {
3121:         for (i = 0, numChildDof = 0; i < numChildren; i++) {
3122:           PetscInt child = children[i];
3123:           PetscInt dof;

3125:           PetscSectionGetFieldDof(section,child,f,&dof);
3126:           numChildDof += dof;
3127:         }
3128:         PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3129:         PetscSectionGetFieldOffset(section,p,f,&selfOff);
3130:         PetscSectionGetFieldComponents(section,f,&fComp);
3131:       }
3132:       else {
3133:         PetscSectionGetOffset(section,p,&selfOff);
3134:         fComp = 1;
3135:       }
3136:       numSelfShapes  = numSelfDof  / fComp;
3137:       numChildShapes = numChildDof / fComp;

3139:       /* find a cell whose closure contains p */
3140:       if (p >= cStart && p < cEnd) {
3141:         parentCell = p;
3142:       }
3143:       else {
3144:         PetscInt *star = NULL;
3145:         PetscInt numStar;

3147:         parentCell = -1;
3148:         DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3149:         for (i = numStar - 1; i >= 0; i--) {
3150:           PetscInt c = star[2 * i];

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

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

3180:           pO = closure[2 * i + 1];
3181:           if (point == p) break;
3182:           DMLabelGetValue(depth,point,&pointDepth);
3183:           cellShapeOff += depthNumDof[pointDepth];
3184:         }
3185:         DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);
3186:       }

3188:       DMGetWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);
3189:       DMGetWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);
3190:       matCols = matRows + numSelfShapes;
3191:       for (i = 0; i < numSelfShapes; i++) {
3192:         matRows[i] = selfOff + i * fComp;
3193:       }
3194:       {
3195:         PetscInt colOff = 0;

3197:         for (i = 0; i < numChildren; i++) {
3198:           PetscInt child = children[i];
3199:           PetscInt dof, off, j;

3201:           if (numSecFields) {
3202:             PetscSectionGetFieldDof(cSection,child,f,&dof);
3203:             PetscSectionGetFieldOffset(cSection,child,f,&off);
3204:           }
3205:           else {
3206:             PetscSectionGetDof(cSection,child,&dof);
3207:             PetscSectionGetOffset(cSection,child,&off);
3208:           }

3210:           for (j = 0; j < dof / fComp; j++) {
3211:             matCols[colOff++] = off + j * fComp;
3212:           }
3213:         }
3214:       }
3215:       if (classId == PETSCFE_CLASSID) {
3216:         PetscFE        fe = (PetscFE) disc;
3217:         PetscInt       fSize;

3219:         PetscFEGetDualSpace(fe,&dsp);
3220:         PetscDualSpaceGetDimension(dsp,&fSize);
3221:         for (i = 0; i < numSelfShapes; i++) { /* for every shape function */
3222:           PetscQuadrature q;
3223:           PetscInt        dim, numPoints, j, k;
3224:           const PetscReal *points;
3225:           const PetscReal *weights;
3226:           PetscInt        *closure = NULL;
3227:           PetscInt        numClosure;
3228:           PetscInt        parentCellShapeDof = cellShapeOff + (pO < 0 ? (numSelfShapes - 1 - i) : i);
3229:           PetscReal       *Bparent;

3231:           PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);
3232:           PetscQuadratureGetData(q,&dim,&numPoints,&points,&weights);
3233:           PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3234:           for (k = 0; k < numChildShapes; k++) {
3235:             pointMat[numChildShapes * i + k] = 0.;
3236:           }
3237:           for (j = 0; j < numPoints; j++) {
3238:             PetscInt          childCell = -1;
3239:             PetscReal         parentValAtPoint;
3240:             const PetscReal   *pointReal = &points[dim * j];
3241:             const PetscScalar *point;
3242:             PetscReal         *Bchild;
3243:             PetscInt          childCellShapeOff, pointMatOff;
3244: #if defined(PETSC_USE_COMPLEX)
3245:             PetscInt          d;

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

3255:             parentValAtPoint = Bparent[(fSize * j + parentCellShapeDof) * fComp];

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

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

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

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

3285:               DMLabelGetValue(depth,child,&childDepth);
3286:               childDof = depthNumDof[childDepth];
3287:               for (l = 0, childCellShapeOff = 0; l < numClosure; l++) {
3288:                 PetscInt point = closure[2 * l];
3289:                 PetscInt pointDepth;

3291:                 childO = closure[2 * l + 1];
3292:                 if (point == child) break;
3293:                 DMLabelGetValue(depth,point,&pointDepth);
3294:                 childCellShapeOff += depthNumDof[pointDepth];
3295:               }
3296:               if (l == numClosure) {
3297:                 pointMatOff += childDof;
3298:                 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3299:               }
3300:               for (l = 0; l < childDof; l++) {
3301:                 PetscInt    childCellDof    = childCellShapeOff + (childO ? (childDof - 1 - l) : l);
3302:                 PetscReal   childValAtPoint = Bchild[childCellDof * fComp];

3304:                 pointMat[i * numChildShapes + pointMatOff + l] += weights[j] * parentValAtPoint * childValAtPoint;
3305:               }
3306:               pointMatOff += childDof;
3307:             }
3308:             DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3309:             PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);
3310:           }
3311:           PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);
3312:         }
3313:       }
3314:       else { /* just the volume-weighted averages of the children */
3315:         PetscInt  childShapeOff;
3316:         PetscReal parentVol;

3318:         DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);
3319:         for (i = 0, childShapeOff = 0; i < numChildren; i++) {
3320:           PetscInt  child = children[i];
3321:           PetscReal childVol;

3323:           if (child < cStart || child >= cEnd) continue;
3324:           DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);
3325:           pointMat[childShapeOff] = childVol / parentVol;
3326:           childShapeOff++;
3327:         }
3328:       }
3329:       /* Insert pointMat into mat */
3330:       for (i = 0; i < fComp; i++) {
3331:         PetscInt j;
3332:         MatSetValues(mat,numSelfShapes,matRows,numChildShapes,matCols,pointMat,INSERT_VALUES);

3334:         for (j = 0; j < numSelfShapes; j++) {
3335:           matRows[j]++;
3336:         }
3337:         for (j = 0; j < numChildShapes; j++) {
3338:           matCols[j]++;
3339:         }
3340:       }
3341:       DMRestoreWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);
3342:       DMRestoreWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);
3343:     }
3344:   }
3345:   PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);
3346:   PetscFree2(pointScalar,pointRef);
3347:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3348:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3349:   *inj = mat;
3350:   return(0);
3351: }

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

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

3376:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
3377:     PetscSectionGetDof(refConSec,p,&pDof);
3378:     PetscSectionGetDof(refSection,parent,&parentDof);
3379:     if (!pDof || !parentDof || parent == p) continue;

3381:     PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);
3382:     for (f = 0; f < numFields; f++) {
3383:       PetscInt cDof, cOff, numCols, r, fComp;
3384:       PetscObject disc;
3385:       PetscClassId id;
3386:       PetscFE fe = NULL;
3387:       PetscFV fv = NULL;

3389:       PetscDSGetDiscretization(ds,f,&disc);
3390:       PetscObjectGetClassId(disc,&id);
3391:       if (id == PETSCFE_CLASSID) {
3392:         fe = (PetscFE) disc;
3393:         PetscFEGetNumComponents(fe,&fComp);
3394:       }
3395:       else if (id == PETSCFV_CLASSID) {
3396:         fv = (PetscFV) disc;
3397:         PetscFVGetNumComponents(fv,&fComp);
3398:       }
3399:       else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);

3401:       if (numFields > 1) {
3402:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3403:         PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
3404:       }
3405:       else {
3406:         PetscSectionGetDof(refConSec,p,&cDof);
3407:         PetscSectionGetOffset(refConSec,p,&cOff);
3408:       }

3410:       if (!cDof) continue;
3411:       for (r = 0; r < cDof; r++) {
3412:         rows[r] = cOff + r;
3413:       }
3414:       numCols = 0;
3415:       {
3416:         PetscInt aDof, aOff, j;

3418:         if (numFields > 1) {
3419:           PetscSectionGetFieldDof(refSection,parent,f,&aDof);
3420:           PetscSectionGetFieldOffset(refSection,parent,f,&aOff);
3421:         }
3422:         else {
3423:           PetscSectionGetDof(refSection,parent,&aDof);
3424:           PetscSectionGetOffset(refSection,parent,&aOff);
3425:         }

3427:         for (j = 0; j < aDof; j++) {
3428:           cols[numCols++] = aOff + j;
3429:         }
3430:       }
3431:       PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
3432:       /* transpose of constraint matrix */
3433:       MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);
3434:     }
3435:   }
3436:   *childrenMats = refPointFieldMats;
3437:   PetscFree(rows);
3438:   PetscFree(cols);
3439:   return(0);
3440: }

3444: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3445: {
3446:   PetscDS        ds;
3447:   PetscScalar    ***refPointFieldMats;
3448:   PetscInt       numFields, pRefStart, pRefEnd, p, f;
3449:   PetscSection   refConSec, refSection;

3453:   refPointFieldMats = *childrenMats;
3454:   *childrenMats = NULL;
3455:   DMGetDS(refTree,&ds);
3456:   DMGetDefaultSection(refTree,&refSection);
3457:   PetscDSGetNumFields(ds,&numFields);
3458:   DMGetDefaultConstraints(refTree,&refConSec,NULL);
3459:   PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3460:   for (p = pRefStart; p < pRefEnd; p++) {
3461:     PetscInt parent, pDof, parentDof;

3463:     DMPlexGetTreeParent(refTree,p,&parent,NULL);
3464:     PetscSectionGetDof(refConSec,p,&pDof);
3465:     PetscSectionGetDof(refSection,parent,&parentDof);
3466:     if (!pDof || !parentDof || parent == p) continue;

3468:     for (f = 0; f < numFields; f++) {
3469:       PetscInt cDof;

3471:       if (numFields > 1) {
3472:         PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3473:       }
3474:       else {
3475:         PetscSectionGetDof(refConSec,p,&cDof);
3476:       }

3478:       if (!cDof) continue;
3479:       PetscFree(refPointFieldMats[p - pRefStart][f]);
3480:     }
3481:     PetscFree(refPointFieldMats[p - pRefStart]);
3482:   }
3483:   PetscFree(refPointFieldMats);
3484:   return(0);
3485: }

3489: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3490: {
3491:   DM             refTree;
3492:   PetscSF        coarseToFineEmbedded;
3493:   const PetscInt *rootDegrees;
3494:   PetscSection   multiRootSec, rootIndicesSec, leafIndicesSec;
3495:   PetscSection   globalCoarse, globalFine;
3496:   PetscSection   localCoarse, localFine;
3497:   PetscSection   cSecRef;
3498:   PetscInt       *rootIndices, *leafIndices, *parentIndices, pRefStart, pRefEnd;
3499:   Mat            cMatRef, injRef;
3500:   PetscInt       numFields, numMulti, maxDof;
3501:   PetscInt       pStartC, pEndC, pStartF, pEndF, p;
3502:   PetscInt       *offsets, *offsetsCopy, *rowOffsets;
3503:   PetscLayout    rowMap, colMap;
3504:   PetscInt       rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3505:   PetscObject    injRefObj;
3506:   PetscScalar    ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */


3511:   /* get the templates for the fine-to-coarse injection from the reference tree */
3512:   DMPlexGetReferenceTree(coarse,&refTree);
3513:   DMGetDefaultConstraints(refTree,&cSecRef,&cMatRef);
3514:   PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
3515:   PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);
3516:   injRef = (Mat) injRefObj;
3517:   if (!injRef) {
3518:     DMPlexComputeInjectorReferenceTree(refTree,&injRef);
3519:     PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)injRef);
3520:     /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3521:     PetscObjectDereference((PetscObject)injRef);
3522:   }

3524:   DMPlexGetChart(fine,&pStartF,&pEndF);
3525:   DMGetDefaultSection(fine,&localFine);
3526:   DMGetDefaultGlobalSection(fine,&globalFine);

3528:   PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
3529:   PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);

3531:   PetscSectionGetNumFields(localFine,&numFields);
3532:   {
3533:     PetscInt maxFields = PetscMax(1,numFields) + 1;
3534:     PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
3535:   }

3537:   { /* winnow fine points that don't have global dofs out of the sf */
3538:     PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;

3540:     for (p = pStartF, numPointsWithDofs = 0; p < pEndF; p++) {
3541:       PetscSectionGetDof(globalFine,p,&dof);
3542:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3543:       if ((dof - cdof) > 0) {
3544:         numPointsWithDofs++;

3546:         PetscSectionGetDof(localFine,p,&dof);
3547:         PetscSectionSetDof(leafIndicesSec,p,dof + 1);
3548:       }
3549:     }
3550:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3551:     PetscSectionSetUp(leafIndicesSec);
3552:     PetscSectionGetStorageSize(leafIndicesSec,&numIndices);
3553:     PetscMalloc1(numIndices,&leafIndices);
3554:     for (p = pStartF, offset = 0; p < pEndF; p++) {
3555:       PetscSectionGetDof(globalFine,p,&dof);
3556:       PetscSectionGetConstraintDof(globalFine,p,&cdof);
3557:       if ((dof - cdof) > 0) {
3558:         PetscInt off, gOff;
3559:         PetscInt *pInd;

3561:         pointsWithDofs[offset++] = p - pStartF;

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

3565:         pInd = &leafIndices[off + 1];
3566:         PetscSectionGetOffset(globalFine,p,&gOff);

3568:         offsets[0] = 0;
3569:         if (numFields) {
3570:           PetscInt f;

3572:           for (f = 0; f < numFields; f++) {
3573:             PetscInt fDof;
3574:             PetscSectionGetFieldDof(localFine,p,f,&fDof);
3575:             offsets[f + 1] = fDof + offsets[f];
3576:           }
3577:           indicesPointFields_private(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);
3578:         }
3579:         else {
3580:           indicesPoint_private(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);
3581:         }
3582:       }
3583:     }
3584:     PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
3585:     PetscFree(pointsWithDofs);
3586:   }

3588:   DMPlexGetChart(coarse,&pStartC,&pEndC);
3589:   DMGetDefaultSection(coarse,&localCoarse);
3590:   DMGetDefaultGlobalSection(coarse,&globalCoarse);
3591:   PetscSectionGetMaxDof(localCoarse,&maxDof);

3593:   { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3594:     MPI_Datatype threeInt;
3595:     PetscMPIInt  rank;
3596:     PetscInt     (*parentNodeAndIdCoarse)[3];
3597:     PetscInt     (*parentNodeAndIdFine)[3];
3598:     PetscInt     p, nleaves, nleavesToParents;
3599:     PetscSF      pointSF, sfToParents;
3600:     const PetscInt *ilocal;
3601:     const PetscSFNode *iremote;
3602:     PetscSFNode  *iremoteToParents;
3603:     PetscInt     *ilocalToParents;

3605:     MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
3606:     MPI_Type_contiguous(3,MPIU_INT,&threeInt);
3607:     MPI_Type_commit(&threeInt);
3608:     PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);
3609:     DMGetPointSF(coarse,&pointSF);
3610:     PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);
3611:     for (p = pStartC; p < pEndC; p++) {
3612:       PetscInt parent, childId;
3613:       DMPlexGetTreeParent(coarse,p,&parent,&childId);
3614:       parentNodeAndIdCoarse[p - pStartC][0] = rank;
3615:       parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3616:       parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3617:       if (nleaves > 0) {
3618:         PetscInt leaf = -1;

3620:         if (ilocal) {
3621:           PetscFindInt(parent,nleaves,ilocal,&leaf);
3622:         }
3623:         else {
3624:           leaf = p - pStartC;
3625:         }
3626:         if (leaf >= 0) {
3627:           parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3628:           parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3629:         }
3630:       }
3631:     }
3632:     for (p = pStartF; p < pEndF; p++) {
3633:       parentNodeAndIdFine[p - pStartF][0] = -1;
3634:       parentNodeAndIdFine[p - pStartF][1] = -1;
3635:       parentNodeAndIdFine[p - pStartF][2] = -1;
3636:     }
3637:     PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3638:     PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3639:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3640:       PetscInt dof;

3642:       PetscSectionGetDof(leafIndicesSec,p,&dof);
3643:       if (dof) {
3644:         PetscInt off;

3646:         PetscSectionGetOffset(leafIndicesSec,p,&off);
3647:         leafIndices[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3648:       }
3649:       if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3650:         nleavesToParents++;
3651:       }
3652:     }
3653:     PetscMalloc1(nleavesToParents,&ilocalToParents);
3654:     PetscMalloc1(nleavesToParents,&iremoteToParents);
3655:     for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3656:       if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3657:         ilocalToParents[nleavesToParents] = p - pStartF;
3658:         iremoteToParents[nleavesToParents].rank  = parentNodeAndIdFine[p-pStartF][0];
3659:         iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1];
3660:         nleavesToParents++;
3661:       }
3662:     }
3663:     PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);
3664:     PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);
3665:     PetscSFDestroy(&coarseToFineEmbedded);

3667:     coarseToFineEmbedded = sfToParents;

3669:     PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);
3670:     MPI_Type_free(&threeInt);
3671:   }

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

3677:     for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3678:       PetscSectionGetDof(globalCoarse,p,&dof);
3679:       PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3680:       if ((dof - cdof) > 0) {
3681:         numPointsWithDofs++;
3682:       }
3683:     }
3684:     PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3685:     for (p = pStartC, offset = 0; p < pEndC; p++) {
3686:       PetscSectionGetDof(globalCoarse,p,&dof);
3687:       PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3688:       if ((dof - cdof) > 0) {
3689:         pointsWithDofs[offset++] = p - pStartC;
3690:       }
3691:     }
3692:     PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);
3693:     PetscSFDestroy(&coarseToFineEmbedded);
3694:     PetscFree(pointsWithDofs);
3695:     coarseToFineEmbedded = sfDofsOnly;
3696:   }

3698:   /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3699:   PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);
3700:   PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);
3701:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);
3702:   PetscSectionSetChart(multiRootSec,pStartC,pEndC);
3703:   for (p = pStartC; p < pEndC; p++) {
3704:     PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);
3705:   }
3706:   PetscSectionSetUp(multiRootSec);
3707:   PetscSectionGetStorageSize(multiRootSec,&numMulti);
3708:   PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
3709:   { /* distribute the leaf section */
3710:     PetscSF multi, multiInv, indicesSF;
3711:     PetscInt *remoteOffsets, numRootIndices;

3713:     PetscSFGetMultiSF(coarseToFineEmbedded,&multi);
3714:     PetscSFCreateInverseSF(multi,&multiInv);
3715:     PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);
3716:     PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);
3717:     PetscFree(remoteOffsets);
3718:     PetscSFDestroy(&multiInv);
3719:     PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
3720:     PetscMalloc1(numRootIndices,&rootIndices);
3721:     PetscSFBcastBegin(indicesSF,MPIU_INT,leafIndices,rootIndices);
3722:     PetscSFBcastEnd(indicesSF,MPIU_INT,leafIndices,rootIndices);
3723:     PetscSFDestroy(&indicesSF);
3724:   }
3725:   PetscSectionDestroy(&leafIndicesSec);
3726:   PetscFree(leafIndices);
3727:   PetscSFDestroy(&coarseToFineEmbedded);

3729:   PetscMalloc1(maxDof,&parentIndices);

3731:   /* count indices */
3732:   MatGetLayouts(mat,&rowMap,&colMap);
3733:   PetscLayoutSetUp(rowMap);
3734:   PetscLayoutSetUp(colMap);
3735:   PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
3736:   PetscLayoutGetRange(colMap,&colStart,&colEnd);
3737:   PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);
3738:   for (p = pStartC; p < pEndC; p++) {
3739:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3741:     PetscSectionGetDof(globalCoarse,p,&dof);
3742:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3743:     if ((dof - cdof) <= 0) continue;
3744:     PetscSectionGetOffset(globalCoarse,p,&gOff);

3746:     rowOffsets[0] = 0;
3747:     offsetsCopy[0] = 0;
3748:     if (numFields) {
3749:       PetscInt f;

3751:       for (f = 0; f < numFields; f++) {
3752:         PetscInt fDof;
3753:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3754:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3755:       }
3756:       indicesPointFields_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,0,parentIndices);
3757:     }
3758:     else {
3759:       indicesPoint_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,0,parentIndices);
3760:       rowOffsets[1] = offsetsCopy[0];
3761:     }

3763:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
3764:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
3765:     leafEnd = leafStart + numLeaves;
3766:     for (l = leafStart; l < leafEnd; l++) {
3767:       PetscInt numIndices, childId, offset;
3768:       const PetscInt *childIndices;

3770:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3771:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
3772:       childId = rootIndices[offset++];
3773:       childIndices = &rootIndices[offset];
3774:       numIndices--;

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

3779:         for (i = 0; i < numIndices; i++) {
3780:           PetscInt colIndex = childIndices[i];
3781:           PetscInt rowIndex = parentIndices[i];
3782:           if (rowIndex < 0) continue;
3783:           if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse");
3784:           if (colIndex >= colStart && colIndex < colEnd) {
3785:             nnzD[rowIndex - rowStart] = 1;
3786:           }
3787:           else {
3788:             nnzO[rowIndex - rowStart] = 1;
3789:           }
3790:         }
3791:       }
3792:       else {
3793:         PetscInt parentId, f, lim;

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

3797:         lim = PetscMax(1,numFields);
3798:         offsets[0] = 0;
3799:         if (numFields) {
3800:           PetscInt f;

3802:           for (f = 0; f < numFields; f++) {
3803:             PetscInt fDof;
3804:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

3806:             offsets[f + 1] = fDof + offsets[f];
3807:           }
3808:         }
3809:         else {
3810:           PetscInt cDof;

3812:           PetscSectionGetDof(cSecRef,childId,&cDof);
3813:           offsets[1] = cDof;
3814:         }
3815:         for (f = 0; f < lim; f++) {
3816:           PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3817:           PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3818:           PetscInt i, numD = 0, numO = 0;

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

3823:             if (colIndex < 0) continue;
3824:             if (colIndex >= colStart && colIndex < colEnd) {
3825:               numD++;
3826:             }
3827:             else {
3828:               numO++;
3829:             }
3830:           }
3831:           for (i = parentStart; i < parentEnd; i++) {
3832:             PetscInt rowIndex = parentIndices[i];

3834:             if (rowIndex < 0) continue;
3835:             nnzD[rowIndex - rowStart] += numD;
3836:             nnzO[rowIndex - rowStart] += numO;
3837:           }
3838:         }
3839:       }
3840:     }
3841:   }
3842:   /* preallocate */
3843:   MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);
3844:   PetscFree2(nnzD,nnzO);
3845:   /* insert values */
3846:   DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3847:   for (p = pStartC; p < pEndC; p++) {
3848:     PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;

3850:     PetscSectionGetDof(globalCoarse,p,&dof);
3851:     PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3852:     if ((dof - cdof) <= 0) continue;
3853:     PetscSectionGetOffset(globalCoarse,p,&gOff);

3855:     rowOffsets[0] = 0;
3856:     offsetsCopy[0] = 0;
3857:     if (numFields) {
3858:       PetscInt f;

3860:       for (f = 0; f < numFields; f++) {
3861:         PetscInt fDof;
3862:         PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3863:         rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3864:       }
3865:       indicesPointFields_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,0,parentIndices);
3866:     }
3867:     else {
3868:       indicesPoint_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,0,parentIndices);
3869:       rowOffsets[1] = offsetsCopy[0];
3870:     }

3872:     PetscSectionGetDof(multiRootSec,p,&numLeaves);
3873:     PetscSectionGetOffset(multiRootSec,p,&leafStart);
3874:     leafEnd = leafStart + numLeaves;
3875:     for (l = leafStart; l < leafEnd; l++) {
3876:       PetscInt numIndices, childId, offset;
3877:       const PetscInt *childIndices;

3879:       PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3880:       PetscSectionGetOffset(rootIndicesSec,l,&offset);
3881:       childId = rootIndices[offset++];
3882:       childIndices = &rootIndices[offset];
3883:       numIndices--;

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

3888:         for (i = 0; i < numIndices; i++) {
3889:           MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);
3890:         }
3891:       }
3892:       else {
3893:         PetscInt parentId, f, lim;

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

3897:         lim = PetscMax(1,numFields);
3898:         offsets[0] = 0;
3899:         if (numFields) {
3900:           PetscInt f;

3902:           for (f = 0; f < numFields; f++) {
3903:             PetscInt fDof;
3904:             PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);

3906:             offsets[f + 1] = fDof + offsets[f];
3907:           }
3908:         }
3909:         else {
3910:           PetscInt cDof;

3912:           PetscSectionGetDof(cSecRef,childId,&cDof);
3913:           offsets[1] = cDof;
3914:         }
3915:         for (f = 0; f < lim; f++) {
3916:           PetscScalar    *childMat   = &childrenMats[childId - pRefStart][f][0];
3917:           PetscInt       *rowIndices = &parentIndices[rowOffsets[f]];
3918:           const PetscInt *colIndices = &childIndices[offsets[f]];

3920:           MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);
3921:         }
3922:       }
3923:     }
3924:   }
3925:   PetscSectionDestroy(&multiRootSec);
3926:   PetscSectionDestroy(&rootIndicesSec);
3927:   PetscFree(parentIndices);
3928:   DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3929:   PetscFree(rootIndices);
3930:   PetscFree3(offsets,offsetsCopy,rowOffsets);

3932:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3933:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3934:   return(0);
3935: }