Actual source code: plex.c

petsc-3.7.6 2017-04-24
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <petsc/private/isimpl.h>
  3: #include <petsc/private/vecimpl.h>
  4: #include <petscsf.h>
  5: #include <petscds.h>

  7: /* Logging support */
  8: PetscLogEvent DMPLEX_Interpolate, PETSCPARTITIONER_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeOverlap, DMPLEX_DistributeField, DMPLEX_DistributeData, DMPLEX_Migrate, DMPLEX_GlobalToNaturalBegin, DMPLEX_GlobalToNaturalEnd, DMPLEX_NaturalToGlobalBegin, DMPLEX_NaturalToGlobalEnd, DMPLEX_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM, DMPLEX_InterpolatorFEM, DMPLEX_InjectorFEM, DMPLEX_IntegralFEM, DMPLEX_CreateGmsh;

 10: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
 11: PETSC_EXTERN PetscErrorCode VecLoad_Default(Vec, PetscViewer);

 15: PetscErrorCode DMPlexGetFieldType_Internal(DM dm, PetscSection section, PetscInt field, PetscInt *sStart, PetscInt *sEnd, PetscViewerVTKFieldType *ft)
 16: {
 17:   PetscInt       dim, pStart, pEnd, vStart, vEnd, cStart, cEnd, cEndInterior, vdof = 0, cdof = 0;

 21:   *ft  = PETSC_VTK_POINT_FIELD;
 22:   DMGetDimension(dm, &dim);
 23:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 24:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 25:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
 26:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
 27:   PetscSectionGetChart(section, &pStart, &pEnd);
 28:   if (field >= 0) {
 29:     if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetFieldDof(section, vStart, field, &vdof);}
 30:     if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetFieldDof(section, cStart, field, &cdof);}
 31:   } else {
 32:     if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetDof(section, vStart, &vdof);}
 33:     if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetDof(section, cStart, &cdof);}
 34:   }
 35:   if (vdof) {
 36:     *sStart = vStart;
 37:     *sEnd   = vEnd;
 38:     if (vdof == dim) *ft = PETSC_VTK_POINT_VECTOR_FIELD;
 39:     else             *ft = PETSC_VTK_POINT_FIELD;
 40:   } else if (cdof) {
 41:     *sStart = cStart;
 42:     *sEnd   = cEnd;
 43:     if (cdof == dim) *ft = PETSC_VTK_CELL_VECTOR_FIELD;
 44:     else             *ft = PETSC_VTK_CELL_FIELD;
 45:   } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Could not classify input Vec for VTK");
 46:   return(0);
 47: }

 51: PetscErrorCode VecView_Plex_Local(Vec v, PetscViewer viewer)
 52: {
 53:   DM             dm;
 54:   PetscBool      isvtk, ishdf5, isseq;

 58:   VecGetDM(v, &dm);
 59:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
 60:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,  &isvtk);
 61:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
 62:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
 63:   if (isvtk || ishdf5) {
 64:     PetscInt  numFields;
 65:     PetscBool fem = PETSC_FALSE;

 67:     DMGetNumFields(dm, &numFields);
 68:     if (numFields) {
 69:       PetscObject fe;

 71:       DMGetField(dm, 0, &fe);
 72:       if (fe->classid == PETSCFE_CLASSID) fem = PETSC_TRUE;
 73:     }
 74:     if (fem) {DMPlexInsertBoundaryValues(dm, PETSC_TRUE, v, 0.0, NULL, NULL, NULL);}
 75:   }
 76:   if (isvtk) {
 77:     PetscSection            section;
 78:     PetscViewerVTKFieldType ft;
 79:     PetscInt                pStart, pEnd;

 81:     DMGetDefaultSection(dm, &section);
 82:     DMPlexGetFieldType_Internal(dm, section, PETSC_DETERMINE, &pStart, &pEnd, &ft);
 83:     PetscObjectReference((PetscObject) v);  /* viewer drops reference */
 84:     PetscViewerVTKAddField(viewer, (PetscObject) dm, DMPlexVTKWriteAll, ft, (PetscObject) v);
 85:   } else if (ishdf5) {
 86: #if defined(PETSC_HAVE_HDF5)
 87:     VecView_Plex_Local_HDF5(v, viewer);
 88: #else
 89:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
 90: #endif
 91:   } else {
 92:     if (isseq) {VecView_Seq(v, viewer);}
 93:     else       {VecView_MPI(v, viewer);}
 94:   }
 95:   return(0);
 96: }

100: PetscErrorCode VecView_Plex(Vec v, PetscViewer viewer)
101: {
102:   DM             dm;
103:   PetscBool      isvtk, ishdf5, isseq;

107:   VecGetDM(v, &dm);
108:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
109:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,  &isvtk);
110:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
111:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
112:   if (isvtk) {
113:     Vec         locv;
114:     const char *name;

116:     DMGetLocalVector(dm, &locv);
117:     PetscObjectGetName((PetscObject) v, &name);
118:     PetscObjectSetName((PetscObject) locv, name);
119:     DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
120:     DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
121:     VecView_Plex_Local(locv, viewer);
122:     DMRestoreLocalVector(dm, &locv);
123:   } else if (ishdf5) {
124: #if defined(PETSC_HAVE_HDF5)
125:     VecView_Plex_HDF5(v, viewer);
126: #else
127:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
128: #endif
129:   } else {
130:     if (isseq) {VecView_Seq(v, viewer);}
131:     else       {VecView_MPI(v, viewer);}
132:   }
133:   return(0);
134: }

138: PetscErrorCode VecView_Plex_Native(Vec originalv, PetscViewer viewer)
139: {
140:   DM                dm;
141:   MPI_Comm          comm;
142:   PetscViewerFormat format;
143:   Vec               v;
144:   PetscBool         isvtk, ishdf5;
145:   PetscErrorCode    ierr;

148:   VecGetDM(originalv, &dm);
149:   PetscObjectGetComm((PetscObject) originalv, &comm);
150:   if (!dm) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
151:   PetscViewerGetFormat(viewer, &format);
152:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
153:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,  &isvtk);
154:   if (format == PETSC_VIEWER_NATIVE) {
155:     const char *vecname;
156:     PetscInt    n, nroots;

158:     if (dm->sfNatural) {
159:       VecGetLocalSize(originalv, &n);
160:       PetscSFGetGraph(dm->sfNatural, &nroots, NULL, NULL, NULL);
161:       if (n == nroots) {
162:         DMGetGlobalVector(dm, &v);
163:         DMPlexGlobalToNaturalBegin(dm, originalv, v);
164:         DMPlexGlobalToNaturalEnd(dm, originalv, v);
165:         PetscObjectGetName((PetscObject) originalv, &vecname);
166:         PetscObjectSetName((PetscObject) v, vecname);
167:       } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "DM global to natural SF only handles global vectors");
168:     } else SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created");
169:   } else {
170:     /* we are viewing a natural DMPlex vec. */
171:     v = originalv;
172:   }
173:   if (ishdf5) {
174: #if defined(PETSC_HAVE_HDF5)
175:     VecView_Plex_HDF5_Native(v, viewer);
176: #else
177:     SETERRQ(comm, PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
178: #endif
179:   } else if (isvtk) {
180:     SETERRQ(comm, PETSC_ERR_SUP, "VTK format does not support viewing in natural order. Please switch to HDF5.");
181:   } else {
182:     PetscBool isseq;

184:     PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
185:     if (isseq) {VecView_Seq(v, viewer);}
186:     else       {VecView_MPI(v, viewer);}
187:   }
188:   if (format == PETSC_VIEWER_NATIVE) {DMRestoreGlobalVector(dm, &v);}
189:   return(0);
190: }

194: PetscErrorCode VecLoad_Plex_Local(Vec v, PetscViewer viewer)
195: {
196:   DM             dm;
197:   PetscBool      ishdf5;

201:   VecGetDM(v, &dm);
202:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
203:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
204:   if (ishdf5) {
205:     DM          dmBC;
206:     Vec         gv;
207:     const char *name;

209:     DMGetOutputDM(dm, &dmBC);
210:     DMGetGlobalVector(dmBC, &gv);
211:     PetscObjectGetName((PetscObject) v, &name);
212:     PetscObjectSetName((PetscObject) gv, name);
213:     VecLoad_Default(gv, viewer);
214:     DMGlobalToLocalBegin(dmBC, gv, INSERT_VALUES, v);
215:     DMGlobalToLocalEnd(dmBC, gv, INSERT_VALUES, v);
216:     DMRestoreGlobalVector(dmBC, &gv);
217:   } else {
218:     VecLoad_Default(v, viewer);
219:   }
220:   return(0);
221: }

225: PetscErrorCode VecLoad_Plex(Vec v, PetscViewer viewer)
226: {
227:   DM             dm;
228:   PetscBool      ishdf5;

232:   VecGetDM(v, &dm);
233:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
234:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
235:   if (ishdf5) {
236: #if defined(PETSC_HAVE_HDF5)
237:     VecLoad_Plex_HDF5(v, viewer);
238: #else
239:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
240: #endif
241:   } else {
242:     VecLoad_Default(v, viewer);
243:   }
244:   return(0);
245: }

249: PetscErrorCode VecLoad_Plex_Native(Vec originalv, PetscViewer viewer)
250: {
251:   DM                dm;
252:   PetscViewerFormat format;
253:   PetscBool         ishdf5;
254:   PetscErrorCode    ierr;

257:   VecGetDM(originalv, &dm);
258:   if (!dm) SETERRQ(PetscObjectComm((PetscObject) originalv), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
259:   PetscViewerGetFormat(viewer, &format);
260:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
261:   if (format == PETSC_VIEWER_NATIVE) {
262:     if (dm->sfNatural) {
263:       if (ishdf5) {
264: #if defined(PETSC_HAVE_HDF5)
265:         Vec         v;
266:         const char *vecname;

268:         DMGetGlobalVector(dm, &v);
269:         PetscObjectGetName((PetscObject) originalv, &vecname);
270:         PetscObjectSetName((PetscObject) v, vecname);
271:         VecLoad_Plex_HDF5_Native(v, viewer);
272:         DMPlexNaturalToGlobalBegin(dm, v, originalv);
273:         DMPlexNaturalToGlobalEnd(dm, v, originalv);
274:         DMRestoreGlobalVector(dm, &v);
275: #else
276:         SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
277: #endif
278:       } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Reading in natural order is not supported for anything but HDF5.");
279:     }
280:   }
281:   return(0);
282: }

286: PetscErrorCode DMPlexView_Ascii_Geometry(DM dm, PetscViewer viewer)
287: {
288:   PetscSection       coordSection;
289:   Vec                coordinates;
290:   DMLabel            depthLabel;
291:   const char        *name[4];
292:   const PetscScalar *a;
293:   PetscInt           dim, pStart, pEnd, cStart, cEnd, c;
294:   PetscErrorCode     ierr;

297:   DMGetDimension(dm, &dim);
298:   DMGetCoordinatesLocal(dm, &coordinates);
299:   DMGetCoordinateSection(dm, &coordSection);
300:   DMPlexGetDepthLabel(dm, &depthLabel);
301:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
302:   PetscSectionGetChart(coordSection, &pStart, &pEnd);
303:   VecGetArrayRead(coordinates, &a);
304:   name[0]     = "vertex";
305:   name[1]     = "edge";
306:   name[dim-1] = "face";
307:   name[dim]   = "cell";
308:   for (c = cStart; c < cEnd; ++c) {
309:     PetscInt *closure = NULL;
310:     PetscInt  closureSize, cl;

312:     PetscViewerASCIIPrintf(viewer, "Geometry for cell %D:\n", c);
313:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
314:     PetscViewerASCIIPushTab(viewer);
315:     for (cl = 0; cl < closureSize*2; cl += 2) {
316:       PetscInt point = closure[cl], depth, dof, off, d, p;

318:       if ((point < pStart) || (point >= pEnd)) continue;
319:       PetscSectionGetDof(coordSection, point, &dof);
320:       if (!dof) continue;
321:       DMLabelGetValue(depthLabel, point, &depth);
322:       PetscSectionGetOffset(coordSection, point, &off);
323:       PetscViewerASCIIPrintf(viewer, "%s %D coords:", name[depth], point);
324:       for (p = 0; p < dof/dim; ++p) {
325:         PetscViewerASCIIPrintf(viewer, " (");
326:         for (d = 0; d < dim; ++d) {
327:           if (d > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
328:           PetscViewerASCIIPrintf(viewer, "%g", PetscRealPart(a[off+p*dim+d]));
329:         }
330:         PetscViewerASCIIPrintf(viewer, ")");
331:       }
332:       PetscViewerASCIIPrintf(viewer, "\n");
333:     }
334:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
335:     PetscViewerASCIIPopTab(viewer);
336:   }
337:   VecRestoreArrayRead(coordinates, &a);
338:   return(0);
339: }

343: PetscErrorCode DMPlexView_Ascii(DM dm, PetscViewer viewer)
344: {
345:   DM_Plex          *mesh = (DM_Plex*) dm->data;
346:   DM                cdm;
347:   DMLabel           markers;
348:   PetscSection      coordSection;
349:   Vec               coordinates;
350:   PetscViewerFormat format;
351:   PetscErrorCode    ierr;

354:   DMGetCoordinateDM(dm, &cdm);
355:   DMGetDefaultSection(cdm, &coordSection);
356:   DMGetCoordinatesLocal(dm, &coordinates);
357:   PetscViewerGetFormat(viewer, &format);
358:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
359:     const char *name;
360:     PetscInt    maxConeSize, maxSupportSize;
361:     PetscInt    pStart, pEnd, p;
362:     PetscMPIInt rank, size;

364:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
365:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
366:     PetscObjectGetName((PetscObject) dm, &name);
367:     DMPlexGetChart(dm, &pStart, &pEnd);
368:     DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
369:     PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
370:     PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
371:     PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
372:     PetscViewerASCIIPushSynchronized(viewer);
373:     PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Max sizes cone: %D support: %D\n", rank,maxConeSize, maxSupportSize);
374:     for (p = pStart; p < pEnd; ++p) {
375:       PetscInt dof, off, s;

377:       PetscSectionGetDof(mesh->supportSection, p, &dof);
378:       PetscSectionGetOffset(mesh->supportSection, p, &off);
379:       for (s = off; s < off+dof; ++s) {
380:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D ----> %D\n", rank, p, mesh->supports[s]);
381:       }
382:     }
383:     PetscViewerFlush(viewer);
384:     PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
385:     for (p = pStart; p < pEnd; ++p) {
386:       PetscInt dof, off, c;

388:       PetscSectionGetDof(mesh->coneSection, p, &dof);
389:       PetscSectionGetOffset(mesh->coneSection, p, &off);
390:       for (c = off; c < off+dof; ++c) {
391:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D <---- %D (%D)\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]);
392:       }
393:     }
394:     PetscViewerFlush(viewer);
395:     PetscViewerASCIIPopSynchronized(viewer);
396:     PetscSectionGetChart(coordSection, &pStart, NULL);
397:     if (pStart >= 0) {PetscSectionVecView(coordSection, coordinates, viewer);}
398:     DMGetLabel(dm, "marker", &markers);
399:     DMLabelView(markers,viewer);
400:     if (size > 1) {
401:       PetscSF sf;

403:       DMGetPointSF(dm, &sf);
404:       PetscSFView(sf, viewer);
405:     }
406:     PetscViewerFlush(viewer);
407:   } else if (format == PETSC_VIEWER_ASCII_LATEX) {
408:     const char  *name, *color;
409:     const char  *defcolors[3]  = {"gray", "orange", "green"};
410:     const char  *deflcolors[4] = {"blue", "cyan", "red", "magenta"};
411:     PetscReal    scale         = 2.0;
412:     PetscBool    useNumbers    = PETSC_TRUE, useLabels, useColors;
413:     double       tcoords[3];
414:     PetscScalar *coords;
415:     PetscInt     numLabels, l, numColors, numLColors, dim, depth, cStart, cEnd, c, vStart, vEnd, v, eStart = 0, eEnd = 0, e, p;
416:     PetscMPIInt  rank, size;
417:     char         **names, **colors, **lcolors;

419:     DMGetDimension(dm, &dim);
420:     DMPlexGetDepth(dm, &depth);
421:     DMGetNumLabels(dm, &numLabels);
422:     numLabels  = PetscMax(numLabels, 10);
423:     numColors  = 10;
424:     numLColors = 10;
425:     PetscCalloc3(numLabels, &names, numColors, &colors, numLColors, &lcolors);
426:     PetscOptionsGetReal(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_scale", &scale, NULL);
427:     PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_numbers", &useNumbers, NULL);
428:     PetscOptionsGetStringArray(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_labels", names, &numLabels, &useLabels);
429:     if (!useLabels) numLabels = 0;
430:     PetscOptionsGetStringArray(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_colors", colors, &numColors, &useColors);
431:     if (!useColors) {
432:       numColors = 3;
433:       for (c = 0; c < numColors; ++c) {PetscStrallocpy(defcolors[c], &colors[c]);}
434:     }
435:     PetscOptionsGetStringArray(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_lcolors", lcolors, &numLColors, &useColors);
436:     if (!useColors) {
437:       numLColors = 4;
438:       for (c = 0; c < numLColors; ++c) {PetscStrallocpy(deflcolors[c], &lcolors[c]);}
439:     }
440:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
441:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
442:     PetscObjectGetName((PetscObject) dm, &name);
443:     PetscViewerASCIIPrintf(viewer, "\
444: \\documentclass[tikz]{standalone}\n\n\
445: \\usepackage{pgflibraryshapes}\n\
446: \\usetikzlibrary{backgrounds}\n\
447: \\usetikzlibrary{arrows}\n\
448: \\begin{document}\n");
449:     if (size > 1) {
450:       PetscViewerASCIIPrintf(viewer, "%s for process ", name);
451:       for (p = 0; p < size; ++p) {
452:         if (p > 0 && p == size-1) {
453:           PetscViewerASCIIPrintf(viewer, ", and ", colors[p%numColors], p);
454:         } else if (p > 0) {
455:           PetscViewerASCIIPrintf(viewer, ", ", colors[p%numColors], p);
456:         }
457:         PetscViewerASCIIPrintf(viewer, "{\\textcolor{%s}%D}", colors[p%numColors], p);
458:       }
459:       PetscViewerASCIIPrintf(viewer, ".\n\n\n");
460:     }
461:     PetscViewerASCIIPrintf(viewer, "\\begin{tikzpicture}[scale = %g,font=\\fontsize{8}{8}\\selectfont]\n", 1.0);
462:     /* Plot vertices */
463:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
464:     VecGetArray(coordinates, &coords);
465:     PetscViewerASCIIPushSynchronized(viewer);
466:     for (v = vStart; v < vEnd; ++v) {
467:       PetscInt  off, dof, d;
468:       PetscBool isLabeled = PETSC_FALSE;

470:       PetscSectionGetDof(coordSection, v, &dof);
471:       PetscSectionGetOffset(coordSection, v, &off);
472:       PetscViewerASCIISynchronizedPrintf(viewer, "\\path (");
473:       if (PetscUnlikely(dof > 3)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"coordSection vertex %D has dof %D > 3",v,dof);
474:       for (d = 0; d < dof; ++d) {
475:         tcoords[d] = (double) (scale*PetscRealPart(coords[off+d]));
476:         tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
477:       }
478:       /* Rotate coordinates since PGF makes z point out of the page instead of up */
479:       if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
480:       for (d = 0; d < dof; ++d) {
481:         if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
482:         PetscViewerASCIISynchronizedPrintf(viewer, "%g", tcoords[d]);
483:       }
484:       color = colors[rank%numColors];
485:       for (l = 0; l < numLabels; ++l) {
486:         PetscInt val;
487:         DMGetLabelValue(dm, names[l], v, &val);
488:         if (val >= 0) {color = lcolors[l%numLColors]; isLabeled = PETSC_TRUE; break;}
489:       }
490:       if (useNumbers) {
491:         PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [draw,shape=circle,color=%s] {%D};\n", v, rank, color, v);
492:       } else {
493:         PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [fill,inner sep=%dpt,shape=circle,color=%s] {};\n", v, rank, !isLabeled ? 1 : 2, color);
494:       }
495:     }
496:     VecRestoreArray(coordinates, &coords);
497:     PetscViewerFlush(viewer);
498:     /* Plot edges */
499:     if (depth > 1) {DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);}
500:     if (dim < 3 && useNumbers) {
501:       VecGetArray(coordinates, &coords);
502:       PetscViewerASCIIPrintf(viewer, "\\path\n");
503:       for (e = eStart; e < eEnd; ++e) {
504:         const PetscInt *cone;
505:         PetscInt        coneSize, offA, offB, dof, d;

507:         DMPlexGetConeSize(dm, e, &coneSize);
508:         if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %D cone should have two vertices, not %D", e, coneSize);
509:         DMPlexGetCone(dm, e, &cone);
510:         PetscSectionGetDof(coordSection, cone[0], &dof);
511:         PetscSectionGetOffset(coordSection, cone[0], &offA);
512:         PetscSectionGetOffset(coordSection, cone[1], &offB);
513:         PetscViewerASCIISynchronizedPrintf(viewer, "(");
514:         for (d = 0; d < dof; ++d) {
515:           tcoords[d] = (double) (scale*PetscRealPart(coords[offA+d]+coords[offB+d]));
516:           tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
517:         }
518:         /* Rotate coordinates since PGF makes z point out of the page instead of up */
519:         if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
520:         for (d = 0; d < dof; ++d) {
521:           if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
522:           PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)tcoords[d]);
523:         }
524:         color = colors[rank%numColors];
525:         for (l = 0; l < numLabels; ++l) {
526:           PetscInt val;
527:           DMGetLabelValue(dm, names[l], v, &val);
528:           if (val >= 0) {color = lcolors[l%numLColors]; break;}
529:         }
530:         PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [draw,shape=circle,color=%s] {%D} --\n", e, rank, color, e);
531:       }
532:       VecRestoreArray(coordinates, &coords);
533:       PetscViewerFlush(viewer);
534:       PetscViewerASCIIPrintf(viewer, "(0,0);\n");
535:     }
536:     /* Plot cells */
537:     if (dim == 3 || !useNumbers) {
538:       for (e = eStart; e < eEnd; ++e) {
539:         const PetscInt *cone;

541:         color = colors[rank%numColors];
542:         for (l = 0; l < numLabels; ++l) {
543:           PetscInt val;
544:           DMGetLabelValue(dm, names[l], e, &val);
545:           if (val >= 0) {color = lcolors[l%numLColors]; break;}
546:         }
547:         DMPlexGetCone(dm, e, &cone);
548:         PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] (%D_%d) -- (%D_%d);\n", color, cone[0], rank, cone[1], rank);
549:       }
550:     } else {
551:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
552:       for (c = cStart; c < cEnd; ++c) {
553:         PetscInt *closure = NULL;
554:         PetscInt  closureSize, firstPoint = -1;

556:         DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
557:         PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] ", colors[rank%numColors]);
558:         for (p = 0; p < closureSize*2; p += 2) {
559:           const PetscInt point = closure[p];

561:           if ((point < vStart) || (point >= vEnd)) continue;
562:           if (firstPoint >= 0) {PetscViewerASCIISynchronizedPrintf(viewer, " -- ");}
563:           PetscViewerASCIISynchronizedPrintf(viewer, "(%D_%d)", point, rank);
564:           if (firstPoint < 0) firstPoint = point;
565:         }
566:         /* Why doesn't this work? PetscViewerASCIISynchronizedPrintf(viewer, " -- cycle;\n"); */
567:         PetscViewerASCIISynchronizedPrintf(viewer, " -- (%D_%d);\n", firstPoint, rank);
568:         DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
569:       }
570:     }
571:     PetscViewerFlush(viewer);
572:     PetscViewerASCIIPopSynchronized(viewer);
573:     PetscViewerASCIIPrintf(viewer, "\\end{tikzpicture}\n");
574:     PetscViewerASCIIPrintf(viewer, "\\end{document}\n", name);
575:     for (l = 0; l < numLabels;  ++l) {PetscFree(names[l]);}
576:     for (c = 0; c < numColors;  ++c) {PetscFree(colors[c]);}
577:     for (c = 0; c < numLColors; ++c) {PetscFree(lcolors[c]);}
578:     PetscFree3(names, colors, lcolors);
579:   } else {
580:     MPI_Comm    comm;
581:     PetscInt   *sizes, *hybsizes;
582:     PetscInt    locDepth, depth, dim, d, pMax[4];
583:     PetscInt    pStart, pEnd, p;
584:     PetscInt    numLabels, l;
585:     const char *name;
586:     PetscMPIInt size;

588:     PetscObjectGetComm((PetscObject)dm,&comm);
589:     MPI_Comm_size(comm, &size);
590:     DMGetDimension(dm, &dim);
591:     PetscObjectGetName((PetscObject) dm, &name);
592:     if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dim);}
593:     else      {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dim);}
594:     DMPlexGetDepth(dm, &locDepth);
595:     MPIU_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
596:     DMPlexGetHybridBounds(dm, &pMax[depth], depth > 0 ? &pMax[depth-1] : NULL, &pMax[1], &pMax[0]);
597:     PetscMalloc2(size,&sizes,size,&hybsizes);
598:     if (depth == 1) {
599:       DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd);
600:       pEnd = pEnd - pStart;
601:       MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
602:       PetscViewerASCIIPrintf(viewer, "  %d-cells:", 0);
603:       for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
604:       PetscViewerASCIIPrintf(viewer, "\n");
605:       DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd);
606:       pEnd = pEnd - pStart;
607:       MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
608:       PetscViewerASCIIPrintf(viewer, "  %D-cells:", dim);
609:       for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
610:       PetscViewerASCIIPrintf(viewer, "\n");
611:     } else {
612:       for (d = 0; d <= dim; d++) {
613:         DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
614:         pEnd    -= pStart;
615:         pMax[d] -= pStart;
616:         MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
617:         MPI_Gather(&pMax[d], 1, MPIU_INT, hybsizes, 1, MPIU_INT, 0, comm);
618:         PetscViewerASCIIPrintf(viewer, "  %D-cells:", d);
619:         for (p = 0; p < size; ++p) {
620:           if (hybsizes[p] >= 0) {PetscViewerASCIIPrintf(viewer, " %D (%D)", sizes[p], sizes[p] - hybsizes[p]);}
621:           else                  {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
622:         }
623:         PetscViewerASCIIPrintf(viewer, "\n");
624:       }
625:     }
626:     PetscFree2(sizes,hybsizes);
627:     DMGetNumLabels(dm, &numLabels);
628:     if (numLabels) {PetscViewerASCIIPrintf(viewer, "Labels:\n");}
629:     for (l = 0; l < numLabels; ++l) {
630:       DMLabel         label;
631:       const char     *name;
632:       IS              valueIS;
633:       const PetscInt *values;
634:       PetscInt        numValues, v;

636:       DMGetLabelName(dm, l, &name);
637:       DMGetLabel(dm, name, &label);
638:       DMLabelGetNumValues(label, &numValues);
639:       PetscViewerASCIIPrintf(viewer, "  %s: %D strata of sizes (", name, numValues);
640:       DMLabelGetValueIS(label, &valueIS);
641:       ISGetIndices(valueIS, &values);
642:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
643:       for (v = 0; v < numValues; ++v) {
644:         PetscInt size;

646:         DMLabelGetStratumSize(label, values[v], &size);
647:         if (v > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
648:         PetscViewerASCIIPrintf(viewer, "%D", size);
649:       }
650:       PetscViewerASCIIPrintf(viewer, ")\n");
651:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
652:       ISRestoreIndices(valueIS, &values);
653:       ISDestroy(&valueIS);
654:     }
655:     DMGetCoarseDM(dm, &cdm);
656:     if (cdm) {
657:       PetscViewerASCIIPushTab(viewer);
658:       DMPlexView_Ascii(cdm, viewer);
659:       PetscViewerASCIIPopTab(viewer);
660:     }
661:   }
662:   return(0);
663: }

667: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
668: {
669:   PetscBool      iascii, ishdf5, isvtk;

675:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
676:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,   &isvtk);
677:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,  &ishdf5);
678:   if (iascii) {
679:     DMPlexView_Ascii(dm, viewer);
680:   } else if (ishdf5) {
681: #if defined(PETSC_HAVE_HDF5)
682:     PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
683:     DMPlexView_HDF5(dm, viewer);
684:     PetscViewerPopFormat(viewer);
685: #else
686:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
687: #endif
688:   }
689:   else if (isvtk) {
690:     DMPlexVTKWriteAll((PetscObject) dm,viewer);
691:   }
692:   return(0);
693: }

697: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
698: {
699:   PetscBool      isbinary, ishdf5;

705:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
706:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,   &ishdf5);
707:   if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");}
708:   else if (ishdf5) {
709: #if defined(PETSC_HAVE_HDF5)
710:     DMPlexLoad_HDF5(dm, viewer);
711: #else
712:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
713: #endif
714:   }
715:   return(0);
716: }


721: PetscErrorCode DMDestroy_Plex(DM dm)
722: {
723:   DM_Plex       *mesh = (DM_Plex*) dm->data;

727:   if (--mesh->refct > 0) return(0);
728:   PetscSectionDestroy(&mesh->coneSection);
729:   PetscFree(mesh->cones);
730:   PetscFree(mesh->coneOrientations);
731:   PetscSectionDestroy(&mesh->supportSection);
732:   PetscFree(mesh->supports);
733:   PetscFree(mesh->facesTmp);
734:   PetscFree(mesh->tetgenOpts);
735:   PetscFree(mesh->triangleOpts);
736:   PetscPartitionerDestroy(&mesh->partitioner);
737:   DMLabelDestroy(&mesh->subpointMap);
738:   ISDestroy(&mesh->globalVertexNumbers);
739:   ISDestroy(&mesh->globalCellNumbers);
740:   PetscSectionDestroy(&mesh->anchorSection);
741:   ISDestroy(&mesh->anchorIS);
742:   PetscSectionDestroy(&mesh->parentSection);
743:   PetscFree(mesh->parents);
744:   PetscFree(mesh->childIDs);
745:   PetscSectionDestroy(&mesh->childSection);
746:   PetscFree(mesh->children);
747:   DMDestroy(&mesh->referenceTree);
748:   PetscGridHashDestroy(&mesh->lbox);
749:   /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
750:   PetscFree(mesh);
751:   return(0);
752: }

756: PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
757: {
758:   PetscSection   sectionGlobal;
759:   PetscInt       bs = -1;
760:   PetscInt       localSize;
761:   PetscBool      isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock;
763:   MatType        mtype;
764:   ISLocalToGlobalMapping ltog;

767:   MatInitializePackage();
768:   mtype = dm->mattype;
769:   DMGetDefaultGlobalSection(dm, &sectionGlobal);
770:   /* PetscSectionGetStorageSize(sectionGlobal, &localSize); */
771:   PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
772:   MatCreate(PetscObjectComm((PetscObject)dm), J);
773:   MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
774:   MatSetType(*J, mtype);
775:   MatSetFromOptions(*J);
776:   PetscStrcmp(mtype, MATSHELL, &isShell);
777:   PetscStrcmp(mtype, MATBAIJ, &isBlock);
778:   PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
779:   PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
780:   PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
781:   PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
782:   PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
783:   if (!isShell) {
784:     PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
785:     PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin;

787:     if (bs < 0) {
788:       if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
789:         PetscInt pStart, pEnd, p, dof, cdof;

791:         PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
792:         for (p = pStart; p < pEnd; ++p) {
793:           PetscInt bdof;

795:           PetscSectionGetDof(sectionGlobal, p, &dof);
796:           PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
797:           bdof = PetscAbs(dof) - cdof;
798:           if (bdof) {
799:             if (bs < 0) {
800:               bs = bdof;
801:             } else if (bs != bdof) {
802:               /* Layout does not admit a pointwise block size */
803:               bs = 1;
804:               break;
805:             }
806:           }
807:         }
808:         /* Must have same blocksize on all procs (some might have no points) */
809:         bsLocal = bs;
810:         MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
811:         bsLocal = bs < 0 ? bsMax : bs;
812:         MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
813:         if (bsMin != bsMax) {
814:           bs = 1;
815:         } else {
816:           bs = bsMax;
817:         }
818:       } else {
819:         bs = 1;
820:       }
821:     }
822:     PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
823:     DMPlexPreallocateOperator(dm, bs, dnz, onz, dnzu, onzu, *J, fillMatrix);
824:     PetscFree4(dnz, onz, dnzu, onzu);

826:     /* Set localtoglobalmapping on the matrix for MatSetValuesLocal() to work */
827:     DMGetLocalToGlobalMapping(dm,&ltog);
828:     MatSetLocalToGlobalMapping(*J,ltog,ltog);
829:   }
830:   return(0);
831: }

835: /*@
836:   DMPlexGetChart - Return the interval for all mesh points [pStart, pEnd)

838:   Not collective

840:   Input Parameter:
841: . mesh - The DMPlex

843:   Output Parameters:
844: + pStart - The first mesh point
845: - pEnd   - The upper bound for mesh points

847:   Level: beginner

849: .seealso: DMPlexCreate(), DMPlexSetChart()
850: @*/
851: PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
852: {
853:   DM_Plex       *mesh = (DM_Plex*) dm->data;

858:   PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
859:   return(0);
860: }

864: /*@
865:   DMPlexSetChart - Set the interval for all mesh points [pStart, pEnd)

867:   Not collective

869:   Input Parameters:
870: + mesh - The DMPlex
871: . pStart - The first mesh point
872: - pEnd   - The upper bound for mesh points

874:   Output Parameters:

876:   Level: beginner

878: .seealso: DMPlexCreate(), DMPlexGetChart()
879: @*/
880: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
881: {
882:   DM_Plex       *mesh = (DM_Plex*) dm->data;

887:   PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
888:   PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
889:   return(0);
890: }

894: /*@
895:   DMPlexGetConeSize - Return the number of in-edges for this point in the Sieve DAG

897:   Not collective

899:   Input Parameters:
900: + mesh - The DMPlex
901: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

903:   Output Parameter:
904: . size - The cone size for point p

906:   Level: beginner

908: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
909: @*/
910: PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size)
911: {
912:   DM_Plex       *mesh = (DM_Plex*) dm->data;

918:   PetscSectionGetDof(mesh->coneSection, p, size);
919:   return(0);
920: }

924: /*@
925:   DMPlexSetConeSize - Set the number of in-edges for this point in the Sieve DAG

927:   Not collective

929:   Input Parameters:
930: + mesh - The DMPlex
931: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
932: - size - The cone size for point p

934:   Output Parameter:

936:   Note:
937:   This should be called after DMPlexSetChart().

939:   Level: beginner

941: .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart()
942: @*/
943: PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size)
944: {
945:   DM_Plex       *mesh = (DM_Plex*) dm->data;

950:   PetscSectionSetDof(mesh->coneSection, p, size);

952:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
953:   return(0);
954: }

958: /*@
959:   DMPlexAddConeSize - Add the given number of in-edges to this point in the Sieve DAG

961:   Not collective

963:   Input Parameters:
964: + mesh - The DMPlex
965: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
966: - size - The additional cone size for point p

968:   Output Parameter:

970:   Note:
971:   This should be called after DMPlexSetChart().

973:   Level: beginner

975: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart()
976: @*/
977: PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size)
978: {
979:   DM_Plex       *mesh = (DM_Plex*) dm->data;
980:   PetscInt       csize;

985:   PetscSectionAddDof(mesh->coneSection, p, size);
986:   PetscSectionGetDof(mesh->coneSection, p, &csize);

988:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
989:   return(0);
990: }

994: /*@C
995:   DMPlexGetCone - Return the points on the in-edges for this point in the Sieve DAG

997:   Not collective

999:   Input Parameters:
1000: + mesh - The DMPlex
1001: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1003:   Output Parameter:
1004: . cone - An array of points which are on the in-edges for point p

1006:   Level: beginner

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

1012:   You must also call DMPlexRestoreCone() after you finish using the returned array.

1014: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart()
1015: @*/
1016: PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
1017: {
1018:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1019:   PetscInt       off;

1025:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1026:   *cone = &mesh->cones[off];
1027:   return(0);
1028: }

1032: /*@
1033:   DMPlexSetCone - Set the points on the in-edges for this point in the Sieve DAG

1035:   Not collective

1037:   Input Parameters:
1038: + mesh - The DMPlex
1039: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1040: - cone - An array of points which are on the in-edges for point p

1042:   Output Parameter:

1044:   Note:
1045:   This should be called after all calls to DMPlexSetConeSize() and DMSetUp().

1047:   Level: beginner

1049: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1050: @*/
1051: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
1052: {
1053:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1054:   PetscInt       pStart, pEnd;
1055:   PetscInt       dof, off, c;

1060:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1061:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1063:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1064:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1065:   for (c = 0; c < dof; ++c) {
1066:     if ((cone[c] < pStart) || (cone[c] >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D, %D)", cone[c], pStart, pEnd);
1067:     mesh->cones[off+c] = cone[c];
1068:   }
1069:   return(0);
1070: }

1074: /*@C
1075:   DMPlexGetConeOrientation - Return the orientations on the in-edges for this point in the Sieve DAG

1077:   Not collective

1079:   Input Parameters:
1080: + mesh - The DMPlex
1081: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1083:   Output Parameter:
1084: . coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1085:                     integer giving the prescription for cone traversal. If it is negative, the cone is
1086:                     traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1087:                     the index of the cone point on which to start.

1089:   Level: beginner

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

1095:   You must also call DMPlexRestoreConeOrientation() after you finish using the returned array.

1097: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
1098: @*/
1099: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
1100: {
1101:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1102:   PetscInt       off;

1107: #if defined(PETSC_USE_DEBUG)
1108:   {
1109:     PetscInt dof;
1110:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1112:   }
1113: #endif
1114:   PetscSectionGetOffset(mesh->coneSection, p, &off);

1116:   *coneOrientation = &mesh->coneOrientations[off];
1117:   return(0);
1118: }

1122: /*@
1123:   DMPlexSetConeOrientation - Set the orientations on the in-edges for this point in the Sieve DAG

1125:   Not collective

1127:   Input Parameters:
1128: + mesh - The DMPlex
1129: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1130: - coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1131:                     integer giving the prescription for cone traversal. If it is negative, the cone is
1132:                     traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1133:                     the index of the cone point on which to start.

1135:   Output Parameter:

1137:   Note:
1138:   This should be called after all calls to DMPlexSetConeSize() and DMSetUp().

1140:   Level: beginner

1142: .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1143: @*/
1144: PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
1145: {
1146:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1147:   PetscInt       pStart, pEnd;
1148:   PetscInt       dof, off, c;

1153:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1154:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1156:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1157:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1158:   for (c = 0; c < dof; ++c) {
1159:     PetscInt cdof, o = coneOrientation[c];

1161:     PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1162:     if (o && ((o < -(cdof+1)) || (o >= cdof))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone orientation %D is not in the valid range [%D. %D)", o, -(cdof+1), cdof);
1163:     mesh->coneOrientations[off+c] = o;
1164:   }
1165:   return(0);
1166: }

1170: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
1171: {
1172:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1173:   PetscInt       pStart, pEnd;
1174:   PetscInt       dof, off;

1179:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1180:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1181:   if ((conePoint < pStart) || (conePoint >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D, %D)", conePoint, pStart, pEnd);
1182:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1183:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1184:   if ((conePos < 0) || (conePos >= dof)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone position %D of point %D is not in the valid range [0, %D)", conePos, p, dof);
1185:   mesh->cones[off+conePos] = conePoint;
1186:   return(0);
1187: }

1191: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1192: {
1193:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1194:   PetscInt       pStart, pEnd;
1195:   PetscInt       dof, off;

1200:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1201:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1202:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1203:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1204:   if ((conePos < 0) || (conePos >= dof)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone position %D of point %D is not in the valid range [0, %D)", conePos, p, dof);
1205:   mesh->coneOrientations[off+conePos] = coneOrientation;
1206:   return(0);
1207: }

1211: /*@
1212:   DMPlexGetSupportSize - Return the number of out-edges for this point in the Sieve DAG

1214:   Not collective

1216:   Input Parameters:
1217: + mesh - The DMPlex
1218: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1220:   Output Parameter:
1221: . size - The support size for point p

1223:   Level: beginner

1225: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1226: @*/
1227: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1228: {
1229:   DM_Plex       *mesh = (DM_Plex*) dm->data;

1235:   PetscSectionGetDof(mesh->supportSection, p, size);
1236:   return(0);
1237: }

1241: /*@
1242:   DMPlexSetSupportSize - Set the number of out-edges for this point in the Sieve DAG

1244:   Not collective

1246:   Input Parameters:
1247: + mesh - The DMPlex
1248: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1249: - size - The support size for point p

1251:   Output Parameter:

1253:   Note:
1254:   This should be called after DMPlexSetChart().

1256:   Level: beginner

1258: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1259: @*/
1260: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1261: {
1262:   DM_Plex       *mesh = (DM_Plex*) dm->data;

1267:   PetscSectionSetDof(mesh->supportSection, p, size);

1269:   mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1270:   return(0);
1271: }

1275: /*@C
1276:   DMPlexGetSupport - Return the points on the out-edges for this point in the Sieve DAG

1278:   Not collective

1280:   Input Parameters:
1281: + mesh - The DMPlex
1282: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1284:   Output Parameter:
1285: . support - An array of points which are on the out-edges for point p

1287:   Level: beginner

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

1293:   You must also call DMPlexRestoreSupport() after you finish using the returned array.

1295: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1296: @*/
1297: PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1298: {
1299:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1300:   PetscInt       off;

1306:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1307:   *support = &mesh->supports[off];
1308:   return(0);
1309: }

1313: /*@
1314:   DMPlexSetSupport - Set the points on the out-edges for this point in the Sieve DAG

1316:   Not collective

1318:   Input Parameters:
1319: + mesh - The DMPlex
1320: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1321: - support - An array of points which are on the in-edges for point p

1323:   Output Parameter:

1325:   Note:
1326:   This should be called after all calls to DMPlexSetSupportSize() and DMSetUp().

1328:   Level: beginner

1330: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1331: @*/
1332: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1333: {
1334:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1335:   PetscInt       pStart, pEnd;
1336:   PetscInt       dof, off, c;

1341:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1342:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1344:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1345:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1346:   for (c = 0; c < dof; ++c) {
1347:     if ((support[c] < pStart) || (support[c] >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support point %D is not in the valid range [%D, %D)", support[c], pStart, pEnd);
1348:     mesh->supports[off+c] = support[c];
1349:   }
1350:   return(0);
1351: }

1355: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1356: {
1357:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1358:   PetscInt       pStart, pEnd;
1359:   PetscInt       dof, off;

1364:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1365:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1366:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1367:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1368:   if ((supportPoint < pStart) || (supportPoint >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support point %D is not in the valid range [%D, %D)", supportPoint, pStart, pEnd);
1369:   if (supportPos >= dof) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support position %D of point %D is not in the valid range [0, %D)", supportPos, p, dof);
1370:   mesh->supports[off+supportPos] = supportPoint;
1371:   return(0);
1372: }

1376: /*@C
1377:   DMPlexGetTransitiveClosure - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG

1379:   Not collective

1381:   Input Parameters:
1382: + mesh - The DMPlex
1383: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1384: . useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
1385: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used

1387:   Output Parameters:
1388: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1389: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]

1391:   Note:
1392:   If using internal storage (points is NULL on input), each call overwrites the last output.

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

1398:   The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.

1400:   Level: beginner

1402: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1403: @*/
1404: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1405: {
1406:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1407:   PetscInt       *closure, *fifo;
1408:   const PetscInt *tmp = NULL, *tmpO = NULL;
1409:   PetscInt        tmpSize, t;
1410:   PetscInt        depth       = 0, maxSize;
1411:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1412:   PetscErrorCode  ierr;

1416:   DMPlexGetDepth(dm, &depth);
1417:   /* This is only 1-level */
1418:   if (useCone) {
1419:     DMPlexGetConeSize(dm, p, &tmpSize);
1420:     DMPlexGetCone(dm, p, &tmp);
1421:     DMPlexGetConeOrientation(dm, p, &tmpO);
1422:   } else {
1423:     DMPlexGetSupportSize(dm, p, &tmpSize);
1424:     DMPlexGetSupport(dm, p, &tmp);
1425:   }
1426:   if (depth == 1) {
1427:     if (*points) {
1428:       closure = *points;
1429:     } else {
1430:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1431:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1432:     }
1433:     closure[0] = p; closure[1] = 0;
1434:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1435:       closure[closureSize]   = tmp[t];
1436:       closure[closureSize+1] = tmpO ? tmpO[t] : 0;
1437:     }
1438:     if (numPoints) *numPoints = closureSize/2;
1439:     if (points)    *points    = closure;
1440:     return(0);
1441:   }
1442:   {
1443:     PetscInt c, coneSeries, s,supportSeries;

1445:     c = mesh->maxConeSize;
1446:     coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1447:     s = mesh->maxSupportSize;
1448:     supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1449:     maxSize = 2*PetscMax(coneSeries,supportSeries);
1450:   }
1451:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1452:   if (*points) {
1453:     closure = *points;
1454:   } else {
1455:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1456:   }
1457:   closure[0] = p; closure[1] = 0;
1458:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1459:     const PetscInt cp = tmp[t];
1460:     const PetscInt co = tmpO ? tmpO[t] : 0;

1462:     closure[closureSize]   = cp;
1463:     closure[closureSize+1] = co;
1464:     fifo[fifoSize]         = cp;
1465:     fifo[fifoSize+1]       = co;
1466:   }
1467:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1468:   while (fifoSize - fifoStart) {
1469:     const PetscInt q   = fifo[fifoStart];
1470:     const PetscInt o   = fifo[fifoStart+1];
1471:     const PetscInt rev = o >= 0 ? 0 : 1;
1472:     const PetscInt off = rev ? -(o+1) : o;

1474:     if (useCone) {
1475:       DMPlexGetConeSize(dm, q, &tmpSize);
1476:       DMPlexGetCone(dm, q, &tmp);
1477:       DMPlexGetConeOrientation(dm, q, &tmpO);
1478:     } else {
1479:       DMPlexGetSupportSize(dm, q, &tmpSize);
1480:       DMPlexGetSupport(dm, q, &tmp);
1481:       tmpO = NULL;
1482:     }
1483:     for (t = 0; t < tmpSize; ++t) {
1484:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1485:       const PetscInt cp = tmp[i];
1486:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1487:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1488:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1489:       PetscInt       co = tmpO ? tmpO[i] : 0;
1490:       PetscInt       c;

1492:       if (rev) {
1493:         PetscInt childSize, coff;
1494:         DMPlexGetConeSize(dm, cp, &childSize);
1495:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1496:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1497:       }
1498:       /* Check for duplicate */
1499:       for (c = 0; c < closureSize; c += 2) {
1500:         if (closure[c] == cp) break;
1501:       }
1502:       if (c == closureSize) {
1503:         closure[closureSize]   = cp;
1504:         closure[closureSize+1] = co;
1505:         fifo[fifoSize]         = cp;
1506:         fifo[fifoSize+1]       = co;
1507:         closureSize           += 2;
1508:         fifoSize              += 2;
1509:       }
1510:     }
1511:     fifoStart += 2;
1512:   }
1513:   if (numPoints) *numPoints = closureSize/2;
1514:   if (points)    *points    = closure;
1515:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1516:   return(0);
1517: }

1521: /*@C
1522:   DMPlexGetTransitiveClosure_Internal - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG with a specified initial orientation

1524:   Not collective

1526:   Input Parameters:
1527: + mesh - The DMPlex
1528: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1529: . orientation - The orientation of the point
1530: . useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
1531: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used

1533:   Output Parameters:
1534: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1535: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]

1537:   Note:
1538:   If using internal storage (points is NULL on input), each call overwrites the last output.

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

1544:   The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.

1546:   Level: beginner

1548: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1549: @*/
1550: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1551: {
1552:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1553:   PetscInt       *closure, *fifo;
1554:   const PetscInt *tmp = NULL, *tmpO = NULL;
1555:   PetscInt        tmpSize, t;
1556:   PetscInt        depth       = 0, maxSize;
1557:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1558:   PetscErrorCode  ierr;

1562:   DMPlexGetDepth(dm, &depth);
1563:   /* This is only 1-level */
1564:   if (useCone) {
1565:     DMPlexGetConeSize(dm, p, &tmpSize);
1566:     DMPlexGetCone(dm, p, &tmp);
1567:     DMPlexGetConeOrientation(dm, p, &tmpO);
1568:   } else {
1569:     DMPlexGetSupportSize(dm, p, &tmpSize);
1570:     DMPlexGetSupport(dm, p, &tmp);
1571:   }
1572:   if (depth == 1) {
1573:     if (*points) {
1574:       closure = *points;
1575:     } else {
1576:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1577:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1578:     }
1579:     closure[0] = p; closure[1] = ornt;
1580:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1581:       const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1582:       closure[closureSize]   = tmp[i];
1583:       closure[closureSize+1] = tmpO ? tmpO[i] : 0;
1584:     }
1585:     if (numPoints) *numPoints = closureSize/2;
1586:     if (points)    *points    = closure;
1587:     return(0);
1588:   }
1589:   {
1590:     PetscInt c, coneSeries, s,supportSeries;

1592:     c = mesh->maxConeSize;
1593:     coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1594:     s = mesh->maxSupportSize;
1595:     supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1596:     maxSize = 2*PetscMax(coneSeries,supportSeries);
1597:   }
1598:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1599:   if (*points) {
1600:     closure = *points;
1601:   } else {
1602:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1603:   }
1604:   closure[0] = p; closure[1] = ornt;
1605:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1606:     const PetscInt i  = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1607:     const PetscInt cp = tmp[i];
1608:     PetscInt       co = tmpO ? tmpO[i] : 0;

1610:     if (ornt < 0) {
1611:       PetscInt childSize, coff;
1612:       DMPlexGetConeSize(dm, cp, &childSize);
1613:       coff = co < 0 ? -(tmpO[i]+1) : tmpO[i];
1614:       co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1615:     }
1616:     closure[closureSize]   = cp;
1617:     closure[closureSize+1] = co;
1618:     fifo[fifoSize]         = cp;
1619:     fifo[fifoSize+1]       = co;
1620:   }
1621:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1622:   while (fifoSize - fifoStart) {
1623:     const PetscInt q   = fifo[fifoStart];
1624:     const PetscInt o   = fifo[fifoStart+1];
1625:     const PetscInt rev = o >= 0 ? 0 : 1;
1626:     const PetscInt off = rev ? -(o+1) : o;

1628:     if (useCone) {
1629:       DMPlexGetConeSize(dm, q, &tmpSize);
1630:       DMPlexGetCone(dm, q, &tmp);
1631:       DMPlexGetConeOrientation(dm, q, &tmpO);
1632:     } else {
1633:       DMPlexGetSupportSize(dm, q, &tmpSize);
1634:       DMPlexGetSupport(dm, q, &tmp);
1635:       tmpO = NULL;
1636:     }
1637:     for (t = 0; t < tmpSize; ++t) {
1638:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1639:       const PetscInt cp = tmp[i];
1640:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1641:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1642:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1643:       PetscInt       co = tmpO ? tmpO[i] : 0;
1644:       PetscInt       c;

1646:       if (rev) {
1647:         PetscInt childSize, coff;
1648:         DMPlexGetConeSize(dm, cp, &childSize);
1649:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1650:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1651:       }
1652:       /* Check for duplicate */
1653:       for (c = 0; c < closureSize; c += 2) {
1654:         if (closure[c] == cp) break;
1655:       }
1656:       if (c == closureSize) {
1657:         closure[closureSize]   = cp;
1658:         closure[closureSize+1] = co;
1659:         fifo[fifoSize]         = cp;
1660:         fifo[fifoSize+1]       = co;
1661:         closureSize           += 2;
1662:         fifoSize              += 2;
1663:       }
1664:     }
1665:     fifoStart += 2;
1666:   }
1667:   if (numPoints) *numPoints = closureSize/2;
1668:   if (points)    *points    = closure;
1669:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1670:   return(0);
1671: }

1675: /*@C
1676:   DMPlexRestoreTransitiveClosure - Restore the array of points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG

1678:   Not collective

1680:   Input Parameters:
1681: + mesh - The DMPlex
1682: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1683: . useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
1684: . numPoints - The number of points in the closure, so points[] is of size 2*numPoints, zeroed on exit
1685: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...], zeroed on exit

1687:   Note:
1688:   If not using internal storage (points is not NULL on input), this call is unnecessary

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

1694:   The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.

1696:   Level: beginner

1698: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1699: @*/
1700: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1701: {

1708:   DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1709:   if (numPoints) *numPoints = 0;
1710:   return(0);
1711: }

1715: /*@
1716:   DMPlexGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG

1718:   Not collective

1720:   Input Parameter:
1721: . mesh - The DMPlex

1723:   Output Parameters:
1724: + maxConeSize - The maximum number of in-edges
1725: - maxSupportSize - The maximum number of out-edges

1727:   Level: beginner

1729: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1730: @*/
1731: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1732: {
1733:   DM_Plex *mesh = (DM_Plex*) dm->data;

1737:   if (maxConeSize)    *maxConeSize    = mesh->maxConeSize;
1738:   if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1739:   return(0);
1740: }

1744: PetscErrorCode DMSetUp_Plex(DM dm)
1745: {
1746:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1747:   PetscInt       size;

1752:   PetscSectionSetUp(mesh->coneSection);
1753:   PetscSectionGetStorageSize(mesh->coneSection, &size);
1754:   PetscMalloc1(size, &mesh->cones);
1755:   PetscCalloc1(size, &mesh->coneOrientations);
1756:   if (mesh->maxSupportSize) {
1757:     PetscSectionSetUp(mesh->supportSection);
1758:     PetscSectionGetStorageSize(mesh->supportSection, &size);
1759:     PetscMalloc1(size, &mesh->supports);
1760:   }
1761:   return(0);
1762: }

1766: PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1767: {

1771:   if (subdm) {DMClone(dm, subdm);}
1772:   DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1773:   return(0);
1774: }

1778: /*@
1779:   DMPlexSymmetrize - Creates support (out-edge) information from cone (in-edge) inoformation

1781:   Not collective

1783:   Input Parameter:
1784: . mesh - The DMPlex

1786:   Output Parameter:

1788:   Note:
1789:   This should be called after all calls to DMPlexSetCone()

1791:   Level: beginner

1793: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1794: @*/
1795: PetscErrorCode DMPlexSymmetrize(DM dm)
1796: {
1797:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1798:   PetscInt      *offsets;
1799:   PetscInt       supportSize;
1800:   PetscInt       pStart, pEnd, p;

1805:   if (mesh->supports) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Supports were already setup in this DMPlex");
1806:   /* Calculate support sizes */
1807:   DMPlexGetChart(dm, &pStart, &pEnd);
1808:   for (p = pStart; p < pEnd; ++p) {
1809:     PetscInt dof, off, c;

1811:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1812:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1813:     for (c = off; c < off+dof; ++c) {
1814:       PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1815:     }
1816:   }
1817:   for (p = pStart; p < pEnd; ++p) {
1818:     PetscInt dof;

1820:     PetscSectionGetDof(mesh->supportSection, p, &dof);

1822:     mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1823:   }
1824:   PetscSectionSetUp(mesh->supportSection);
1825:   /* Calculate supports */
1826:   PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1827:   PetscMalloc1(supportSize, &mesh->supports);
1828:   PetscCalloc1(pEnd - pStart, &offsets);
1829:   for (p = pStart; p < pEnd; ++p) {
1830:     PetscInt dof, off, c;

1832:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1833:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1834:     for (c = off; c < off+dof; ++c) {
1835:       const PetscInt q = mesh->cones[c];
1836:       PetscInt       offS;

1838:       PetscSectionGetOffset(mesh->supportSection, q, &offS);

1840:       mesh->supports[offS+offsets[q]] = p;
1841:       ++offsets[q];
1842:     }
1843:   }
1844:   PetscFree(offsets);
1845:   return(0);
1846: }

1850: /*@
1851:   DMPlexStratify - The Sieve DAG for most topologies is a graded poset (http://en.wikipedia.org/wiki/Graded_poset), and
1852:   can be illustrated by Hasse Diagram (a http://en.wikipedia.org/wiki/Hasse_diagram). The strata group all points of the
1853:   same grade, and this function calculates the strata. This grade can be seen as the height (or depth) of the point in
1854:   the DAG.

1856:   Collective on dm

1858:   Input Parameter:
1859: . mesh - The DMPlex

1861:   Output Parameter:

1863:   Notes:
1864:   Concretely, DMPlexStratify() creates a new label named "depth" containing the dimension of each element: 0 for vertices,
1865:   1 for edges, and so on.  The depth label can be accessed through DMPlexGetDepthLabel() or DMPlexGetDepthStratum(), or
1866:   manually via DMGetLabel().  The height is defined implicitly by height = maxDimension - depth, and can be accessed
1867:   via DMPlexGetHeightStratum().  For example, cells have height 0 and faces have height 1.

1869:   DMPlexStratify() should be called after all calls to DMPlexSymmetrize()

1871:   Level: beginner

1873: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1874: @*/
1875: PetscErrorCode DMPlexStratify(DM dm)
1876: {
1877:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1878:   DMLabel        label;
1879:   PetscInt       pStart, pEnd, p;
1880:   PetscInt       numRoots = 0, numLeaves = 0;

1885:   PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1886:   /* Calculate depth */
1887:   DMPlexGetChart(dm, &pStart, &pEnd);
1888:   DMCreateLabel(dm, "depth");
1889:   DMPlexGetDepthLabel(dm, &label);
1890:   /* Initialize roots and count leaves */
1891:   for (p = pStart; p < pEnd; ++p) {
1892:     PetscInt coneSize, supportSize;

1894:     DMPlexGetConeSize(dm, p, &coneSize);
1895:     DMPlexGetSupportSize(dm, p, &supportSize);
1896:     if (!coneSize && supportSize) {
1897:       ++numRoots;
1898:       DMLabelSetValue(label, p, 0);
1899:     } else if (!supportSize && coneSize) {
1900:       ++numLeaves;
1901:     } else if (!supportSize && !coneSize) {
1902:       /* Isolated points */
1903:       DMLabelSetValue(label, p, 0);
1904:     }
1905:   }
1906:   if (numRoots + numLeaves == (pEnd - pStart)) {
1907:     for (p = pStart; p < pEnd; ++p) {
1908:       PetscInt coneSize, supportSize;

1910:       DMPlexGetConeSize(dm, p, &coneSize);
1911:       DMPlexGetSupportSize(dm, p, &supportSize);
1912:       if (!supportSize && coneSize) {
1913:         DMLabelSetValue(label, p, 1);
1914:       }
1915:     }
1916:   } else {
1917:     IS       pointIS;
1918:     PetscInt numPoints = 0, level = 0;

1920:     DMLabelGetStratumIS(label, level, &pointIS);
1921:     if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1922:     while (numPoints) {
1923:       const PetscInt *points;
1924:       const PetscInt  newLevel = level+1;

1926:       ISGetIndices(pointIS, &points);
1927:       for (p = 0; p < numPoints; ++p) {
1928:         const PetscInt  point = points[p];
1929:         const PetscInt *support;
1930:         PetscInt        supportSize, s;

1932:         DMPlexGetSupportSize(dm, point, &supportSize);
1933:         DMPlexGetSupport(dm, point, &support);
1934:         for (s = 0; s < supportSize; ++s) {
1935:           DMLabelSetValue(label, support[s], newLevel);
1936:         }
1937:       }
1938:       ISRestoreIndices(pointIS, &points);
1939:       ++level;
1940:       ISDestroy(&pointIS);
1941:       DMLabelGetStratumIS(label, level, &pointIS);
1942:       if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1943:       else         {numPoints = 0;}
1944:     }
1945:     ISDestroy(&pointIS);
1946:   }
1947:   { /* just in case there is an empty process */
1948:     PetscInt numValues, maxValues = 0, v;

1950:     DMLabelGetNumValues(label,&numValues);
1951:     MPI_Allreduce(&numValues,&maxValues,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
1952:     for (v = numValues; v < maxValues; v++) {
1953:       DMLabelAddStratum(label,v);
1954:     }
1955:   }

1957:   DMLabelGetState(label, &mesh->depthState);
1958:   PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1959:   return(0);
1960: }

1964: /*@C
1965:   DMPlexGetJoin - Get an array for the join of the set of points

1967:   Not Collective

1969:   Input Parameters:
1970: + dm - The DMPlex object
1971: . numPoints - The number of input points for the join
1972: - points - The input points

1974:   Output Parameters:
1975: + numCoveredPoints - The number of points in the join
1976: - coveredPoints - The points in the join

1978:   Level: intermediate

1980:   Note: Currently, this is restricted to a single level join

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

1986:   The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.

1988: .keywords: mesh
1989: .seealso: DMPlexRestoreJoin(), DMPlexGetMeet()
1990: @*/
1991: PetscErrorCode DMPlexGetJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1992: {
1993:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1994:   PetscInt      *join[2];
1995:   PetscInt       joinSize, i = 0;
1996:   PetscInt       dof, off, p, c, m;

2004:   DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[0]);
2005:   DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1]);
2006:   /* Copy in support of first point */
2007:   PetscSectionGetDof(mesh->supportSection, points[0], &dof);
2008:   PetscSectionGetOffset(mesh->supportSection, points[0], &off);
2009:   for (joinSize = 0; joinSize < dof; ++joinSize) {
2010:     join[i][joinSize] = mesh->supports[off+joinSize];
2011:   }
2012:   /* Check each successive support */
2013:   for (p = 1; p < numPoints; ++p) {
2014:     PetscInt newJoinSize = 0;

2016:     PetscSectionGetDof(mesh->supportSection, points[p], &dof);
2017:     PetscSectionGetOffset(mesh->supportSection, points[p], &off);
2018:     for (c = 0; c < dof; ++c) {
2019:       const PetscInt point = mesh->supports[off+c];

2021:       for (m = 0; m < joinSize; ++m) {
2022:         if (point == join[i][m]) {
2023:           join[1-i][newJoinSize++] = point;
2024:           break;
2025:         }
2026:       }
2027:     }
2028:     joinSize = newJoinSize;
2029:     i        = 1-i;
2030:   }
2031:   *numCoveredPoints = joinSize;
2032:   *coveredPoints    = join[i];
2033:   DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2034:   return(0);
2035: }

2039: /*@C
2040:   DMPlexRestoreJoin - Restore an array for the join of the set of points

2042:   Not Collective

2044:   Input Parameters:
2045: + dm - The DMPlex object
2046: . numPoints - The number of input points for the join
2047: - points - The input points

2049:   Output Parameters:
2050: + numCoveredPoints - The number of points in the join
2051: - coveredPoints - The points in the join

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

2057:   The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.

2059:   Level: intermediate

2061: .keywords: mesh
2062: .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet()
2063: @*/
2064: PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2065: {

2073:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2074:   if (numCoveredPoints) *numCoveredPoints = 0;
2075:   return(0);
2076: }

2080: /*@C
2081:   DMPlexGetFullJoin - Get an array for the join of the set of points

2083:   Not Collective

2085:   Input Parameters:
2086: + dm - The DMPlex object
2087: . numPoints - The number of input points for the join
2088: - points - The input points

2090:   Output Parameters:
2091: + numCoveredPoints - The number of points in the join
2092: - coveredPoints - The points in the join

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

2098:   The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.

2100:   Level: intermediate

2102: .keywords: mesh
2103: .seealso: DMPlexGetJoin(), DMPlexRestoreJoin(), DMPlexGetMeet()
2104: @*/
2105: PetscErrorCode DMPlexGetFullJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2106: {
2107:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2108:   PetscInt      *offsets, **closures;
2109:   PetscInt      *join[2];
2110:   PetscInt       depth = 0, maxSize, joinSize = 0, i = 0;
2111:   PetscInt       p, d, c, m, ms;


2120:   DMPlexGetDepth(dm, &depth);
2121:   PetscCalloc1(numPoints, &closures);
2122:   DMGetWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2123:   ms      = mesh->maxSupportSize;
2124:   maxSize = (ms > 1) ? ((PetscPowInt(ms,depth+1)-1)/(ms-1)) : depth + 1;
2125:   DMGetWorkArray(dm, maxSize, PETSC_INT, &join[0]);
2126:   DMGetWorkArray(dm, maxSize, PETSC_INT, &join[1]);

2128:   for (p = 0; p < numPoints; ++p) {
2129:     PetscInt closureSize;

2131:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closures[p]);

2133:     offsets[p*(depth+2)+0] = 0;
2134:     for (d = 0; d < depth+1; ++d) {
2135:       PetscInt pStart, pEnd, i;

2137:       DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
2138:       for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) {
2139:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2140:           offsets[p*(depth+2)+d+1] = i;
2141:           break;
2142:         }
2143:       }
2144:       if (i == closureSize) offsets[p*(depth+2)+d+1] = i;
2145:     }
2146:     if (offsets[p*(depth+2)+depth+1] != closureSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Total size of closure %D should be %D", offsets[p*(depth+2)+depth+1], closureSize);
2147:   }
2148:   for (d = 0; d < depth+1; ++d) {
2149:     PetscInt dof;

2151:     /* Copy in support of first point */
2152:     dof = offsets[d+1] - offsets[d];
2153:     for (joinSize = 0; joinSize < dof; ++joinSize) {
2154:       join[i][joinSize] = closures[0][(offsets[d]+joinSize)*2];
2155:     }
2156:     /* Check each successive cone */
2157:     for (p = 1; p < numPoints && joinSize; ++p) {
2158:       PetscInt newJoinSize = 0;

2160:       dof = offsets[p*(depth+2)+d+1] - offsets[p*(depth+2)+d];
2161:       for (c = 0; c < dof; ++c) {
2162:         const PetscInt point = closures[p][(offsets[p*(depth+2)+d]+c)*2];

2164:         for (m = 0; m < joinSize; ++m) {
2165:           if (point == join[i][m]) {
2166:             join[1-i][newJoinSize++] = point;
2167:             break;
2168:           }
2169:         }
2170:       }
2171:       joinSize = newJoinSize;
2172:       i        = 1-i;
2173:     }
2174:     if (joinSize) break;
2175:   }
2176:   *numCoveredPoints = joinSize;
2177:   *coveredPoints    = join[i];
2178:   for (p = 0; p < numPoints; ++p) {
2179:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, NULL, &closures[p]);
2180:   }
2181:   PetscFree(closures);
2182:   DMRestoreWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2183:   DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2184:   return(0);
2185: }

2189: /*@C
2190:   DMPlexGetMeet - Get an array for the meet of the set of points

2192:   Not Collective

2194:   Input Parameters:
2195: + dm - The DMPlex object
2196: . numPoints - The number of input points for the meet
2197: - points - The input points

2199:   Output Parameters:
2200: + numCoveredPoints - The number of points in the meet
2201: - coveredPoints - The points in the meet

2203:   Level: intermediate

2205:   Note: Currently, this is restricted to a single level meet

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

2211:   The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.

2213: .keywords: mesh
2214: .seealso: DMPlexRestoreMeet(), DMPlexGetJoin()
2215: @*/
2216: PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2217: {
2218:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2219:   PetscInt      *meet[2];
2220:   PetscInt       meetSize, i = 0;
2221:   PetscInt       dof, off, p, c, m;

2229:   DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[0]);
2230:   DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1]);
2231:   /* Copy in cone of first point */
2232:   PetscSectionGetDof(mesh->coneSection, points[0], &dof);
2233:   PetscSectionGetOffset(mesh->coneSection, points[0], &off);
2234:   for (meetSize = 0; meetSize < dof; ++meetSize) {
2235:     meet[i][meetSize] = mesh->cones[off+meetSize];
2236:   }
2237:   /* Check each successive cone */
2238:   for (p = 1; p < numPoints; ++p) {
2239:     PetscInt newMeetSize = 0;

2241:     PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2242:     PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2243:     for (c = 0; c < dof; ++c) {
2244:       const PetscInt point = mesh->cones[off+c];

2246:       for (m = 0; m < meetSize; ++m) {
2247:         if (point == meet[i][m]) {
2248:           meet[1-i][newMeetSize++] = point;
2249:           break;
2250:         }
2251:       }
2252:     }
2253:     meetSize = newMeetSize;
2254:     i        = 1-i;
2255:   }
2256:   *numCoveringPoints = meetSize;
2257:   *coveringPoints    = meet[i];
2258:   DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2259:   return(0);
2260: }

2264: /*@C
2265:   DMPlexRestoreMeet - Restore an array for the meet of the set of points

2267:   Not Collective

2269:   Input Parameters:
2270: + dm - The DMPlex object
2271: . numPoints - The number of input points for the meet
2272: - points - The input points

2274:   Output Parameters:
2275: + numCoveredPoints - The number of points in the meet
2276: - coveredPoints - The points in the meet

2278:   Level: intermediate

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

2284:   The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.

2286: .keywords: mesh
2287: .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin()
2288: @*/
2289: PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2290: {

2298:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2299:   if (numCoveredPoints) *numCoveredPoints = 0;
2300:   return(0);
2301: }

2305: /*@C
2306:   DMPlexGetFullMeet - Get an array for the meet of the set of points

2308:   Not Collective

2310:   Input Parameters:
2311: + dm - The DMPlex object
2312: . numPoints - The number of input points for the meet
2313: - points - The input points

2315:   Output Parameters:
2316: + numCoveredPoints - The number of points in the meet
2317: - coveredPoints - The points in the meet

2319:   Level: intermediate

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

2325:   The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.

2327: .keywords: mesh
2328: .seealso: DMPlexGetMeet(), DMPlexRestoreMeet(), DMPlexGetJoin()
2329: @*/
2330: PetscErrorCode DMPlexGetFullMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2331: {
2332:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2333:   PetscInt      *offsets, **closures;
2334:   PetscInt      *meet[2];
2335:   PetscInt       height = 0, maxSize, meetSize = 0, i = 0;
2336:   PetscInt       p, h, c, m, mc;


2345:   DMPlexGetDepth(dm, &height);
2346:   PetscMalloc1(numPoints, &closures);
2347:   DMGetWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2348:   mc      = mesh->maxConeSize;
2349:   maxSize = (mc > 1) ? ((PetscPowInt(mc,height+1)-1)/(mc-1)) : height + 1;
2350:   DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[0]);
2351:   DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[1]);

2353:   for (p = 0; p < numPoints; ++p) {
2354:     PetscInt closureSize;

2356:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closures[p]);

2358:     offsets[p*(height+2)+0] = 0;
2359:     for (h = 0; h < height+1; ++h) {
2360:       PetscInt pStart, pEnd, i;

2362:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
2363:       for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) {
2364:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2365:           offsets[p*(height+2)+h+1] = i;
2366:           break;
2367:         }
2368:       }
2369:       if (i == closureSize) offsets[p*(height+2)+h+1] = i;
2370:     }
2371:     if (offsets[p*(height+2)+height+1] != closureSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Total size of closure %D should be %D", offsets[p*(height+2)+height+1], closureSize);
2372:   }
2373:   for (h = 0; h < height+1; ++h) {
2374:     PetscInt dof;

2376:     /* Copy in cone of first point */
2377:     dof = offsets[h+1] - offsets[h];
2378:     for (meetSize = 0; meetSize < dof; ++meetSize) {
2379:       meet[i][meetSize] = closures[0][(offsets[h]+meetSize)*2];
2380:     }
2381:     /* Check each successive cone */
2382:     for (p = 1; p < numPoints && meetSize; ++p) {
2383:       PetscInt newMeetSize = 0;

2385:       dof = offsets[p*(height+2)+h+1] - offsets[p*(height+2)+h];
2386:       for (c = 0; c < dof; ++c) {
2387:         const PetscInt point = closures[p][(offsets[p*(height+2)+h]+c)*2];

2389:         for (m = 0; m < meetSize; ++m) {
2390:           if (point == meet[i][m]) {
2391:             meet[1-i][newMeetSize++] = point;
2392:             break;
2393:           }
2394:         }
2395:       }
2396:       meetSize = newMeetSize;
2397:       i        = 1-i;
2398:     }
2399:     if (meetSize) break;
2400:   }
2401:   *numCoveredPoints = meetSize;
2402:   *coveredPoints    = meet[i];
2403:   for (p = 0; p < numPoints; ++p) {
2404:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, NULL, &closures[p]);
2405:   }
2406:   PetscFree(closures);
2407:   DMRestoreWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2408:   DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2409:   return(0);
2410: }

2414: /*@C
2415:   DMPlexEqual - Determine if two DMs have the same topology

2417:   Not Collective

2419:   Input Parameters:
2420: + dmA - A DMPlex object
2421: - dmB - A DMPlex object

2423:   Output Parameters:
2424: . equal - PETSC_TRUE if the topologies are identical

2426:   Level: intermediate

2428:   Notes:
2429:   We are not solving graph isomorphism, so we do not permutation.

2431: .keywords: mesh
2432: .seealso: DMPlexGetCone()
2433: @*/
2434: PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal)
2435: {
2436:   PetscInt       depth, depthB, pStart, pEnd, pStartB, pEndB, p;


2444:   *equal = PETSC_FALSE;
2445:   DMPlexGetDepth(dmA, &depth);
2446:   DMPlexGetDepth(dmB, &depthB);
2447:   if (depth != depthB) return(0);
2448:   DMPlexGetChart(dmA, &pStart,  &pEnd);
2449:   DMPlexGetChart(dmB, &pStartB, &pEndB);
2450:   if ((pStart != pStartB) || (pEnd != pEndB)) return(0);
2451:   for (p = pStart; p < pEnd; ++p) {
2452:     const PetscInt *cone, *coneB, *ornt, *orntB, *support, *supportB;
2453:     PetscInt        coneSize, coneSizeB, c, supportSize, supportSizeB, s;

2455:     DMPlexGetConeSize(dmA, p, &coneSize);
2456:     DMPlexGetCone(dmA, p, &cone);
2457:     DMPlexGetConeOrientation(dmA, p, &ornt);
2458:     DMPlexGetConeSize(dmB, p, &coneSizeB);
2459:     DMPlexGetCone(dmB, p, &coneB);
2460:     DMPlexGetConeOrientation(dmB, p, &orntB);
2461:     if (coneSize != coneSizeB) return(0);
2462:     for (c = 0; c < coneSize; ++c) {
2463:       if (cone[c] != coneB[c]) return(0);
2464:       if (ornt[c] != orntB[c]) return(0);
2465:     }
2466:     DMPlexGetSupportSize(dmA, p, &supportSize);
2467:     DMPlexGetSupport(dmA, p, &support);
2468:     DMPlexGetSupportSize(dmB, p, &supportSizeB);
2469:     DMPlexGetSupport(dmB, p, &supportB);
2470:     if (supportSize != supportSizeB) return(0);
2471:     for (s = 0; s < supportSize; ++s) {
2472:       if (support[s] != supportB[s]) return(0);
2473:     }
2474:   }
2475:   *equal = PETSC_TRUE;
2476:   return(0);
2477: }

2481: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2482: {
2483:   MPI_Comm       comm;

2487:   PetscObjectGetComm((PetscObject)dm,&comm);
2489:   switch (cellDim) {
2490:   case 0:
2491:     *numFaceVertices = 0;
2492:     break;
2493:   case 1:
2494:     *numFaceVertices = 1;
2495:     break;
2496:   case 2:
2497:     switch (numCorners) {
2498:     case 3: /* triangle */
2499:       *numFaceVertices = 2; /* Edge has 2 vertices */
2500:       break;
2501:     case 4: /* quadrilateral */
2502:       *numFaceVertices = 2; /* Edge has 2 vertices */
2503:       break;
2504:     case 6: /* quadratic triangle, tri and quad cohesive Lagrange cells */
2505:       *numFaceVertices = 3; /* Edge has 3 vertices */
2506:       break;
2507:     case 9: /* quadratic quadrilateral, quadratic quad cohesive Lagrange cells */
2508:       *numFaceVertices = 3; /* Edge has 3 vertices */
2509:       break;
2510:     default:
2511:       SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %D for dimension %D", numCorners, cellDim);
2512:     }
2513:     break;
2514:   case 3:
2515:     switch (numCorners) {
2516:     case 4: /* tetradehdron */
2517:       *numFaceVertices = 3; /* Face has 3 vertices */
2518:       break;
2519:     case 6: /* tet cohesive cells */
2520:       *numFaceVertices = 4; /* Face has 4 vertices */
2521:       break;
2522:     case 8: /* hexahedron */
2523:       *numFaceVertices = 4; /* Face has 4 vertices */
2524:       break;
2525:     case 9: /* tet cohesive Lagrange cells */
2526:       *numFaceVertices = 6; /* Face has 6 vertices */
2527:       break;
2528:     case 10: /* quadratic tetrahedron */
2529:       *numFaceVertices = 6; /* Face has 6 vertices */
2530:       break;
2531:     case 12: /* hex cohesive Lagrange cells */
2532:       *numFaceVertices = 6; /* Face has 6 vertices */
2533:       break;
2534:     case 18: /* quadratic tet cohesive Lagrange cells */
2535:       *numFaceVertices = 6; /* Face has 6 vertices */
2536:       break;
2537:     case 27: /* quadratic hexahedron, quadratic hex cohesive Lagrange cells */
2538:       *numFaceVertices = 9; /* Face has 9 vertices */
2539:       break;
2540:     default:
2541:       SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %D for dimension %D", numCorners, cellDim);
2542:     }
2543:     break;
2544:   default:
2545:     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid cell dimension %D", cellDim);
2546:   }
2547:   return(0);
2548: }

2552: /*@
2553:   DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point

2555:   Not Collective

2557:   Input Parameter:
2558: . dm    - The DMPlex object

2560:   Output Parameter:
2561: . depthLabel - The DMLabel recording point depth

2563:   Level: developer

2565: .keywords: mesh, points
2566: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2567: @*/
2568: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
2569: {

2575:   if (!dm->depthLabel) {DMGetLabel(dm, "depth", &dm->depthLabel);}
2576:   *depthLabel = dm->depthLabel;
2577:   return(0);
2578: }

2582: /*@
2583:   DMPlexGetDepth - Get the depth of the DAG representing this mesh

2585:   Not Collective

2587:   Input Parameter:
2588: . dm    - The DMPlex object

2590:   Output Parameter:
2591: . depth - The number of strata (breadth first levels) in the DAG

2593:   Level: developer

2595: .keywords: mesh, points
2596: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2597: @*/
2598: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
2599: {
2600:   DMLabel        label;
2601:   PetscInt       d = 0;

2607:   DMPlexGetDepthLabel(dm, &label);
2608:   if (label) {DMLabelGetNumValues(label, &d);}
2609:   *depth = d-1;
2610:   return(0);
2611: }

2615: /*@
2616:   DMPlexGetDepthStratum - Get the bounds [start, end) for all points at a certain depth.

2618:   Not Collective

2620:   Input Parameters:
2621: + dm           - The DMPlex object
2622: - stratumValue - The requested depth

2624:   Output Parameters:
2625: + start - The first point at this depth
2626: - end   - One beyond the last point at this depth

2628:   Level: developer

2630: .keywords: mesh, points
2631: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
2632: @*/
2633: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2634: {
2635:   DMLabel        label;
2636:   PetscInt       pStart, pEnd;

2643:   DMPlexGetChart(dm, &pStart, &pEnd);
2644:   if (pStart == pEnd) return(0);
2645:   if (stratumValue < 0) {
2646:     if (start) *start = pStart;
2647:     if (end)   *end   = pEnd;
2648:     return(0);
2649:   }
2650:   DMPlexGetDepthLabel(dm, &label);
2651:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2652:   DMLabelGetStratumBounds(label, stratumValue, start, end);
2653:   return(0);
2654: }

2658: /*@
2659:   DMPlexGetHeightStratum - Get the bounds [start, end) for all points at a certain height.

2661:   Not Collective

2663:   Input Parameters:
2664: + dm           - The DMPlex object
2665: - stratumValue - The requested height

2667:   Output Parameters:
2668: + start - The first point at this height
2669: - end   - One beyond the last point at this height

2671:   Level: developer

2673: .keywords: mesh, points
2674: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
2675: @*/
2676: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2677: {
2678:   DMLabel        label;
2679:   PetscInt       depth, pStart, pEnd;

2686:   DMPlexGetChart(dm, &pStart, &pEnd);
2687:   if (pStart == pEnd) return(0);
2688:   if (stratumValue < 0) {
2689:     if (start) *start = pStart;
2690:     if (end)   *end   = pEnd;
2691:     return(0);
2692:   }
2693:   DMPlexGetDepthLabel(dm, &label);
2694:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2695:   DMLabelGetNumValues(label, &depth);
2696:   DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
2697:   return(0);
2698: }

2702: /* Set the number of dof on each point and separate by fields */
2703: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
2704: {
2705:   PetscInt      *pMax;
2706:   PetscInt       depth, pStart = 0, pEnd = 0;
2707:   PetscInt       Nf, p, d, dep, f;
2708:   PetscBool     *isFE;

2712:   PetscMalloc1(numFields, &isFE);
2713:   DMGetNumFields(dm, &Nf);
2714:   for (f = 0; f < numFields; ++f) {
2715:     PetscObject  obj;
2716:     PetscClassId id;

2718:     isFE[f] = PETSC_FALSE;
2719:     if (f >= Nf) continue;
2720:     DMGetField(dm, f, &obj);
2721:     PetscObjectGetClassId(obj, &id);
2722:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
2723:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
2724:   }
2725:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
2726:   if (numFields > 0) {
2727:     PetscSectionSetNumFields(*section, numFields);
2728:     if (numComp) {
2729:       for (f = 0; f < numFields; ++f) {
2730:         PetscSectionSetFieldComponents(*section, f, numComp[f]);
2731:       }
2732:     }
2733:   }
2734:   DMPlexGetChart(dm, &pStart, &pEnd);
2735:   PetscSectionSetChart(*section, pStart, pEnd);
2736:   DMPlexGetDepth(dm, &depth);
2737:   PetscMalloc1(depth+1,&pMax);
2738:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2739:   for (dep = 0; dep <= depth; ++dep) {
2740:     d    = dim == depth ? dep : (!dep ? 0 : dim);
2741:     DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
2742:     pMax[dep] = pMax[dep] < 0 ? pEnd : pMax[dep];
2743:     for (p = pStart; p < pEnd; ++p) {
2744:       PetscInt tot = 0;

2746:       for (f = 0; f < numFields; ++f) {
2747:         if (isFE[f] && p >= pMax[dep]) continue;
2748:         PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
2749:         tot += numDof[f*(dim+1)+d];
2750:       }
2751:       PetscSectionSetDof(*section, p, tot);
2752:     }
2753:   }
2754:   PetscFree(pMax);
2755:   PetscFree(isFE);
2756:   return(0);
2757: }

2761: /* Set the number of dof on each point and separate by fields
2762:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
2763: */
2764: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2765: {
2766:   PetscInt       numFields;
2767:   PetscInt       bc;
2768:   PetscSection   aSec;

2772:   PetscSectionGetNumFields(section, &numFields);
2773:   for (bc = 0; bc < numBC; ++bc) {
2774:     PetscInt        field = 0;
2775:     const PetscInt *comp;
2776:     const PetscInt *idx;
2777:     PetscInt        Nc = -1, n, i;

2779:     if (numFields) field = bcField[bc];
2780:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2781:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2782:     ISGetLocalSize(bcPoints[bc], &n);
2783:     ISGetIndices(bcPoints[bc], &idx);
2784:     for (i = 0; i < n; ++i) {
2785:       const PetscInt p = idx[i];
2786:       PetscInt       numConst;

2788:       if (numFields) {
2789:         PetscSectionGetFieldDof(section, p, field, &numConst);
2790:       } else {
2791:         PetscSectionGetDof(section, p, &numConst);
2792:       }
2793:       /* If Nc < 0, constrain every dof on the point */
2794:       if (Nc > 0) numConst = PetscMin(numConst, Nc);
2795:       if (numFields) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
2796:       PetscSectionAddConstraintDof(section, p, numConst);
2797:     }
2798:     ISRestoreIndices(bcPoints[bc], &idx);
2799:     if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2800:   }
2801:   DMPlexGetAnchors(dm, &aSec, NULL);
2802:   if (aSec) {
2803:     PetscInt aStart, aEnd, a;

2805:     PetscSectionGetChart(aSec, &aStart, &aEnd);
2806:     for (a = aStart; a < aEnd; a++) {
2807:       PetscInt dof, f;

2809:       PetscSectionGetDof(aSec, a, &dof);
2810:       if (dof) {
2811:         /* if there are point-to-point constraints, then all dofs are constrained */
2812:         PetscSectionGetDof(section, a, &dof);
2813:         PetscSectionSetConstraintDof(section, a, dof);
2814:         for (f = 0; f < numFields; f++) {
2815:           PetscSectionGetFieldDof(section, a, f, &dof);
2816:           PetscSectionSetFieldConstraintDof(section, a, f, dof);
2817:         }
2818:       }
2819:     }
2820:   }
2821:   return(0);
2822: }

2826: /* Set the constrained field indices on each point
2827:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
2828: */
2829: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2830: {
2831:   PetscSection   aSec;
2832:   PetscInt      *indices;
2833:   PetscInt       numFields, maxDof, pStart, pEnd, p, bc, f, d;

2837:   PetscSectionGetNumFields(section, &numFields);
2838:   if (!numFields) return(0);
2839:   /* Initialize all field indices to -1 */
2840:   PetscSectionGetChart(section, &pStart, &pEnd);
2841:   PetscSectionGetMaxDof(section, &maxDof);
2842:   PetscMalloc1(maxDof, &indices);
2843:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
2844:   for (p = pStart; p < pEnd; ++p) for (f = 0; f < numFields; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
2845:   /* Handle BC constraints */
2846:   for (bc = 0; bc < numBC; ++bc) {
2847:     const PetscInt  field = bcField[bc];
2848:     const PetscInt *comp, *idx;
2849:     PetscInt        Nc = -1, n, i;

2851:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2852:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2853:     ISGetLocalSize(bcPoints[bc], &n);
2854:     ISGetIndices(bcPoints[bc], &idx);
2855:     for (i = 0; i < n; ++i) {
2856:       const PetscInt  p = idx[i];
2857:       const PetscInt *find;
2858:       PetscInt        fcdof, c;

2860:       PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
2861:       if (Nc < 0) {
2862:         for (d = 0; d < fcdof; ++d) indices[d] = d;
2863:       } else {
2864:         PetscSectionGetFieldConstraintIndices(section, p, field, &find);
2865:         for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
2866:         for (c = 0; c < Nc; ++c) indices[d+c] = comp[c];
2867:         PetscSortInt(d+Nc, indices);
2868:         for (c = d+Nc; c < fcdof; ++c) indices[c] = -1;
2869:       }
2870:       PetscSectionSetFieldConstraintIndices(section, p, field, indices);
2871:     }
2872:     if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2873:     ISRestoreIndices(bcPoints[bc], &idx);
2874:   }
2875:   /* Handle anchors */
2876:   DMPlexGetAnchors(dm, &aSec, NULL);
2877:   if (aSec) {
2878:     PetscInt aStart, aEnd, a;

2880:     for (d = 0; d < maxDof; ++d) indices[d] = d;
2881:     PetscSectionGetChart(aSec, &aStart, &aEnd);
2882:     for (a = aStart; a < aEnd; a++) {
2883:       PetscInt dof, fdof, f;

2885:       PetscSectionGetDof(aSec, a, &dof);
2886:       if (dof) {
2887:         /* if there are point-to-point constraints, then all dofs are constrained */
2888:         for (f = 0; f < numFields; f++) {
2889:           PetscSectionGetFieldDof(section, a, f, &fdof);
2890:           PetscSectionSetFieldConstraintIndices(section, a, f, indices);
2891:         }
2892:       }
2893:     }
2894:   }
2895:   PetscFree(indices);
2896:   return(0);
2897: }

2901: /* Set the constrained indices on each point */
2902: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
2903: {
2904:   PetscInt      *indices;
2905:   PetscInt       numFields, maxDof, pStart, pEnd, p, f, d;

2909:   PetscSectionGetNumFields(section, &numFields);
2910:   PetscSectionGetMaxDof(section, &maxDof);
2911:   PetscSectionGetChart(section, &pStart, &pEnd);
2912:   PetscMalloc1(maxDof, &indices);
2913:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
2914:   for (p = pStart; p < pEnd; ++p) {
2915:     PetscInt cdof, d;

2917:     PetscSectionGetConstraintDof(section, p, &cdof);
2918:     if (cdof) {
2919:       if (numFields) {
2920:         PetscInt numConst = 0, foff = 0;

2922:         for (f = 0; f < numFields; ++f) {
2923:           const PetscInt *find;
2924:           PetscInt        fcdof, fdof;

2926:           PetscSectionGetFieldDof(section, p, f, &fdof);
2927:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
2928:           /* Change constraint numbering from field component to local dof number */
2929:           PetscSectionGetFieldConstraintIndices(section, p, f, &find);
2930:           for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
2931:           numConst += fcdof;
2932:           foff     += fdof;
2933:         }
2934:         if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
2935:       } else {
2936:         for (d = 0; d < cdof; ++d) indices[d] = d;
2937:       }
2938:       PetscSectionSetConstraintIndices(section, p, indices);
2939:     }
2940:   }
2941:   PetscFree(indices);
2942:   return(0);
2943: }

2947: /*@C
2948:   DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided.

2950:   Not Collective

2952:   Input Parameters:
2953: + dm        - The DMPlex object
2954: . dim       - The spatial dimension of the problem
2955: . numFields - The number of fields in the problem
2956: . numComp   - An array of size numFields that holds the number of components for each field
2957: . numDof    - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
2958: . numBC     - The number of boundary conditions
2959: . bcField   - An array of size numBC giving the field number for each boundry condition
2960: . bcComps   - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
2961: . bcPoints  - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
2962: - perm      - Optional permutation of the chart, or NULL

2964:   Output Parameter:
2965: . section - The PetscSection object

2967:   Notes: numDof[f*(dim+1)+d] gives the number of dof for field f on sieve points of dimension d. For instance, numDof[1] is the
2968:   number of dof for field 0 on each edge.

2970:   The chart permutation is the same one set using PetscSectionSetPermutation()

2972:   Level: developer

2974:   Fortran Notes:
2975:   A Fortran 90 version is available as DMPlexCreateSectionF90()

2977: .keywords: mesh, elements
2978: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
2979: @*/
2980: PetscErrorCode DMPlexCreateSection(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], IS perm, PetscSection *section)
2981: {
2982:   PetscSection   aSec;

2986:   DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
2987:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
2988:   if (perm) {PetscSectionSetPermutation(*section, perm);}
2989:   PetscSectionSetUp(*section);
2990:   DMPlexGetAnchors(dm,&aSec,NULL);
2991:   if (numBC || aSec) {
2992:     DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
2993:     DMPlexCreateSectionBCIndices(dm, *section);
2994:   }
2995:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
2996:   return(0);
2997: }

3001: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
3002: {
3003:   PetscSection   section, s;
3004:   Mat            m;

3008:   DMClone(dm, cdm);
3009:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
3010:   DMSetDefaultSection(*cdm, section);
3011:   PetscSectionDestroy(&section);
3012:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3013:   MatCreate(PETSC_COMM_SELF, &m);
3014:   DMSetDefaultConstraints(*cdm, s, m);
3015:   PetscSectionDestroy(&s);
3016:   MatDestroy(&m);
3017:   return(0);
3018: }

3022: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
3023: {
3024:   DM_Plex *mesh = (DM_Plex*) dm->data;

3028:   if (section) *section = mesh->coneSection;
3029:   return(0);
3030: }

3034: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
3035: {
3036:   DM_Plex *mesh = (DM_Plex*) dm->data;

3040:   if (section) *section = mesh->supportSection;
3041:   return(0);
3042: }

3046: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
3047: {
3048:   DM_Plex *mesh = (DM_Plex*) dm->data;

3052:   if (cones) *cones = mesh->cones;
3053:   return(0);
3054: }

3058: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
3059: {
3060:   DM_Plex *mesh = (DM_Plex*) dm->data;

3064:   if (coneOrientations) *coneOrientations = mesh->coneOrientations;
3065:   return(0);
3066: }

3068: /******************************** FEM Support **********************************/

3072: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3073: {
3074:   PetscScalar    *array, *vArray;
3075:   const PetscInt *cone, *coneO;
3076:   PetscInt        pStart, pEnd, p, numPoints, size = 0, offset = 0;
3077:   PetscErrorCode  ierr;

3080:   PetscSectionGetChart(section, &pStart, &pEnd);
3081:   DMPlexGetConeSize(dm, point, &numPoints);
3082:   DMPlexGetCone(dm, point, &cone);
3083:   DMPlexGetConeOrientation(dm, point, &coneO);
3084:   if (!values || !*values) {
3085:     if ((point >= pStart) && (point < pEnd)) {
3086:       PetscInt dof;

3088:       PetscSectionGetDof(section, point, &dof);
3089:       size += dof;
3090:     }
3091:     for (p = 0; p < numPoints; ++p) {
3092:       const PetscInt cp = cone[p];
3093:       PetscInt       dof;

3095:       if ((cp < pStart) || (cp >= pEnd)) continue;
3096:       PetscSectionGetDof(section, cp, &dof);
3097:       size += dof;
3098:     }
3099:     if (!values) {
3100:       if (csize) *csize = size;
3101:       return(0);
3102:     }
3103:     DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
3104:   } else {
3105:     array = *values;
3106:   }
3107:   size = 0;
3108:   VecGetArray(v, &vArray);
3109:   if ((point >= pStart) && (point < pEnd)) {
3110:     PetscInt     dof, off, d;
3111:     PetscScalar *varr;

3113:     PetscSectionGetDof(section, point, &dof);
3114:     PetscSectionGetOffset(section, point, &off);
3115:     varr = &vArray[off];
3116:     for (d = 0; d < dof; ++d, ++offset) {
3117:       array[offset] = varr[d];
3118:     }
3119:     size += dof;
3120:   }
3121:   for (p = 0; p < numPoints; ++p) {
3122:     const PetscInt cp = cone[p];
3123:     PetscInt       o  = coneO[p];
3124:     PetscInt       dof, off, d;
3125:     PetscScalar   *varr;

3127:     if ((cp < pStart) || (cp >= pEnd)) continue;
3128:     PetscSectionGetDof(section, cp, &dof);
3129:     PetscSectionGetOffset(section, cp, &off);
3130:     varr = &vArray[off];
3131:     if (o >= 0) {
3132:       for (d = 0; d < dof; ++d, ++offset) {
3133:         array[offset] = varr[d];
3134:       }
3135:     } else {
3136:       for (d = dof-1; d >= 0; --d, ++offset) {
3137:         array[offset] = varr[d];
3138:       }
3139:     }
3140:     size += dof;
3141:   }
3142:   VecRestoreArray(v, &vArray);
3143:   if (!*values) {
3144:     if (csize) *csize = size;
3145:     *values = array;
3146:   } else {
3147:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3148:     *csize = size;
3149:   }
3150:   return(0);
3151: }

3155: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3156: {
3157:   PetscInt       offset = 0, p;

3161:   *size = 0;
3162:   for (p = 0; p < numPoints*2; p += 2) {
3163:     const PetscInt point = points[p];
3164:     const PetscInt o     = points[p+1];
3165:     PetscInt       dof, off, d;
3166:     const PetscScalar *varr;

3168:     PetscSectionGetDof(section, point, &dof);
3169:     PetscSectionGetOffset(section, point, &off);
3170:     varr = &vArray[off];
3171:     if (o >= 0) {
3172:       for (d = 0; d < dof; ++d, ++offset)    array[offset] = varr[d];
3173:     } else {
3174:       for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
3175:     }
3176:   }
3177:   *size = offset;
3178:   return(0);
3179: }

3183: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Fields_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], PetscInt numFields, const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3184: {
3185:   PetscInt       offset = 0, f;

3189:   *size = 0;
3190:   for (f = 0; f < numFields; ++f) {
3191:     PetscInt fcomp, p;

3193:     PetscSectionGetFieldComponents(section, f, &fcomp);
3194:     for (p = 0; p < numPoints*2; p += 2) {
3195:       const PetscInt point = points[p];
3196:       const PetscInt o     = points[p+1];
3197:       PetscInt       fdof, foff, d, c;
3198:       const PetscScalar *varr;

3200:       PetscSectionGetFieldDof(section, point, f, &fdof);
3201:       PetscSectionGetFieldOffset(section, point, f, &foff);
3202:       varr = &vArray[foff];
3203:       if (o >= 0) {
3204:         for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
3205:       } else {
3206:         for (d = fdof/fcomp-1; d >= 0; --d) {
3207:           for (c = 0; c < fcomp; ++c, ++offset) {
3208:             array[offset] = varr[d*fcomp+c];
3209:           }
3210:         }
3211:       }
3212:     }
3213:   }
3214:   *size = offset;
3215:   return(0);
3216: }

3220: /*@C
3221:   DMPlexVecGetClosure - Get an array of the values on the closure of 'point'

3223:   Not collective

3225:   Input Parameters:
3226: + dm - The DM
3227: . section - The section describing the layout in v, or NULL to use the default section
3228: . v - The local vector
3229: - point - The sieve point in the DM

3231:   Output Parameters:
3232: + csize - The number of values in the closure, or NULL
3233: - values - The array of values, which is a borrowed array and should not be freed

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

3239:   The csize argument is not present in the Fortran 90 binding since it is internal to the array.

3241:   Level: intermediate

3243: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3244: @*/
3245: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3246: {
3247:   PetscSection    clSection;
3248:   IS              clPoints;
3249:   PetscScalar    *array, *vArray;
3250:   PetscInt       *points = NULL;
3251:   const PetscInt *clp;
3252:   PetscInt        depth, numFields, numPoints, size;
3253:   PetscErrorCode  ierr;

3257:   if (!section) {DMGetDefaultSection(dm, &section);}
3260:   DMPlexGetDepth(dm, &depth);
3261:   PetscSectionGetNumFields(section, &numFields);
3262:   if (depth == 1 && numFields < 2) {
3263:     DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
3264:     return(0);
3265:   }
3266:   /* Get points */
3267:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3268:   if (!clPoints) {
3269:     PetscInt pStart, pEnd, p, q;

3271:     PetscSectionGetChart(section, &pStart, &pEnd);
3272:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3273:     /* Compress out points not in the section */
3274:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3275:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3276:         points[q*2]   = points[p];
3277:         points[q*2+1] = points[p+1];
3278:         ++q;
3279:       }
3280:     }
3281:     numPoints = q;
3282:   } else {
3283:     PetscInt dof, off;

3285:     PetscSectionGetDof(clSection, point, &dof);
3286:     PetscSectionGetOffset(clSection, point, &off);
3287:     ISGetIndices(clPoints, &clp);
3288:     numPoints = dof/2;
3289:     points    = (PetscInt *) &clp[off];
3290:   }
3291:   /* Get array */
3292:   if (!values || !*values) {
3293:     PetscInt asize = 0, dof, p;

3295:     for (p = 0; p < numPoints*2; p += 2) {
3296:       PetscSectionGetDof(section, points[p], &dof);
3297:       asize += dof;
3298:     }
3299:     if (!values) {
3300:       if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3301:       else           {ISRestoreIndices(clPoints, &clp);}
3302:       if (csize) *csize = asize;
3303:       return(0);
3304:     }
3305:     DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
3306:   } else {
3307:     array = *values;
3308:   }
3309:   VecGetArray(v, &vArray);
3310:   /* Get values */
3311:   if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);}
3312:   else               {DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);}
3313:   /* Cleanup points */
3314:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3315:   else           {ISRestoreIndices(clPoints, &clp);}
3316:   /* Cleanup array */
3317:   VecRestoreArray(v, &vArray);
3318:   if (!*values) {
3319:     if (csize) *csize = size;
3320:     *values = array;
3321:   } else {
3322:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %D < actual size %D", *csize, size);
3323:     *csize = size;
3324:   }
3325:   return(0);
3326: }

3330: /*@C
3331:   DMPlexVecRestoreClosure - Restore the array of the values on the closure of 'point'

3333:   Not collective

3335:   Input Parameters:
3336: + dm - The DM
3337: . section - The section describing the layout in v, or NULL to use the default section
3338: . v - The local vector
3339: . point - The sieve point in the DM
3340: . csize - The number of values in the closure, or NULL
3341: - values - The array of values, which is a borrowed array and should not be freed

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

3347:   The csize argument is not present in the Fortran 90 binding since it is internal to the array.

3349:   Level: intermediate

3351: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3352: @*/
3353: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3354: {
3355:   PetscInt       size = 0;

3359:   /* Should work without recalculating size */
3360:   DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
3361:   return(0);
3362: }

3364: PETSC_STATIC_INLINE void add   (PetscScalar *x, PetscScalar y) {*x += y;}
3365: PETSC_STATIC_INLINE void insert(PetscScalar *x, PetscScalar y) {*x  = y;}

3369: PETSC_STATIC_INLINE PetscErrorCode updatePoint_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3370: {
3371:   PetscInt        cdof;   /* The number of constraints on this point */
3372:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3373:   PetscScalar    *a;
3374:   PetscInt        off, cind = 0, k;
3375:   PetscErrorCode  ierr;

3378:   PetscSectionGetConstraintDof(section, point, &cdof);
3379:   PetscSectionGetOffset(section, point, &off);
3380:   a    = &array[off];
3381:   if (!cdof || setBC) {
3382:     if (orientation >= 0) {
3383:       for (k = 0; k < dof; ++k) {
3384:         fuse(&a[k], values[k]);
3385:       }
3386:     } else {
3387:       for (k = 0; k < dof; ++k) {
3388:         fuse(&a[k], values[dof-k-1]);
3389:       }
3390:     }
3391:   } else {
3392:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3393:     if (orientation >= 0) {
3394:       for (k = 0; k < dof; ++k) {
3395:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3396:         fuse(&a[k], values[k]);
3397:       }
3398:     } else {
3399:       for (k = 0; k < dof; ++k) {
3400:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3401:         fuse(&a[k], values[dof-k-1]);
3402:       }
3403:     }
3404:   }
3405:   return(0);
3406: }

3410: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3411: {
3412:   PetscInt        cdof;   /* The number of constraints on this point */
3413:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3414:   PetscScalar    *a;
3415:   PetscInt        off, cind = 0, k;
3416:   PetscErrorCode  ierr;

3419:   PetscSectionGetConstraintDof(section, point, &cdof);
3420:   PetscSectionGetOffset(section, point, &off);
3421:   a    = &array[off];
3422:   if (cdof) {
3423:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3424:     if (orientation >= 0) {
3425:       for (k = 0; k < dof; ++k) {
3426:         if ((cind < cdof) && (k == cdofs[cind])) {
3427:           fuse(&a[k], values[k]);
3428:           ++cind;
3429:         }
3430:       }
3431:     } else {
3432:       for (k = 0; k < dof; ++k) {
3433:         if ((cind < cdof) && (k == cdofs[cind])) {
3434:           fuse(&a[k], values[dof-k-1]);
3435:           ++cind;
3436:         }
3437:       }
3438:     }
3439:   }
3440:   return(0);
3441: }

3445: PETSC_STATIC_INLINE PetscErrorCode updatePointFields_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, const PetscScalar values[], PetscInt *offset, PetscScalar array[])
3446: {
3447:   PetscScalar    *a;
3448:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3449:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3450:   PetscInt        cind = 0, k, c;
3451:   PetscErrorCode  ierr;

3454:   PetscSectionGetFieldDof(section, point, f, &fdof);
3455:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3456:   PetscSectionGetFieldOffset(section, point, f, &foff);
3457:   a    = &array[foff];
3458:   if (!fcdof || setBC) {
3459:     if (o >= 0) {
3460:       for (k = 0; k < fdof; ++k) fuse(&a[k], values[foffset+k]);
3461:     } else {
3462:       for (k = fdof/fcomp-1; k >= 0; --k) {
3463:         for (c = 0; c < fcomp; ++c) {
3464:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3465:         }
3466:       }
3467:     }
3468:   } else {
3469:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3470:     if (o >= 0) {
3471:       for (k = 0; k < fdof; ++k) {
3472:         if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
3473:         fuse(&a[k], values[foffset+k]);
3474:       }
3475:     } else {
3476:       for (k = fdof/fcomp-1; k >= 0; --k) {
3477:         for (c = 0; c < fcomp; ++c) {
3478:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
3479:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3480:         }
3481:       }
3482:     }
3483:   }
3484:   *offset += fdof;
3485:   return(0);
3486: }

3490: PETSC_STATIC_INLINE PetscErrorCode updatePointFieldsBC_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), const PetscScalar values[], PetscInt *offset, PetscScalar array[])
3491: {
3492:   PetscScalar    *a;
3493:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3494:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3495:   PetscInt        cind = 0, k, c;
3496:   PetscErrorCode  ierr;

3499:   PetscSectionGetFieldDof(section, point, f, &fdof);
3500:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3501:   PetscSectionGetFieldOffset(section, point, f, &foff);
3502:   a    = &array[foff];
3503:   if (fcdof) {
3504:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3505:     if (o >= 0) {
3506:       for (k = 0; k < fdof; ++k) {
3507:         if ((cind < fcdof) && (k == fcdofs[cind])) {
3508:           fuse(&a[k], values[foffset+k]);
3509:           ++cind;
3510:         }
3511:       }
3512:     } else {
3513:       for (k = fdof/fcomp-1; k >= 0; --k) {
3514:         for (c = 0; c < fcomp; ++c) {
3515:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {
3516:             fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3517:             ++cind;
3518:           }
3519:         }
3520:       }
3521:     }
3522:   }
3523:   *offset += fdof;
3524:   return(0);
3525: }

3529: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3530: {
3531:   PetscScalar    *array;
3532:   const PetscInt *cone, *coneO;
3533:   PetscInt        pStart, pEnd, p, numPoints, off, dof;
3534:   PetscErrorCode  ierr;

3537:   PetscSectionGetChart(section, &pStart, &pEnd);
3538:   DMPlexGetConeSize(dm, point, &numPoints);
3539:   DMPlexGetCone(dm, point, &cone);
3540:   DMPlexGetConeOrientation(dm, point, &coneO);
3541:   VecGetArray(v, &array);
3542:   for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
3543:     const PetscInt cp = !p ? point : cone[p-1];
3544:     const PetscInt o  = !p ? 0     : coneO[p-1];

3546:     if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
3547:     PetscSectionGetDof(section, cp, &dof);
3548:     /* ADD_VALUES */
3549:     {
3550:       const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3551:       PetscScalar    *a;
3552:       PetscInt        cdof, coff, cind = 0, k;

3554:       PetscSectionGetConstraintDof(section, cp, &cdof);
3555:       PetscSectionGetOffset(section, cp, &coff);
3556:       a    = &array[coff];
3557:       if (!cdof) {
3558:         if (o >= 0) {
3559:           for (k = 0; k < dof; ++k) {
3560:             a[k] += values[off+k];
3561:           }
3562:         } else {
3563:           for (k = 0; k < dof; ++k) {
3564:             a[k] += values[off+dof-k-1];
3565:           }
3566:         }
3567:       } else {
3568:         PetscSectionGetConstraintIndices(section, cp, &cdofs);
3569:         if (o >= 0) {
3570:           for (k = 0; k < dof; ++k) {
3571:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3572:             a[k] += values[off+k];
3573:           }
3574:         } else {
3575:           for (k = 0; k < dof; ++k) {
3576:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3577:             a[k] += values[off+dof-k-1];
3578:           }
3579:         }
3580:       }
3581:     }
3582:   }
3583:   VecRestoreArray(v, &array);
3584:   return(0);
3585: }

3589: /*@C
3590:   DMPlexVecSetClosure - Set an array of the values on the closure of 'point'

3592:   Not collective

3594:   Input Parameters:
3595: + dm - The DM
3596: . section - The section describing the layout in v, or NULL to use the default section
3597: . v - The local vector
3598: . point - The sieve point in the DM
3599: . values - The array of values
3600: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

3602:   Fortran Notes:
3603:   This routine is only available in Fortran 90, and you must include petsc.h90 in your code.

3605:   Level: intermediate

3607: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
3608: @*/
3609: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3610: {
3611:   PetscSection    clSection;
3612:   IS              clPoints;
3613:   PetscScalar    *array;
3614:   PetscInt       *points = NULL;
3615:   const PetscInt *clp;
3616:   PetscInt        depth, numFields, numPoints, p;
3617:   PetscErrorCode  ierr;

3621:   if (!section) {DMGetDefaultSection(dm, &section);}
3624:   DMPlexGetDepth(dm, &depth);
3625:   PetscSectionGetNumFields(section, &numFields);
3626:   if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
3627:     DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
3628:     return(0);
3629:   }
3630:   /* Get points */
3631:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3632:   if (!clPoints) {
3633:     PetscInt pStart, pEnd, q;

3635:     PetscSectionGetChart(section, &pStart, &pEnd);
3636:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3637:     /* Compress out points not in the section */
3638:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3639:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3640:         points[q*2]   = points[p];
3641:         points[q*2+1] = points[p+1];
3642:         ++q;
3643:       }
3644:     }
3645:     numPoints = q;
3646:   } else {
3647:     PetscInt dof, off;

3649:     PetscSectionGetDof(clSection, point, &dof);
3650:     PetscSectionGetOffset(clSection, point, &off);
3651:     ISGetIndices(clPoints, &clp);
3652:     numPoints = dof/2;
3653:     points    = (PetscInt *) &clp[off];
3654:   }
3655:   /* Get array */
3656:   VecGetArray(v, &array);
3657:   /* Get values */
3658:   if (numFields > 0) {
3659:     PetscInt offset = 0, fcomp, f;
3660:     for (f = 0; f < numFields; ++f) {
3661:       PetscSectionGetFieldComponents(section, f, &fcomp);
3662:       switch (mode) {
3663:       case INSERT_VALUES:
3664:         for (p = 0; p < numPoints*2; p += 2) {
3665:           const PetscInt point = points[p];
3666:           const PetscInt o     = points[p+1];
3667:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3668:         } break;
3669:       case INSERT_ALL_VALUES:
3670:         for (p = 0; p < numPoints*2; p += 2) {
3671:           const PetscInt point = points[p];
3672:           const PetscInt o     = points[p+1];
3673:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3674:         } break;
3675:       case INSERT_BC_VALUES:
3676:         for (p = 0; p < numPoints*2; p += 2) {
3677:           const PetscInt point = points[p];
3678:           const PetscInt o     = points[p+1];
3679:           updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3680:         } break;
3681:       case ADD_VALUES:
3682:         for (p = 0; p < numPoints*2; p += 2) {
3683:           const PetscInt point = points[p];
3684:           const PetscInt o     = points[p+1];
3685:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3686:         } break;
3687:       case ADD_ALL_VALUES:
3688:         for (p = 0; p < numPoints*2; p += 2) {
3689:           const PetscInt point = points[p];
3690:           const PetscInt o     = points[p+1];
3691:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3692:         } break;
3693:       default:
3694:         SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3695:       }
3696:     }
3697:   } else {
3698:     PetscInt dof, off;

3700:     switch (mode) {
3701:     case INSERT_VALUES:
3702:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3703:         PetscInt o = points[p+1];
3704:         PetscSectionGetDof(section, points[p], &dof);
3705:         updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
3706:       } break;
3707:     case INSERT_ALL_VALUES:
3708:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3709:         PetscInt o = points[p+1];
3710:         PetscSectionGetDof(section, points[p], &dof);
3711:         updatePoint_private(section, points[p], dof, insert, PETSC_TRUE,  o, &values[off], array);
3712:       } break;
3713:     case INSERT_BC_VALUES:
3714:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3715:         PetscInt o = points[p+1];
3716:         PetscSectionGetDof(section, points[p], &dof);
3717:         updatePointBC_private(section, points[p], dof, insert,  o, &values[off], array);
3718:       } break;
3719:     case ADD_VALUES:
3720:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3721:         PetscInt o = points[p+1];
3722:         PetscSectionGetDof(section, points[p], &dof);
3723:         updatePoint_private(section, points[p], dof, add,    PETSC_FALSE, o, &values[off], array);
3724:       } break;
3725:     case ADD_ALL_VALUES:
3726:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3727:         PetscInt o = points[p+1];
3728:         PetscSectionGetDof(section, points[p], &dof);
3729:         updatePoint_private(section, points[p], dof, add,    PETSC_TRUE,  o, &values[off], array);
3730:       } break;
3731:     default:
3732:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3733:     }
3734:   }
3735:   /* Cleanup points */
3736:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3737:   else           {ISRestoreIndices(clPoints, &clp);}
3738:   /* Cleanup array */
3739:   VecRestoreArray(v, &array);
3740:   return(0);
3741: }

3745: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
3746: {
3747:   PetscSection    clSection;
3748:   IS              clPoints;
3749:   PetscScalar    *array;
3750:   PetscInt       *points = NULL;
3751:   const PetscInt *clp;
3752:   PetscInt        numFields, numPoints, p;
3753:   PetscInt        offset = 0, fcomp, f;
3754:   PetscErrorCode  ierr;

3758:   if (!section) {DMGetDefaultSection(dm, &section);}
3761:   PetscSectionGetNumFields(section, &numFields);
3762:   /* Get points */
3763:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3764:   if (!clPoints) {
3765:     PetscInt pStart, pEnd, q;

3767:     PetscSectionGetChart(section, &pStart, &pEnd);
3768:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3769:     /* Compress out points not in the section */
3770:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3771:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3772:         points[q*2]   = points[p];
3773:         points[q*2+1] = points[p+1];
3774:         ++q;
3775:       }
3776:     }
3777:     numPoints = q;
3778:   } else {
3779:     PetscInt dof, off;

3781:     PetscSectionGetDof(clSection, point, &dof);
3782:     PetscSectionGetOffset(clSection, point, &off);
3783:     ISGetIndices(clPoints, &clp);
3784:     numPoints = dof/2;
3785:     points    = (PetscInt *) &clp[off];
3786:   }
3787:   /* Get array */
3788:   VecGetArray(v, &array);
3789:   /* Get values */
3790:   for (f = 0; f < numFields; ++f) {
3791:     PetscSectionGetFieldComponents(section, f, &fcomp);
3792:     if (!fieldActive[f]) {
3793:       for (p = 0; p < numPoints*2; p += 2) {
3794:         PetscInt fdof;
3795:         PetscSectionGetFieldDof(section, points[p], f, &fdof);
3796:         offset += fdof;
3797:       }
3798:       continue;
3799:     }
3800:     switch (mode) {
3801:     case INSERT_VALUES:
3802:       for (p = 0; p < numPoints*2; p += 2) {
3803:         const PetscInt point = points[p];
3804:         const PetscInt o     = points[p+1];
3805:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3806:       } break;
3807:     case INSERT_ALL_VALUES:
3808:       for (p = 0; p < numPoints*2; p += 2) {
3809:         const PetscInt point = points[p];
3810:         const PetscInt o     = points[p+1];
3811:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3812:         } break;
3813:     case INSERT_BC_VALUES:
3814:       for (p = 0; p < numPoints*2; p += 2) {
3815:         const PetscInt point = points[p];
3816:         const PetscInt o     = points[p+1];
3817:         updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3818:       } break;
3819:     case ADD_VALUES:
3820:       for (p = 0; p < numPoints*2; p += 2) {
3821:         const PetscInt point = points[p];
3822:         const PetscInt o     = points[p+1];
3823:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3824:       } break;
3825:     case ADD_ALL_VALUES:
3826:       for (p = 0; p < numPoints*2; p += 2) {
3827:         const PetscInt point = points[p];
3828:         const PetscInt o     = points[p+1];
3829:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3830:       } break;
3831:     default:
3832:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3833:     }
3834:   }
3835:   /* Cleanup points */
3836:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3837:   else           {ISRestoreIndices(clPoints, &clp);}
3838:   /* Cleanup array */
3839:   VecRestoreArray(v, &array);
3840:   return(0);
3841: }

3845: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
3846: {
3847:   PetscMPIInt    rank;
3848:   PetscInt       i, j;

3852:   MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
3853:   PetscViewerASCIIPrintf(viewer, "[%d]mat for sieve point %D\n", rank, point);
3854:   for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
3855:   for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
3856:   numCIndices = numCIndices ? numCIndices : numRIndices;
3857:   for (i = 0; i < numRIndices; i++) {
3858:     PetscViewerASCIIPrintf(viewer, "[%d]", rank);
3859:     for (j = 0; j < numCIndices; j++) {
3860: #if defined(PETSC_USE_COMPLEX)
3861:       PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
3862: #else
3863:       PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
3864: #endif
3865:     }
3866:     PetscViewerASCIIPrintf(viewer, "\n");
3867:   }
3868:   return(0);
3869: }

3873: /* . off - The global offset of this point */
3874: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
3875: {
3876:   PetscInt        dof;    /* The number of unknowns on this point */
3877:   PetscInt        cdof;   /* The number of constraints on this point */
3878:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3879:   PetscInt        cind = 0, k;
3880:   PetscErrorCode  ierr;

3883:   PetscSectionGetDof(section, point, &dof);
3884:   PetscSectionGetConstraintDof(section, point, &cdof);
3885:   if (!cdof || setBC) {
3886:     if (orientation >= 0) {
3887:       for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
3888:     } else {
3889:       for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
3890:     }
3891:   } else {
3892:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3893:     if (orientation >= 0) {
3894:       for (k = 0; k < dof; ++k) {
3895:         if ((cind < cdof) && (k == cdofs[cind])) {
3896:           /* Insert check for returning constrained indices */
3897:           indices[*loff+k] = -(off+k+1);
3898:           ++cind;
3899:         } else {
3900:           indices[*loff+k] = off+k-cind;
3901:         }
3902:       }
3903:     } else {
3904:       for (k = 0; k < dof; ++k) {
3905:         if ((cind < cdof) && (k == cdofs[cind])) {
3906:           /* Insert check for returning constrained indices */
3907:           indices[*loff+dof-k-1] = -(off+k+1);
3908:           ++cind;
3909:         } else {
3910:           indices[*loff+dof-k-1] = off+k-cind;
3911:         }
3912:       }
3913:     }
3914:   }
3915:   *loff += dof;
3916:   return(0);
3917: }

3921: /* . off - The global offset of this point */
3922: PetscErrorCode indicesPointFields_private(PetscSection section, PetscInt point, PetscInt off, PetscInt foffs[], PetscBool setBC, PetscInt orientation, PetscInt indices[])
3923: {
3924:   PetscInt       numFields, foff, f;

3928:   PetscSectionGetNumFields(section, &numFields);
3929:   for (f = 0, foff = 0; f < numFields; ++f) {
3930:     PetscInt        fdof, fcomp, cfdof;
3931:     const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3932:     PetscInt        cind = 0, k, c;

3934:     PetscSectionGetFieldComponents(section, f, &fcomp);
3935:     PetscSectionGetFieldDof(section, point, f, &fdof);
3936:     PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
3937:     if (!cfdof || setBC) {
3938:       if (orientation >= 0) {
3939:         for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
3940:       } else {
3941:         for (k = fdof/fcomp-1; k >= 0; --k) {
3942:           for (c = 0; c < fcomp; ++c) {
3943:             indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
3944:           }
3945:         }
3946:       }
3947:     } else {
3948:       PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3949:       if (orientation >= 0) {
3950:         for (k = 0; k < fdof; ++k) {
3951:           if ((cind < cfdof) && (k == fcdofs[cind])) {
3952:             indices[foffs[f]+k] = -(off+foff+k+1);
3953:             ++cind;
3954:           } else {
3955:             indices[foffs[f]+k] = off+foff+k-cind;
3956:           }
3957:         }
3958:       } else {
3959:         for (k = fdof/fcomp-1; k >= 0; --k) {
3960:           for (c = 0; c < fcomp; ++c) {
3961:             if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
3962:               indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
3963:               ++cind;
3964:             } else {
3965:               indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
3966:             }
3967:           }
3968:         }
3969:       }
3970:     }
3971:     foff     += (setBC ? fdof : (fdof - cfdof));
3972:     foffs[f] += fdof;
3973:   }
3974:   return(0);
3975: }

3979: PetscErrorCode DMPlexAnchorsModifyMat(DM dm, PetscSection section, PetscInt numPoints, PetscInt numIndices, const PetscInt points[], const PetscScalar values[], PetscInt *outNumPoints, PetscInt *outNumIndices, PetscInt *outPoints[], PetscScalar *outValues[], PetscInt offsets[], PetscBool multiplyLeft)
3980: {
3981:   Mat             cMat;
3982:   PetscSection    aSec, cSec;
3983:   IS              aIS;
3984:   PetscInt        aStart = -1, aEnd = -1;
3985:   const PetscInt  *anchors;
3986:   PetscInt        numFields, f, p, q, newP = 0;
3987:   PetscInt        newNumPoints = 0, newNumIndices = 0;
3988:   PetscInt        *newPoints, *indices, *newIndices;
3989:   PetscInt        maxAnchor, maxDof;
3990:   PetscInt        newOffsets[32];
3991:   PetscInt        *pointMatOffsets[32];
3992:   PetscInt        *newPointOffsets[32];
3993:   PetscScalar     *pointMat[32];
3994:   PetscScalar     *newValues=NULL,*tmpValues;
3995:   PetscBool       anyConstrained = PETSC_FALSE;
3996:   PetscErrorCode  ierr;

4001:   PetscSectionGetNumFields(section, &numFields);

4003:   DMPlexGetAnchors(dm,&aSec,&aIS);
4004:   /* if there are point-to-point constraints */
4005:   if (aSec) {
4006:     PetscMemzero(newOffsets, 32 * sizeof(PetscInt));
4007:     ISGetIndices(aIS,&anchors);
4008:     PetscSectionGetChart(aSec,&aStart,&aEnd);
4009:     /* figure out how many points are going to be in the new element matrix
4010:      * (we allow double counting, because it's all just going to be summed
4011:      * into the global matrix anyway) */
4012:     for (p = 0; p < 2*numPoints; p+=2) {
4013:       PetscInt b    = points[p];
4014:       PetscInt bDof = 0, bSecDof;

4016:       PetscSectionGetDof(section,b,&bSecDof);
4017:       if (!bSecDof) {
4018:         continue;
4019:       }
4020:       if (b >= aStart && b < aEnd) {
4021:         PetscSectionGetDof(aSec,b,&bDof);
4022:       }
4023:       if (bDof) {
4024:         /* this point is constrained */
4025:         /* it is going to be replaced by its anchors */
4026:         PetscInt bOff, q;

4028:         anyConstrained = PETSC_TRUE;
4029:         newNumPoints  += bDof;
4030:         PetscSectionGetOffset(aSec,b,&bOff);
4031:         for (q = 0; q < bDof; q++) {
4032:           PetscInt a = anchors[bOff + q];
4033:           PetscInt aDof;

4035:           PetscSectionGetDof(section,a,&aDof);
4036:           newNumIndices += aDof;
4037:           for (f = 0; f < numFields; ++f) {
4038:             PetscInt fDof;

4040:             PetscSectionGetFieldDof(section, a, f, &fDof);
4041:             newOffsets[f+1] += fDof;
4042:           }
4043:         }
4044:       }
4045:       else {
4046:         /* this point is not constrained */
4047:         newNumPoints++;
4048:         newNumIndices += bSecDof;
4049:         for (f = 0; f < numFields; ++f) {
4050:           PetscInt fDof;

4052:           PetscSectionGetFieldDof(section, b, f, &fDof);
4053:           newOffsets[f+1] += fDof;
4054:         }
4055:       }
4056:     }
4057:   }
4058:   if (!anyConstrained) {
4059:     if (outNumPoints)  *outNumPoints  = 0;
4060:     if (outNumIndices) *outNumIndices = 0;
4061:     if (outPoints)     *outPoints     = NULL;
4062:     if (outValues)     *outValues     = NULL;
4063:     if (aSec) {ISRestoreIndices(aIS,&anchors);}
4064:     return(0);
4065:   }

4067:   if (outNumPoints)  *outNumPoints  = newNumPoints;
4068:   if (outNumIndices) *outNumIndices = newNumIndices;

4070:   for (f = 1; f < numFields; ++f) newOffsets[f+1] += newOffsets[f];

4072:   if (!outPoints && !outValues) {
4073:     if (offsets) {
4074:       for (f = 0; f <= numFields; f++) {
4075:         offsets[f] = newOffsets[f];
4076:       }
4077:     }
4078:     if (aSec) {ISRestoreIndices(aIS,&anchors);}
4079:     return(0);
4080:   }

4082:   if (numFields && newOffsets[numFields] != newNumIndices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid size for closure %D should be %D", newOffsets[numFields], newNumIndices);

4084:   DMGetDefaultConstraints(dm, &cSec, &cMat);

4086:   /* workspaces */
4087:   if (numFields) {
4088:     for (f = 0; f < numFields; f++) {
4089:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4090:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4091:     }
4092:   }
4093:   else {
4094:     DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4095:     DMGetWorkArray(dm,numPoints,PETSC_INT,&newPointOffsets[0]);
4096:   }

4098:   /* get workspaces for the point-to-point matrices */
4099:   if (numFields) {
4100:     PetscInt totalOffset, totalMatOffset;

4102:     for (p = 0; p < numPoints; p++) {
4103:       PetscInt b    = points[2*p];
4104:       PetscInt bDof = 0, bSecDof;

4106:       PetscSectionGetDof(section,b,&bSecDof);
4107:       if (!bSecDof) {
4108:         for (f = 0; f < numFields; f++) {
4109:           newPointOffsets[f][p + 1] = 0;
4110:           pointMatOffsets[f][p + 1] = 0;
4111:         }
4112:         continue;
4113:       }
4114:       if (b >= aStart && b < aEnd) {
4115:         PetscSectionGetDof(aSec, b, &bDof);
4116:       }
4117:       if (bDof) {
4118:         for (f = 0; f < numFields; f++) {
4119:           PetscInt fDof, q, bOff, allFDof = 0;

4121:           PetscSectionGetFieldDof(section, b, f, &fDof);
4122:           PetscSectionGetOffset(aSec, b, &bOff);
4123:           for (q = 0; q < bDof; q++) {
4124:             PetscInt a = anchors[bOff + q];
4125:             PetscInt aFDof;

4127:             PetscSectionGetFieldDof(section, a, f, &aFDof);
4128:             allFDof += aFDof;
4129:           }
4130:           newPointOffsets[f][p+1] = allFDof;
4131:           pointMatOffsets[f][p+1] = fDof * allFDof;
4132:         }
4133:       }
4134:       else {
4135:         for (f = 0; f < numFields; f++) {
4136:           PetscInt fDof;

4138:           PetscSectionGetFieldDof(section, b, f, &fDof);
4139:           newPointOffsets[f][p+1] = fDof;
4140:           pointMatOffsets[f][p+1] = 0;
4141:         }
4142:       }
4143:     }
4144:     for (f = 0, totalOffset = 0, totalMatOffset = 0; f < numFields; f++) {
4145:       newPointOffsets[f][0] = totalOffset;
4146:       pointMatOffsets[f][0] = totalMatOffset;
4147:       for (p = 0; p < numPoints; p++) {
4148:         newPointOffsets[f][p+1] += newPointOffsets[f][p];
4149:         pointMatOffsets[f][p+1] += pointMatOffsets[f][p];
4150:       }
4151:       totalOffset    = newPointOffsets[f][numPoints];
4152:       totalMatOffset = pointMatOffsets[f][numPoints];
4153:       DMGetWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4154:     }
4155:   }
4156:   else {
4157:     for (p = 0; p < numPoints; p++) {
4158:       PetscInt b    = points[2*p];
4159:       PetscInt bDof = 0, bSecDof;

4161:       PetscSectionGetDof(section,b,&bSecDof);
4162:       if (!bSecDof) {
4163:         newPointOffsets[0][p + 1] = 0;
4164:         pointMatOffsets[0][p + 1] = 0;
4165:         continue;
4166:       }
4167:       if (b >= aStart && b < aEnd) {
4168:         PetscSectionGetDof(aSec, b, &bDof);
4169:       }
4170:       if (bDof) {
4171:         PetscInt bOff, q, allDof = 0;

4173:         PetscSectionGetOffset(aSec, b, &bOff);
4174:         for (q = 0; q < bDof; q++) {
4175:           PetscInt a = anchors[bOff + q], aDof;

4177:           PetscSectionGetDof(section, a, &aDof);
4178:           allDof += aDof;
4179:         }
4180:         newPointOffsets[0][p+1] = allDof;
4181:         pointMatOffsets[0][p+1] = bSecDof * allDof;
4182:       }
4183:       else {
4184:         newPointOffsets[0][p+1] = bSecDof;
4185:         pointMatOffsets[0][p+1] = 0;
4186:       }
4187:     }
4188:     newPointOffsets[0][0] = 0;
4189:     pointMatOffsets[0][0] = 0;
4190:     for (p = 0; p < numPoints; p++) {
4191:       newPointOffsets[0][p+1] += newPointOffsets[0][p];
4192:       pointMatOffsets[0][p+1] += pointMatOffsets[0][p];
4193:     }
4194:     DMGetWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4195:   }

4197:   /* output arrays */
4198:   DMGetWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);

4200:   /* get the point-to-point matrices; construct newPoints */
4201:   PetscSectionGetMaxDof(aSec, &maxAnchor);
4202:   PetscSectionGetMaxDof(section, &maxDof);
4203:   DMGetWorkArray(dm,maxDof,PETSC_INT,&indices);
4204:   DMGetWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4205:   if (numFields) {
4206:     for (p = 0, newP = 0; p < numPoints; p++) {
4207:       PetscInt b    = points[2*p];
4208:       PetscInt o    = points[2*p+1];
4209:       PetscInt bDof = 0, bSecDof;

4211:       PetscSectionGetDof(section, b, &bSecDof);
4212:       if (!bSecDof) {
4213:         continue;
4214:       }
4215:       if (b >= aStart && b < aEnd) {
4216:         PetscSectionGetDof(aSec, b, &bDof);
4217:       }
4218:       if (bDof) {
4219:         PetscInt fStart[32], fEnd[32], fAnchorStart[32], fAnchorEnd[32], bOff, q;

4221:         fStart[0] = 0;
4222:         fEnd[0]   = 0;
4223:         for (f = 0; f < numFields; f++) {
4224:           PetscInt fDof;

4226:           PetscSectionGetFieldDof(cSec, b, f, &fDof);
4227:           fStart[f+1] = fStart[f] + fDof;
4228:           fEnd[f+1]   = fStart[f+1];
4229:         }
4230:         PetscSectionGetOffset(cSec, b, &bOff);
4231:         indicesPointFields_private(cSec, b, bOff, fEnd, PETSC_TRUE, o, indices);

4233:         fAnchorStart[0] = 0;
4234:         fAnchorEnd[0]   = 0;
4235:         for (f = 0; f < numFields; f++) {
4236:           PetscInt fDof = newPointOffsets[f][p + 1] - newPointOffsets[f][p];

4238:           fAnchorStart[f+1] = fAnchorStart[f] + fDof;
4239:           fAnchorEnd[f+1]   = fAnchorStart[f + 1];
4240:         }
4241:         PetscSectionGetOffset(aSec, b, &bOff);
4242:         for (q = 0; q < bDof; q++) {
4243:           PetscInt a = anchors[bOff + q], aOff;

4245:           /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4246:           newPoints[2*(newP + q)]     = a;
4247:           newPoints[2*(newP + q) + 1] = 0;
4248:           PetscSectionGetOffset(section, a, &aOff);
4249:           indicesPointFields_private(section, a, aOff, fAnchorEnd, PETSC_TRUE, 0, newIndices);
4250:         }
4251:         newP += bDof;

4253:         if (outValues) {
4254:           /* get the point-to-point submatrix */
4255:           for (f = 0; f < numFields; f++) {
4256:             MatGetValues(cMat,fEnd[f]-fStart[f],indices + fStart[f],fAnchorEnd[f] - fAnchorStart[f],newIndices + fAnchorStart[f],pointMat[f] + pointMatOffsets[f][p]);
4257:           }
4258:         }
4259:       }
4260:       else {
4261:         newPoints[2 * newP]     = b;
4262:         newPoints[2 * newP + 1] = o;
4263:         newP++;
4264:       }
4265:     }
4266:   } else {
4267:     for (p = 0; p < numPoints; p++) {
4268:       PetscInt b    = points[2*p];
4269:       PetscInt o    = points[2*p+1];
4270:       PetscInt bDof = 0, bSecDof;

4272:       PetscSectionGetDof(section, b, &bSecDof);
4273:       if (!bSecDof) {
4274:         continue;
4275:       }
4276:       if (b >= aStart && b < aEnd) {
4277:         PetscSectionGetDof(aSec, b, &bDof);
4278:       }
4279:       if (bDof) {
4280:         PetscInt bEnd = 0, bAnchorEnd = 0, bOff;

4282:         PetscSectionGetOffset(cSec, b, &bOff);
4283:         indicesPoint_private(cSec, b, bOff, &bEnd, PETSC_TRUE, o, indices);

4285:         PetscSectionGetOffset (aSec, b, &bOff);
4286:         for (q = 0; q < bDof; q++) {
4287:           PetscInt a = anchors[bOff + q], aOff;

4289:           /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */

4291:           newPoints[2*(newP + q)]     = a;
4292:           newPoints[2*(newP + q) + 1] = 0;
4293:           PetscSectionGetOffset(section, a, &aOff);
4294:           indicesPoint_private(section, a, aOff, &bAnchorEnd, PETSC_TRUE, 0, newIndices);
4295:         }
4296:         newP += bDof;

4298:         /* get the point-to-point submatrix */
4299:         if (outValues) {
4300:           MatGetValues(cMat,bEnd,indices,bAnchorEnd,newIndices,pointMat[0] + pointMatOffsets[0][p]);
4301:         }
4302:       }
4303:       else {
4304:         newPoints[2 * newP]     = b;
4305:         newPoints[2 * newP + 1] = o;
4306:         newP++;
4307:       }
4308:     }
4309:   }

4311:   if (outValues) {
4312:     DMGetWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4313:     PetscMemzero(tmpValues,newNumIndices*numIndices*sizeof(*tmpValues));
4314:     /* multiply constraints on the right */
4315:     if (numFields) {
4316:       for (f = 0; f < numFields; f++) {
4317:         PetscInt oldOff = offsets[f];

4319:         for (p = 0; p < numPoints; p++) {
4320:           PetscInt cStart = newPointOffsets[f][p];
4321:           PetscInt b      = points[2 * p];
4322:           PetscInt c, r, k;
4323:           PetscInt dof;

4325:           PetscSectionGetFieldDof(section,b,f,&dof);
4326:           if (!dof) {
4327:             continue;
4328:           }
4329:           if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4330:             PetscInt nCols         = newPointOffsets[f][p+1]-cStart;
4331:             const PetscScalar *mat = pointMat[f] + pointMatOffsets[f][p];

4333:             for (r = 0; r < numIndices; r++) {
4334:               for (c = 0; c < nCols; c++) {
4335:                 for (k = 0; k < dof; k++) {
4336:                   tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4337:                 }
4338:               }
4339:             }
4340:           }
4341:           else {
4342:             /* copy this column as is */
4343:             for (r = 0; r < numIndices; r++) {
4344:               for (c = 0; c < dof; c++) {
4345:                 tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4346:               }
4347:             }
4348:           }
4349:           oldOff += dof;
4350:         }
4351:       }
4352:     }
4353:     else {
4354:       PetscInt oldOff = 0;
4355:       for (p = 0; p < numPoints; p++) {
4356:         PetscInt cStart = newPointOffsets[0][p];
4357:         PetscInt b      = points[2 * p];
4358:         PetscInt c, r, k;
4359:         PetscInt dof;

4361:         PetscSectionGetDof(section,b,&dof);
4362:         if (!dof) {
4363:           continue;
4364:         }
4365:         if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4366:           PetscInt nCols         = newPointOffsets[0][p+1]-cStart;
4367:           const PetscScalar *mat = pointMat[0] + pointMatOffsets[0][p];

4369:           for (r = 0; r < numIndices; r++) {
4370:             for (c = 0; c < nCols; c++) {
4371:               for (k = 0; k < dof; k++) {
4372:                 tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4373:               }
4374:             }
4375:           }
4376:         }
4377:         else {
4378:           /* copy this column as is */
4379:           for (r = 0; r < numIndices; r++) {
4380:             for (c = 0; c < dof; c++) {
4381:               tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4382:             }
4383:           }
4384:         }
4385:         oldOff += dof;
4386:       }
4387:     }

4389:     if (multiplyLeft) {
4390:       DMGetWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4391:       PetscMemzero(newValues,newNumIndices*newNumIndices*sizeof(*newValues));
4392:       /* multiply constraints transpose on the left */
4393:       if (numFields) {
4394:         for (f = 0; f < numFields; f++) {
4395:           PetscInt oldOff = offsets[f];

4397:           for (p = 0; p < numPoints; p++) {
4398:             PetscInt rStart = newPointOffsets[f][p];
4399:             PetscInt b      = points[2 * p];
4400:             PetscInt c, r, k;
4401:             PetscInt dof;

4403:             PetscSectionGetFieldDof(section,b,f,&dof);
4404:             if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4405:               PetscInt nRows                        = newPointOffsets[f][p+1]-rStart;
4406:               const PetscScalar *PETSC_RESTRICT mat = pointMat[f] + pointMatOffsets[f][p];

4408:               for (r = 0; r < nRows; r++) {
4409:                 for (c = 0; c < newNumIndices; c++) {
4410:                   for (k = 0; k < dof; k++) {
4411:                     newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4412:                   }
4413:                 }
4414:               }
4415:             }
4416:             else {
4417:               /* copy this row as is */
4418:               for (r = 0; r < dof; r++) {
4419:                 for (c = 0; c < newNumIndices; c++) {
4420:                   newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4421:                 }
4422:               }
4423:             }
4424:             oldOff += dof;
4425:           }
4426:         }
4427:       }
4428:       else {
4429:         PetscInt oldOff = 0;

4431:         for (p = 0; p < numPoints; p++) {
4432:           PetscInt rStart = newPointOffsets[0][p];
4433:           PetscInt b      = points[2 * p];
4434:           PetscInt c, r, k;
4435:           PetscInt dof;

4437:           PetscSectionGetDof(section,b,&dof);
4438:           if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4439:             PetscInt nRows                        = newPointOffsets[0][p+1]-rStart;
4440:             const PetscScalar *PETSC_RESTRICT mat = pointMat[0] + pointMatOffsets[0][p];

4442:             for (r = 0; r < nRows; r++) {
4443:               for (c = 0; c < newNumIndices; c++) {
4444:                 for (k = 0; k < dof; k++) {
4445:                   newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4446:                 }
4447:               }
4448:             }
4449:           }
4450:           else {
4451:             /* copy this row as is */
4452:             for (r = 0; r < dof; r++) {
4453:               for (c = 0; c < newNumIndices; c++) {
4454:                 newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4455:               }
4456:             }
4457:           }
4458:           oldOff += dof;
4459:         }
4460:       }

4462:       DMRestoreWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4463:     }
4464:     else {
4465:       newValues = tmpValues;
4466:     }
4467:   }

4469:   /* clean up */
4470:   DMRestoreWorkArray(dm,maxDof,PETSC_INT,&indices);
4471:   DMRestoreWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);

4473:   if (numFields) {
4474:     for (f = 0; f < numFields; f++) {
4475:       DMRestoreWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4476:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4477:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4478:     }
4479:   }
4480:   else {
4481:     DMRestoreWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4482:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4483:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[0]);
4484:   }
4485:   ISRestoreIndices(aIS,&anchors);

4487:   /* output */
4488:   if (outPoints) {
4489:     *outPoints = newPoints;
4490:   }
4491:   else {
4492:     DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4493:   }
4494:   if (outValues) {
4495:     *outValues = newValues;
4496:   }
4497:   else {
4498:     DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4499:   }
4500:   for (f = 0; f <= numFields; f++) {
4501:     offsets[f] = newOffsets[f];
4502:   }
4503:   return(0);
4504: }

4508: PetscErrorCode DMPlexGetClosureIndices(DM dm, PetscSection section, PetscSection globalSection, PetscInt point, PetscInt *numIndices, PetscInt **indices, PetscInt *outOffsets)
4509: {
4510:   PetscSection    clSection;
4511:   IS              clPoints;
4512:   const PetscInt *clp;
4513:   PetscInt       *points = NULL, *pointsNew;
4514:   PetscInt        numPoints, numPointsNew;
4515:   PetscInt        offsets[32];
4516:   PetscInt        Nf, Nind, NindNew, off, globalOff, f, p;
4517:   PetscErrorCode  ierr;

4525:   PetscSectionGetNumFields(section, &Nf);
4526:   if (Nf > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", Nf);
4527:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
4528:   /* Get points in closure */
4529:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4530:   if (!clPoints) {
4531:     PetscInt pStart, pEnd, q;

4533:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4534:     /* Compress out points not in the section */
4535:     PetscSectionGetChart(section, &pStart, &pEnd);
4536:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4537:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4538:         points[q*2]   = points[p];
4539:         points[q*2+1] = points[p+1];
4540:         ++q;
4541:       }
4542:     }
4543:     numPoints = q;
4544:   } else {
4545:     PetscInt dof, off;

4547:     PetscSectionGetDof(clSection, point, &dof);
4548:     numPoints = dof/2;
4549:     PetscSectionGetOffset(clSection, point, &off);
4550:     ISGetIndices(clPoints, &clp);
4551:     points = (PetscInt *) &clp[off];
4552:   }
4553:   /* Get number of indices and indices per field */
4554:   for (p = 0, Nind = 0; p < numPoints*2; p += 2) {
4555:     PetscInt dof, fdof;

4557:     PetscSectionGetDof(section, points[p], &dof);
4558:     for (f = 0; f < Nf; ++f) {
4559:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
4560:       offsets[f+1] += fdof;
4561:     }
4562:     Nind += dof;
4563:   }
4564:   for (f = 1; f < Nf; ++f) offsets[f+1] += offsets[f];
4565:   if (Nf && offsets[Nf] != Nind) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[Nf], Nind);
4566:   /* Correct for hanging node constraints */
4567:   {
4568:     DMPlexAnchorsModifyMat(dm, section, numPoints, Nind, points, NULL, &numPointsNew, &NindNew, &pointsNew, NULL, offsets, PETSC_TRUE);
4569:     if (numPointsNew) {
4570:       if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4571:       else           {ISRestoreIndices(clPoints, &clp);}
4572:       numPoints = numPointsNew;
4573:       Nind      = NindNew;
4574:       points    = pointsNew;
4575:     }
4576:   }
4577:   /* Calculate indices */
4578:   DMGetWorkArray(dm, Nind, PETSC_INT, indices);
4579:   if (Nf) {
4580:     if (outOffsets) {
4581:       PetscInt f;

4583:       for (f = 0; f <= Nf; f++) {
4584:         outOffsets[f] = offsets[f];
4585:       }
4586:     }
4587:     for (p = 0; p < numPoints*2; p += 2) {
4588:       PetscInt o = points[p+1];
4589:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4590:       indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, *indices);
4591:     }
4592:   } else {
4593:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
4594:       PetscInt o = points[p+1];
4595:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4596:       indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, *indices);
4597:     }
4598:   }
4599:   /* Cleanup points */
4600:   if (numPointsNew) {
4601:     DMRestoreWorkArray(dm, 2*numPointsNew, PETSC_INT, &pointsNew);
4602:   } else {
4603:     if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4604:     else           {ISRestoreIndices(clPoints, &clp);}
4605:   }
4606:   if (numIndices) *numIndices = Nind;
4607:   return(0);
4608: }

4612: PetscErrorCode DMPlexRestoreClosureIndices(DM dm, PetscSection section, PetscSection globalSection, PetscInt point, PetscInt *numIndices, PetscInt **indices,PetscInt *outOffsets)
4613: {

4619:   DMRestoreWorkArray(dm, 0, PETSC_INT, indices);
4620:   return(0);
4621: }

4625: /*@C
4626:   DMPlexMatSetClosure - Set an array of the values on the closure of 'point'

4628:   Not collective

4630:   Input Parameters:
4631: + dm - The DM
4632: . section - The section describing the layout in v, or NULL to use the default section
4633: . globalSection - The section describing the layout in v, or NULL to use the default global section
4634: . A - The matrix
4635: . point - The sieve point in the DM
4636: . values - The array of values
4637: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

4639:   Fortran Notes:
4640:   This routine is only available in Fortran 90, and you must include petsc.h90 in your code.

4642:   Level: intermediate

4644: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
4645: @*/
4646: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4647: {
4648:   DM_Plex        *mesh   = (DM_Plex*) dm->data;
4649:   PetscSection    clSection;
4650:   IS              clPoints;
4651:   PetscInt       *points = NULL, *newPoints;
4652:   const PetscInt *clp;
4653:   PetscInt       *indices;
4654:   PetscInt        offsets[32];
4655:   PetscInt        numFields, numPoints, newNumPoints, numIndices, newNumIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
4656:   PetscScalar    *newValues;
4657:   PetscErrorCode  ierr;

4661:   if (!section) {DMGetDefaultSection(dm, &section);}
4663:   if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
4666:   PetscSectionGetNumFields(section, &numFields);
4667:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4668:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
4669:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4670:   if (!clPoints) {
4671:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4672:     /* Compress out points not in the section */
4673:     PetscSectionGetChart(section, &pStart, &pEnd);
4674:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4675:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4676:         points[q*2]   = points[p];
4677:         points[q*2+1] = points[p+1];
4678:         ++q;
4679:       }
4680:     }
4681:     numPoints = q;
4682:   } else {
4683:     PetscInt dof, off;

4685:     PetscSectionGetDof(clSection, point, &dof);
4686:     numPoints = dof/2;
4687:     PetscSectionGetOffset(clSection, point, &off);
4688:     ISGetIndices(clPoints, &clp);
4689:     points = (PetscInt *) &clp[off];
4690:   }
4691:   for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
4692:     PetscInt fdof;

4694:     PetscSectionGetDof(section, points[p], &dof);
4695:     for (f = 0; f < numFields; ++f) {
4696:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
4697:       offsets[f+1] += fdof;
4698:     }
4699:     numIndices += dof;
4700:   }
4701:   for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];

4703:   if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
4704:   DMPlexAnchorsModifyMat(dm,section,numPoints,numIndices,points,values,&newNumPoints,&newNumIndices,&newPoints,&newValues,offsets,PETSC_TRUE);
4705:   if (newNumPoints) {
4706:     if (!clPoints) {
4707:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4708:     } else {
4709:       ISRestoreIndices(clPoints, &clp);
4710:     }
4711:     numPoints  = newNumPoints;
4712:     numIndices = newNumIndices;
4713:     points     = newPoints;
4714:     values     = newValues;
4715:   }
4716:   DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
4717:   if (numFields) {
4718:     for (p = 0; p < numPoints*2; p += 2) {
4719:       PetscInt o = points[p+1];
4720:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4721:       indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
4722:     }
4723:   } else {
4724:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
4725:       PetscInt o = points[p+1];
4726:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4727:       indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
4728:     }
4729:   }
4730:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
4731:   MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
4732:   if (mesh->printFEM > 1) {
4733:     PetscInt i;
4734:     PetscPrintf(PETSC_COMM_SELF, "  Indices:");
4735:     for (i = 0; i < numIndices; ++i) {PetscPrintf(PETSC_COMM_SELF, " %d", indices[i]);}
4736:     PetscPrintf(PETSC_COMM_SELF, "\n");
4737:   }
4738:   if (ierr) {
4739:     PetscMPIInt    rank;
4740:     PetscErrorCode ierr2;

4742:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4743:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4744:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
4745:     ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
4746: 
4747:   }
4748:   if (newNumPoints) {
4749:     DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4750:     DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4751:   }
4752:   else {
4753:     if (!clPoints) {
4754:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4755:     } else {
4756:       ISRestoreIndices(clPoints, &clp);
4757:     }
4758:   }
4759:   DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
4760:   return(0);
4761: }

4765: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4766: {
4767:   DM_Plex        *mesh   = (DM_Plex*) dmf->data;
4768:   PetscInt       *fpoints = NULL, *ftotpoints = NULL;
4769:   PetscInt       *cpoints = NULL;
4770:   PetscInt       *findices, *cindices;
4771:   PetscInt        foffsets[32], coffsets[32];
4772:   CellRefiner     cellRefiner;
4773:   PetscInt        numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4774:   PetscErrorCode  ierr;

4779:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4781:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
4783:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4785:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4788:   PetscSectionGetNumFields(fsection, &numFields);
4789:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4790:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4791:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4792:   /* Column indices */
4793:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4794:   maxFPoints = numCPoints;
4795:   /* Compress out points not in the section */
4796:   /*   TODO: Squeeze out points with 0 dof as well */
4797:   PetscSectionGetChart(csection, &pStart, &pEnd);
4798:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4799:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4800:       cpoints[q*2]   = cpoints[p];
4801:       cpoints[q*2+1] = cpoints[p+1];
4802:       ++q;
4803:     }
4804:   }
4805:   numCPoints = q;
4806:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4807:     PetscInt fdof;

4809:     PetscSectionGetDof(csection, cpoints[p], &dof);
4810:     if (!dof) continue;
4811:     for (f = 0; f < numFields; ++f) {
4812:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4813:       coffsets[f+1] += fdof;
4814:     }
4815:     numCIndices += dof;
4816:   }
4817:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4818:   /* Row indices */
4819:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4820:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4821:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4822:   for (r = 0, q = 0; r < numSubcells; ++r) {
4823:     /* TODO Map from coarse to fine cells */
4824:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4825:     /* Compress out points not in the section */
4826:     PetscSectionGetChart(fsection, &pStart, &pEnd);
4827:     for (p = 0; p < numFPoints*2; p += 2) {
4828:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4829:         PetscSectionGetDof(fsection, fpoints[p], &dof);
4830:         if (!dof) continue;
4831:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4832:         if (s < q) continue;
4833:         ftotpoints[q*2]   = fpoints[p];
4834:         ftotpoints[q*2+1] = fpoints[p+1];
4835:         ++q;
4836:       }
4837:     }
4838:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4839:   }
4840:   numFPoints = q;
4841:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4842:     PetscInt fdof;

4844:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4845:     if (!dof) continue;
4846:     for (f = 0; f < numFields; ++f) {
4847:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4848:       foffsets[f+1] += fdof;
4849:     }
4850:     numFIndices += dof;
4851:   }
4852:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

4854:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4855:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4856:   DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4857:   DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4858:   if (numFields) {
4859:     for (p = 0; p < numFPoints*2; p += 2) {
4860:       PetscInt o = ftotpoints[p+1];
4861:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4862:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4863:     }
4864:     for (p = 0; p < numCPoints*2; p += 2) {
4865:       PetscInt o = cpoints[p+1];
4866:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4867:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4868:     }
4869:   } else {
4870:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4871:       PetscInt o = ftotpoints[p+1];
4872:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4873:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4874:     }
4875:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4876:       PetscInt o = cpoints[p+1];
4877:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4878:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4879:     }
4880:   }
4881:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
4882:   MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
4883:   if (ierr) {
4884:     PetscMPIInt    rank;
4885:     PetscErrorCode ierr2;

4887:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4888:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4889:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
4890:     ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
4891:     ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
4892: 
4893:   }
4894:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4895:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4896:   DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4897:   DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4898:   return(0);
4899: }

4903: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
4904: {
4905:   PetscInt      *fpoints = NULL, *ftotpoints = NULL;
4906:   PetscInt      *cpoints = NULL;
4907:   PetscInt       foffsets[32], coffsets[32];
4908:   CellRefiner    cellRefiner;
4909:   PetscInt       numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;

4915:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4917:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
4919:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4921:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4923:   PetscSectionGetNumFields(fsection, &numFields);
4924:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4925:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4926:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4927:   /* Column indices */
4928:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4929:   maxFPoints = numCPoints;
4930:   /* Compress out points not in the section */
4931:   /*   TODO: Squeeze out points with 0 dof as well */
4932:   PetscSectionGetChart(csection, &pStart, &pEnd);
4933:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4934:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4935:       cpoints[q*2]   = cpoints[p];
4936:       cpoints[q*2+1] = cpoints[p+1];
4937:       ++q;
4938:     }
4939:   }
4940:   numCPoints = q;
4941:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4942:     PetscInt fdof;

4944:     PetscSectionGetDof(csection, cpoints[p], &dof);
4945:     if (!dof) continue;
4946:     for (f = 0; f < numFields; ++f) {
4947:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4948:       coffsets[f+1] += fdof;
4949:     }
4950:     numCIndices += dof;
4951:   }
4952:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4953:   /* Row indices */
4954:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4955:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4956:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4957:   for (r = 0, q = 0; r < numSubcells; ++r) {
4958:     /* TODO Map from coarse to fine cells */
4959:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4960:     /* Compress out points not in the section */
4961:     PetscSectionGetChart(fsection, &pStart, &pEnd);
4962:     for (p = 0; p < numFPoints*2; p += 2) {
4963:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4964:         PetscSectionGetDof(fsection, fpoints[p], &dof);
4965:         if (!dof) continue;
4966:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4967:         if (s < q) continue;
4968:         ftotpoints[q*2]   = fpoints[p];
4969:         ftotpoints[q*2+1] = fpoints[p+1];
4970:         ++q;
4971:       }
4972:     }
4973:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4974:   }
4975:   numFPoints = q;
4976:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4977:     PetscInt fdof;

4979:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4980:     if (!dof) continue;
4981:     for (f = 0; f < numFields; ++f) {
4982:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4983:       foffsets[f+1] += fdof;
4984:     }
4985:     numFIndices += dof;
4986:   }
4987:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

4989:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4990:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4991:   if (numFields) {
4992:     for (p = 0; p < numFPoints*2; p += 2) {
4993:       PetscInt o = ftotpoints[p+1];
4994:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4995:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4996:     }
4997:     for (p = 0; p < numCPoints*2; p += 2) {
4998:       PetscInt o = cpoints[p+1];
4999:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5000:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5001:     }
5002:   } else {
5003:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5004:       PetscInt o = ftotpoints[p+1];
5005:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5006:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5007:     }
5008:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5009:       PetscInt o = cpoints[p+1];
5010:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5011:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5012:     }
5013:   }
5014:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5015:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5016:   return(0);
5017: }

5021: /*@
5022:   DMPlexGetHybridBounds - Get the first mesh point of each dimension which is a hybrid

5024:   Input Parameter:
5025: . dm - The DMPlex object

5027:   Output Parameters:
5028: + cMax - The first hybrid cell
5029: . fMax - The first hybrid face
5030: . eMax - The first hybrid edge
5031: - vMax - The first hybrid vertex

5033:   Level: developer

5035: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
5036: @*/
5037: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
5038: {
5039:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5040:   PetscInt       dim;

5045:   DMGetDimension(dm, &dim);
5046:   if (cMax) *cMax = mesh->hybridPointMax[dim];
5047:   if (fMax) *fMax = mesh->hybridPointMax[dim-1];
5048:   if (eMax) *eMax = mesh->hybridPointMax[1];
5049:   if (vMax) *vMax = mesh->hybridPointMax[0];
5050:   return(0);
5051: }

5055: /*@
5056:   DMPlexSetHybridBounds - Set the first mesh point of each dimension which is a hybrid

5058:   Input Parameters:
5059: . dm   - The DMPlex object
5060: . cMax - The first hybrid cell
5061: . fMax - The first hybrid face
5062: . eMax - The first hybrid edge
5063: - vMax - The first hybrid vertex

5065:   Level: developer

5067: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5068: @*/
5069: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5070: {
5071:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5072:   PetscInt       dim;

5077:   DMGetDimension(dm, &dim);
5078:   if (cMax >= 0) mesh->hybridPointMax[dim]   = cMax;
5079:   if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5080:   if (eMax >= 0) mesh->hybridPointMax[1]     = eMax;
5081:   if (vMax >= 0) mesh->hybridPointMax[0]     = vMax;
5082:   return(0);
5083: }

5087: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5088: {
5089:   DM_Plex *mesh = (DM_Plex*) dm->data;

5094:   *cellHeight = mesh->vtkCellHeight;
5095:   return(0);
5096: }

5100: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5101: {
5102:   DM_Plex *mesh = (DM_Plex*) dm->data;

5106:   mesh->vtkCellHeight = cellHeight;
5107:   return(0);
5108: }

5112: /* We can easily have a form that takes an IS instead */
5113: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5114: {
5115:   PetscSection   section, globalSection;
5116:   PetscInt      *numbers, p;

5120:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
5121:   PetscSectionSetChart(section, pStart, pEnd);
5122:   for (p = pStart; p < pEnd; ++p) {
5123:     PetscSectionSetDof(section, p, 1);
5124:   }
5125:   PetscSectionSetUp(section);
5126:   PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
5127:   PetscMalloc1(pEnd - pStart, &numbers);
5128:   for (p = pStart; p < pEnd; ++p) {
5129:     PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5130:     if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5131:     else                       numbers[p-pStart] += shift;
5132:   }
5133:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5134:   if (globalSize) {
5135:     PetscLayout layout;
5136:     PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5137:     PetscLayoutGetSize(layout, globalSize);
5138:     PetscLayoutDestroy(&layout);
5139:   }
5140:   PetscSectionDestroy(&section);
5141:   PetscSectionDestroy(&globalSection);
5142:   return(0);
5143: }

5147: PetscErrorCode DMPlexCreateCellNumbering_Internal(DM dm, PetscBool includeHybrid, IS *globalCellNumbers)
5148: {
5149:   PetscInt       cellHeight, cStart, cEnd, cMax;

5153:   DMPlexGetVTKCellHeight(dm, &cellHeight);
5154:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5155:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5156:   if (cMax >= 0 && !includeHybrid) cEnd = PetscMin(cEnd, cMax);
5157:   DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, globalCellNumbers);
5158:   return(0);
5159: }

5163: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5164: {
5165:   DM_Plex       *mesh = (DM_Plex*) dm->data;

5170:   if (!mesh->globalCellNumbers) {DMPlexCreateCellNumbering_Internal(dm, PETSC_FALSE, &mesh->globalCellNumbers);}
5171:   *globalCellNumbers = mesh->globalCellNumbers;
5172:   return(0);
5173: }

5177: PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM dm, PetscBool includeHybrid, IS *globalVertexNumbers)
5178: {
5179:   PetscInt       vStart, vEnd, vMax;

5184:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5185:   DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5186:   if (vMax >= 0 && !includeHybrid) vEnd = PetscMin(vEnd, vMax);
5187:   DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, globalVertexNumbers);
5188:   return(0);
5189: }

5193: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5194: {
5195:   DM_Plex       *mesh = (DM_Plex*) dm->data;

5200:   if (!mesh->globalVertexNumbers) {DMPlexCreateVertexNumbering_Internal(dm, PETSC_FALSE, &mesh->globalVertexNumbers);}
5201:   *globalVertexNumbers = mesh->globalVertexNumbers;
5202:   return(0);
5203: }

5207: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5208: {
5209:   IS             nums[4];
5210:   PetscInt       depths[4];
5211:   PetscInt       depth, d, shift = 0;

5216:   DMPlexGetDepth(dm, &depth);
5217:   /* For unstratified meshes use dim instead of depth */
5218:   if (depth < 0) {DMGetDimension(dm, &depth);}
5219:   depths[0] = depth; depths[1] = 0;
5220:   for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5221:   for (d = 0; d <= depth; ++d) {
5222:     PetscInt pStart, pEnd, gsize;

5224:     DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5225:     DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5226:     shift += gsize;
5227:   }
5228:   ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5229:   for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5230:   return(0);
5231: }

5235: /*@
5236:   DMPlexCheckSymmetry - Check that the adjacency information in the mesh is symmetric.

5238:   Input Parameters:
5239:   + dm - The DMPlex object

5241:   Note: This is a useful diagnostic when creating meshes programmatically.

5243:   Level: developer

5245: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5246: @*/
5247: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5248: {
5249:   PetscSection    coneSection, supportSection;
5250:   const PetscInt *cone, *support;
5251:   PetscInt        coneSize, c, supportSize, s;
5252:   PetscInt        pStart, pEnd, p, csize, ssize;
5253:   PetscErrorCode  ierr;

5257:   DMPlexGetConeSection(dm, &coneSection);
5258:   DMPlexGetSupportSection(dm, &supportSection);
5259:   /* Check that point p is found in the support of its cone points, and vice versa */
5260:   DMPlexGetChart(dm, &pStart, &pEnd);
5261:   for (p = pStart; p < pEnd; ++p) {
5262:     DMPlexGetConeSize(dm, p, &coneSize);
5263:     DMPlexGetCone(dm, p, &cone);
5264:     for (c = 0; c < coneSize; ++c) {
5265:       PetscBool dup = PETSC_FALSE;
5266:       PetscInt  d;
5267:       for (d = c-1; d >= 0; --d) {
5268:         if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5269:       }
5270:       DMPlexGetSupportSize(dm, cone[c], &supportSize);
5271:       DMPlexGetSupport(dm, cone[c], &support);
5272:       for (s = 0; s < supportSize; ++s) {
5273:         if (support[s] == p) break;
5274:       }
5275:       if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5276:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5277:         for (s = 0; s < coneSize; ++s) {
5278:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
5279:         }
5280:         PetscPrintf(PETSC_COMM_SELF, "\n");
5281:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
5282:         for (s = 0; s < supportSize; ++s) {
5283:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
5284:         }
5285:         PetscPrintf(PETSC_COMM_SELF, "\n");
5286:         if (dup) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
5287:         else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
5288:       }
5289:     }
5290:     DMPlexGetSupportSize(dm, p, &supportSize);
5291:     DMPlexGetSupport(dm, p, &support);
5292:     for (s = 0; s < supportSize; ++s) {
5293:       DMPlexGetConeSize(dm, support[s], &coneSize);
5294:       DMPlexGetCone(dm, support[s], &cone);
5295:       for (c = 0; c < coneSize; ++c) {
5296:         if (cone[c] == p) break;
5297:       }
5298:       if (c >= coneSize) {
5299:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
5300:         for (c = 0; c < supportSize; ++c) {
5301:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
5302:         }
5303:         PetscPrintf(PETSC_COMM_SELF, "\n");
5304:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
5305:         for (c = 0; c < coneSize; ++c) {
5306:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
5307:         }
5308:         PetscPrintf(PETSC_COMM_SELF, "\n");
5309:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
5310:       }
5311:     }
5312:   }
5313:   PetscSectionGetStorageSize(coneSection, &csize);
5314:   PetscSectionGetStorageSize(supportSection, &ssize);
5315:   if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
5316:   return(0);
5317: }

5321: /*@
5322:   DMPlexCheckSkeleton - Check that each cell has the correct number of vertices

5324:   Input Parameters:
5325: + dm - The DMPlex object
5326: . isSimplex - Are the cells simplices or tensor products
5327: - cellHeight - Normally 0

5329:   Note: This is a useful diagnostic when creating meshes programmatically.

5331:   Level: developer

5333: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
5334: @*/
5335: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5336: {
5337:   PetscInt       dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;

5342:   DMGetDimension(dm, &dim);
5343:   switch (dim) {
5344:   case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
5345:   case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
5346:   case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
5347:   default:
5348:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
5349:   }
5350:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5351:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5352:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5353:   cMax = cMax >= 0 ? cMax : cEnd;
5354:   for (c = cStart; c < cMax; ++c) {
5355:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5357:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5358:     for (cl = 0; cl < closureSize*2; cl += 2) {
5359:       const PetscInt p = closure[cl];
5360:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5361:     }
5362:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5363:     if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has  %d vertices != %d", c, coneSize, numCorners);
5364:   }
5365:   for (c = cMax; c < cEnd; ++c) {
5366:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5368:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5369:     for (cl = 0; cl < closureSize*2; cl += 2) {
5370:       const PetscInt p = closure[cl];
5371:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5372:     }
5373:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5374:     if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has  %d vertices > %d", c, coneSize, numHybridCorners);
5375:   }
5376:   return(0);
5377: }

5381: /*@
5382:   DMPlexCheckFaces - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type

5384:   Input Parameters:
5385: + dm - The DMPlex object
5386: . isSimplex - Are the cells simplices or tensor products
5387: - cellHeight - Normally 0

5389:   Note: This is a useful diagnostic when creating meshes programmatically.

5391:   Level: developer

5393: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
5394: @*/
5395: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5396: {
5397:   PetscInt       pMax[4];
5398:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, h;

5403:   DMGetDimension(dm, &dim);
5404:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5405:   DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
5406:   for (h = cellHeight; h < dim; ++h) {
5407:     DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
5408:     for (c = cStart; c < cEnd; ++c) {
5409:       const PetscInt *cone, *ornt, *faces;
5410:       PetscInt        numFaces, faceSize, coneSize,f;
5411:       PetscInt       *closure = NULL, closureSize, cl, numCorners = 0;

5413:       if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
5414:       DMPlexGetConeSize(dm, c, &coneSize);
5415:       DMPlexGetCone(dm, c, &cone);
5416:       DMPlexGetConeOrientation(dm, c, &ornt);
5417:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5418:       for (cl = 0; cl < closureSize*2; cl += 2) {
5419:         const PetscInt p = closure[cl];
5420:         if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
5421:       }
5422:       DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
5423:       if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
5424:       for (f = 0; f < numFaces; ++f) {
5425:         PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;

5427:         DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
5428:         for (cl = 0; cl < fclosureSize*2; cl += 2) {
5429:           const PetscInt p = fclosure[cl];
5430:           if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
5431:         }
5432:         if (fnumCorners != faceSize) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d (%d) of cell %d has %d vertices but should have %d", cone[f], f, c, fnumCorners, faceSize);
5433:         for (v = 0; v < fnumCorners; ++v) {
5434:           if (fclosure[v] != faces[f*faceSize+v]) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d (%d) of cell %d vertex %d, %d != %d", cone[f], f, c, v, fclosure[v], faces[f*faceSize+v]);
5435:         }
5436:         DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
5437:       }
5438:       DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
5439:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5440:     }
5441:   }
5442:   return(0);
5443: }

5447: /* Pointwise interpolation
5448:      Just code FEM for now
5449:      u^f = I u^c
5450:      sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
5451:      u^f_i = sum_j psi^f_i I phi^c_j u^c_j
5452:      I_{ij} = psi^f_i phi^c_j
5453: */
5454: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
5455: {
5456:   PetscSection   gsc, gsf;
5457:   PetscInt       m, n;
5458:   void          *ctx;
5459:   DM             cdm;
5460:   PetscBool      regular;

5464:   DMGetDefaultGlobalSection(dmFine, &gsf);
5465:   PetscSectionGetConstrainedStorageSize(gsf, &m);
5466:   DMGetDefaultGlobalSection(dmCoarse, &gsc);
5467:   PetscSectionGetConstrainedStorageSize(gsc, &n);

5469:   MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
5470:   MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
5471:   MatSetType(*interpolation, dmCoarse->mattype);
5472:   DMGetApplicationContext(dmFine, &ctx);

5474:   DMGetCoarseDM(dmFine, &cdm);
5475:   DMPlexGetRegularRefinement(dmFine, &regular);
5476:   if (regular && cdm == dmCoarse) {DMPlexComputeInterpolatorNested(dmCoarse, dmFine, *interpolation, ctx);}
5477:   else                            {DMPlexComputeInterpolatorGeneral(dmCoarse, dmFine, *interpolation, ctx);}
5478:   MatViewFromOptions(*interpolation, NULL, "-interp_mat_view");
5479:   /* Use naive scaling */
5480:   DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
5481:   return(0);
5482: }

5486: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat)
5487: {
5489:   VecScatter     ctx;

5492:   DMPlexComputeInjectorFEM(dmCoarse, dmFine, &ctx, NULL);
5493:   MatCreateScatter(PetscObjectComm((PetscObject)ctx), ctx, mat);
5494:   VecScatterDestroy(&ctx);
5495:   return(0);
5496: }

5500: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
5501: {
5502:   PetscSection   section;
5503:   IS            *bcPoints, *bcComps;
5504:   PetscBool     *isFE;
5505:   PetscInt      *bcFields, *numComp, *numDof;
5506:   PetscInt       depth, dim, numBd, numBC = 0, numFields, bd, bc = 0, f;
5507:   PetscInt       cStart, cEnd, cEndInterior;

5511:   DMGetNumFields(dm, &numFields);
5512:   if (!numFields) return(0);
5513:   /* FE and FV boundary conditions are handled slightly differently */
5514:   PetscMalloc1(numFields, &isFE);
5515:   for (f = 0; f < numFields; ++f) {
5516:     PetscObject  obj;
5517:     PetscClassId id;

5519:     DMGetField(dm, f, &obj);
5520:     PetscObjectGetClassId(obj, &id);
5521:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
5522:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
5523:     else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
5524:   }
5525:   /* Allocate boundary point storage for FEM boundaries */
5526:   DMPlexGetDepth(dm, &depth);
5527:   DMGetDimension(dm, &dim);
5528:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
5529:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
5530:   DMGetNumBoundary(dm, &numBd);
5531:   for (bd = 0; bd < numBd; ++bd) {
5532:     PetscInt  field;
5533:     PetscBool isEssential;

5535:     DMGetBoundary(dm, bd, &isEssential, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL, NULL);
5536:     if (isFE[field] && isEssential) ++numBC;
5537:   }
5538:   /* Add ghost cell boundaries for FVM */
5539:   for (f = 0; f < numFields; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
5540:   PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);
5541:   /* Constrain ghost cells for FV */
5542:   for (f = 0; f < numFields; ++f) {
5543:     PetscInt *newidx, c;

5545:     if (isFE[f] || cEndInterior < 0) continue;
5546:     PetscMalloc1(cEnd-cEndInterior,&newidx);
5547:     for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
5548:     bcFields[bc] = f;
5549:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5550:   }
5551:   /* Handle FEM Dirichlet boundaries */
5552:   for (bd = 0; bd < numBd; ++bd) {
5553:     const char     *bdLabel;
5554:     DMLabel         label;
5555:     const PetscInt *comps;
5556:     const PetscInt *values;
5557:     PetscInt        bd2, field, numComps, numValues;
5558:     PetscBool       isEssential, duplicate = PETSC_FALSE;

5560:     DMGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
5561:     if (!isFE[field]) continue;
5562:     DMGetLabel(dm, bdLabel, &label);
5563:     /* Only want to modify label once */
5564:     for (bd2 = 0; bd2 < bd; ++bd2) {
5565:       const char *bdname;
5566:       DMGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
5567:       PetscStrcmp(bdname, bdLabel, &duplicate);
5568:       if (duplicate) break;
5569:     }
5570:     if (!duplicate && (isFE[field])) {
5571:       /* don't complete cells, which are just present to give orientation to the boundary */
5572:       DMPlexLabelComplete(dm, label);
5573:     }
5574:     /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
5575:     if (isEssential) {
5576:       PetscInt       *newidx;
5577:       PetscInt        n, newn = 0, p, v;

5579:       bcFields[bc] = field;
5580:       if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
5581:       for (v = 0; v < numValues; ++v) {
5582:         IS              tmp;
5583:         const PetscInt *idx;

5585:         DMGetStratumIS(dm, bdLabel, values[v], &tmp);
5586:         if (!tmp) continue;
5587:         ISGetLocalSize(tmp, &n);
5588:         ISGetIndices(tmp, &idx);
5589:         if (isFE[field]) {
5590:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
5591:         } else {
5592:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
5593:         }
5594:         ISRestoreIndices(tmp, &idx);
5595:         ISDestroy(&tmp);
5596:       }
5597:       PetscMalloc1(newn,&newidx);
5598:       newn = 0;
5599:       for (v = 0; v < numValues; ++v) {
5600:         IS              tmp;
5601:         const PetscInt *idx;

5603:         DMGetStratumIS(dm, bdLabel, values[v], &tmp);
5604:         if (!tmp) continue;
5605:         ISGetLocalSize(tmp, &n);
5606:         ISGetIndices(tmp, &idx);
5607:         if (isFE[field]) {
5608:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
5609:         } else {
5610:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
5611:         }
5612:         ISRestoreIndices(tmp, &idx);
5613:         ISDestroy(&tmp);
5614:       }
5615:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5616:     }
5617:   }
5618:   /* Handle discretization */
5619:   PetscCalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
5620:   for (f = 0; f < numFields; ++f) {
5621:     PetscObject obj;

5623:     DMGetField(dm, f, &obj);
5624:     if (isFE[f]) {
5625:       PetscFE         fe = (PetscFE) obj;
5626:       const PetscInt *numFieldDof;
5627:       PetscInt        d;

5629:       PetscFEGetNumComponents(fe, &numComp[f]);
5630:       PetscFEGetNumDof(fe, &numFieldDof);
5631:       for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
5632:     } else {
5633:       PetscFV fv = (PetscFV) obj;

5635:       PetscFVGetNumComponents(fv, &numComp[f]);
5636:       numDof[f*(dim+1)+dim] = numComp[f];
5637:     }
5638:   }
5639:   for (f = 0; f < numFields; ++f) {
5640:     PetscInt d;
5641:     for (d = 1; d < dim; ++d) {
5642:       if ((numDof[f*(dim+1)+d] > 0) && (depth < dim)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
5643:     }
5644:   }
5645:   DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section);
5646:   for (f = 0; f < numFields; ++f) {
5647:     PetscFE     fe;
5648:     const char *name;

5650:     DMGetField(dm, f, (PetscObject *) &fe);
5651:     PetscObjectGetName((PetscObject) fe, &name);
5652:     PetscSectionSetFieldName(section, f, name);
5653:   }
5654:   DMSetDefaultSection(dm, section);
5655:   PetscSectionDestroy(&section);
5656:   for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
5657:   PetscFree3(bcFields,bcPoints,bcComps);
5658:   PetscFree2(numComp,numDof);
5659:   PetscFree(isFE);
5660:   return(0);
5661: }

5665: /*@
5666:   DMPlexGetRegularRefinement - Get the flag indicating that this mesh was obtained by regular refinement from its coarse mesh

5668:   Input Parameter:
5669: . dm - The DMPlex object

5671:   Output Parameter:
5672: . regular - The flag

5674:   Level: intermediate

5676: .seealso: DMPlexSetRegularRefinement()
5677: @*/
5678: PetscErrorCode DMPlexGetRegularRefinement(DM dm, PetscBool *regular)
5679: {
5683:   *regular = ((DM_Plex *) dm->data)->regularRefinement;
5684:   return(0);
5685: }

5689: /*@
5690:   DMPlexSetRegularRefinement - Set the flag indicating that this mesh was obtained by regular refinement from its coarse mesh

5692:   Input Parameters:
5693: + dm - The DMPlex object
5694: - regular - The flag

5696:   Level: intermediate

5698: .seealso: DMPlexGetRegularRefinement()
5699: @*/
5700: PetscErrorCode DMPlexSetRegularRefinement(DM dm, PetscBool regular)
5701: {
5704:   ((DM_Plex *) dm->data)->regularRefinement = regular;
5705:   return(0);
5706: }

5708: /* anchors */
5711: /*@
5712:   DMPlexGetAnchors - Get the layout of the anchor (point-to-point) constraints.  Typically, the user will not have to
5713:   call DMPlexGetAnchors() directly: if there are anchors, then DMPlexGetAnchors() is called during DMGetConstraints().

5715:   not collective

5717:   Input Parameters:
5718: . dm - The DMPlex object

5720:   Output Parameters:
5721: + anchorSection - If not NULL, set to the section describing which points anchor the constrained points.
5722: - anchorIS - If not NULL, set to the list of anchors indexed by anchorSection


5725:   Level: intermediate

5727: .seealso: DMPlexSetAnchors(), DMGetConstraints(), DMSetConstraints()
5728: @*/
5729: PetscErrorCode DMPlexGetAnchors(DM dm, PetscSection *anchorSection, IS *anchorIS)
5730: {
5731:   DM_Plex *plex = (DM_Plex *)dm->data;

5736:   if (!plex->anchorSection && !plex->anchorIS && plex->createanchors) {(*plex->createanchors)(dm);}
5737:   if (anchorSection) *anchorSection = plex->anchorSection;
5738:   if (anchorIS) *anchorIS = plex->anchorIS;
5739:   return(0);
5740: }

5744: /*@
5745:   DMPlexSetAnchors - Set the layout of the local anchor (point-to-point) constraints.  Unlike boundary conditions,
5746:   when a point's degrees of freedom in a section are constrained to an outside value, the anchor constraints set a
5747:   point's degrees of freedom to be a linear combination of other points' degrees of freedom.

5749:   After specifying the layout of constraints with DMPlexSetAnchors(), one specifies the constraints by calling
5750:   DMGetConstraints() and filling in the entries in the constraint matrix.

5752:   collective on dm

5754:   Input Parameters:
5755: + dm - The DMPlex object
5756: . anchorSection - The section that describes the mapping from constrained points to the anchor points listed in anchorIS.  Must have a local communicator (PETSC_COMM_SELF or derivative).
5757: - anchorIS - The list of all anchor points.  Must have a local communicator (PETSC_COMM_SELF or derivative).

5759:   The reference counts of anchorSection and anchorIS are incremented.

5761:   Level: intermediate

5763: .seealso: DMPlexGetAnchors(), DMGetConstraints(), DMSetConstraints()
5764: @*/
5765: PetscErrorCode DMPlexSetAnchors(DM dm, PetscSection anchorSection, IS anchorIS)
5766: {
5767:   DM_Plex        *plex = (DM_Plex *)dm->data;
5768:   PetscMPIInt    result;

5773:   if (anchorSection) {
5775:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorSection),&result);
5776:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor section must have local communicator");
5777:   }
5778:   if (anchorIS) {
5780:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorIS),&result);
5781:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor IS must have local communicator");
5782:   }

5784:   PetscObjectReference((PetscObject)anchorSection);
5785:   PetscSectionDestroy(&plex->anchorSection);
5786:   plex->anchorSection = anchorSection;

5788:   PetscObjectReference((PetscObject)anchorIS);
5789:   ISDestroy(&plex->anchorIS);
5790:   plex->anchorIS = anchorIS;

5792: #if defined(PETSC_USE_DEBUG)
5793:   if (anchorIS && anchorSection) {
5794:     PetscInt size, a, pStart, pEnd;
5795:     const PetscInt *anchors;

5797:     PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5798:     ISGetLocalSize(anchorIS,&size);
5799:     ISGetIndices(anchorIS,&anchors);
5800:     for (a = 0; a < size; a++) {
5801:       PetscInt p;

5803:       p = anchors[a];
5804:       if (p >= pStart && p < pEnd) {
5805:         PetscInt dof;

5807:         PetscSectionGetDof(anchorSection,p,&dof);
5808:         if (dof) {
5809:           PetscErrorCode ierr2;

5811:           ierr2 = ISRestoreIndices(anchorIS,&anchors);CHKERRQ(ierr2);
5812:           SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point %d cannot be constrained and an anchor",p);
5813:         }
5814:       }
5815:     }
5816:     ISRestoreIndices(anchorIS,&anchors);
5817:   }
5818: #endif
5819:   /* reset the generic constraints */
5820:   DMSetDefaultConstraints(dm,NULL,NULL);
5821:   return(0);
5822: }

5826: static PetscErrorCode DMPlexCreateConstraintSection_Anchors(DM dm, PetscSection section, PetscSection *cSec)
5827: {
5828:   PetscSection anchorSection;
5829:   PetscInt pStart, pEnd, sStart, sEnd, p, dof, numFields, f;

5834:   DMPlexGetAnchors(dm,&anchorSection,NULL);
5835:   PetscSectionCreate(PETSC_COMM_SELF,cSec);
5836:   PetscSectionGetNumFields(section,&numFields);
5837:   if (numFields) {
5838:     PetscInt f;
5839:     PetscSectionSetNumFields(*cSec,numFields);

5841:     for (f = 0; f < numFields; f++) {
5842:       PetscInt numComp;

5844:       PetscSectionGetFieldComponents(section,f,&numComp);
5845:       PetscSectionSetFieldComponents(*cSec,f,numComp);
5846:     }
5847:   }
5848:   PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5849:   PetscSectionGetChart(section,&sStart,&sEnd);
5850:   pStart = PetscMax(pStart,sStart);
5851:   pEnd   = PetscMin(pEnd,sEnd);
5852:   pEnd   = PetscMax(pStart,pEnd);
5853:   PetscSectionSetChart(*cSec,pStart,pEnd);
5854:   for (p = pStart; p < pEnd; p++) {
5855:     PetscSectionGetDof(anchorSection,p,&dof);
5856:     if (dof) {
5857:       PetscSectionGetDof(section,p,&dof);
5858:       PetscSectionSetDof(*cSec,p,dof);
5859:       for (f = 0; f < numFields; f++) {
5860:         PetscSectionGetFieldDof(section,p,f,&dof);
5861:         PetscSectionSetFieldDof(*cSec,p,f,dof);
5862:       }
5863:     }
5864:   }
5865:   PetscSectionSetUp(*cSec);
5866:   return(0);
5867: }

5871: static PetscErrorCode DMPlexCreateConstraintMatrix_Anchors(DM dm, PetscSection section, PetscSection cSec, Mat *cMat)
5872: {
5873:   PetscSection aSec;
5874:   PetscInt pStart, pEnd, p, dof, aDof, aOff, off, nnz, annz, m, n, q, a, offset, *i, *j;
5875:   const PetscInt *anchors;
5876:   PetscInt numFields, f;
5877:   IS aIS;

5882:   PetscSectionGetStorageSize(cSec, &m);
5883:   PetscSectionGetStorageSize(section, &n);
5884:   MatCreate(PETSC_COMM_SELF,cMat);
5885:   MatSetSizes(*cMat,m,n,m,n);
5886:   MatSetType(*cMat,MATSEQAIJ);
5887:   DMPlexGetAnchors(dm,&aSec,&aIS);
5888:   ISGetIndices(aIS,&anchors);
5889:   /* cSec will be a subset of aSec and section */
5890:   PetscSectionGetChart(cSec,&pStart,&pEnd);
5891:   PetscMalloc1(m+1,&i);
5892:   i[0] = 0;
5893:   PetscSectionGetNumFields(section,&numFields);
5894:   for (p = pStart; p < pEnd; p++) {
5895:     PetscInt rDof, rOff, r;

5897:     PetscSectionGetDof(aSec,p,&rDof);
5898:     if (!rDof) continue;
5899:     PetscSectionGetOffset(aSec,p,&rOff);
5900:     if (numFields) {
5901:       for (f = 0; f < numFields; f++) {
5902:         annz = 0;
5903:         for (r = 0; r < rDof; r++) {
5904:           a = anchors[rOff + r];
5905:           PetscSectionGetFieldDof(section,a,f,&aDof);
5906:           annz += aDof;
5907:         }
5908:         PetscSectionGetFieldDof(cSec,p,f,&dof);
5909:         PetscSectionGetFieldOffset(cSec,p,f,&off);
5910:         for (q = 0; q < dof; q++) {
5911:           i[off + q + 1] = i[off + q] + annz;
5912:         }
5913:       }
5914:     }
5915:     else {
5916:       annz = 0;
5917:       for (q = 0; q < dof; q++) {
5918:         a = anchors[off + q];
5919:         PetscSectionGetDof(section,a,&aDof);
5920:         annz += aDof;
5921:       }
5922:       PetscSectionGetDof(cSec,p,&dof);
5923:       PetscSectionGetOffset(cSec,p,&off);
5924:       for (q = 0; q < dof; q++) {
5925:         i[off + q + 1] = i[off + q] + annz;
5926:       }
5927:     }
5928:   }
5929:   nnz = i[m];
5930:   PetscMalloc1(nnz,&j);
5931:   offset = 0;
5932:   for (p = pStart; p < pEnd; p++) {
5933:     if (numFields) {
5934:       for (f = 0; f < numFields; f++) {
5935:         PetscSectionGetFieldDof(cSec,p,f,&dof);
5936:         for (q = 0; q < dof; q++) {
5937:           PetscInt rDof, rOff, r;
5938:           PetscSectionGetDof(aSec,p,&rDof);
5939:           PetscSectionGetOffset(aSec,p,&rOff);
5940:           for (r = 0; r < rDof; r++) {
5941:             PetscInt s;

5943:             a = anchors[rOff + r];
5944:             PetscSectionGetFieldDof(section,a,f,&aDof);
5945:             PetscSectionGetFieldOffset(section,a,f,&aOff);
5946:             for (s = 0; s < aDof; s++) {
5947:               j[offset++] = aOff + s;
5948:             }
5949:           }
5950:         }
5951:       }
5952:     }
5953:     else {
5954:       PetscSectionGetDof(cSec,p,&dof);
5955:       for (q = 0; q < dof; q++) {
5956:         PetscInt rDof, rOff, r;
5957:         PetscSectionGetDof(aSec,p,&rDof);
5958:         PetscSectionGetOffset(aSec,p,&rOff);
5959:         for (r = 0; r < rDof; r++) {
5960:           PetscInt s;

5962:           a = anchors[rOff + r];
5963:           PetscSectionGetDof(section,a,&aDof);
5964:           PetscSectionGetOffset(section,a,&aOff);
5965:           for (s = 0; s < aDof; s++) {
5966:             j[offset++] = aOff + s;
5967:           }
5968:         }
5969:       }
5970:     }
5971:   }
5972:   MatSeqAIJSetPreallocationCSR(*cMat,i,j,NULL);
5973:   PetscFree(i);
5974:   PetscFree(j);
5975:   ISRestoreIndices(aIS,&anchors);
5976:   return(0);
5977: }

5981: PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm)
5982: {
5983:   DM_Plex        *plex = (DM_Plex *)dm->data;
5984:   PetscSection   anchorSection, section, cSec;
5985:   Mat            cMat;

5990:   DMPlexGetAnchors(dm,&anchorSection,NULL);
5991:   if (anchorSection) {
5992:     PetscDS  ds;
5993:     PetscInt nf;

5995:     DMGetDefaultSection(dm,&section);
5996:     DMPlexCreateConstraintSection_Anchors(dm,section,&cSec);
5997:     DMPlexCreateConstraintMatrix_Anchors(dm,section,cSec,&cMat);
5998:     DMGetDS(dm,&ds);
5999:     PetscDSGetNumFields(ds,&nf);
6000:     if (nf && plex->computeanchormatrix) {(*plex->computeanchormatrix)(dm,section,cSec,cMat);}
6001:     DMSetDefaultConstraints(dm,cSec,cMat);
6002:     PetscSectionDestroy(&cSec);
6003:     MatDestroy(&cMat);
6004:   }
6005:   return(0);
6006: }