Actual source code: plex.c

petsc-3.7.0 2016-04-25
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) dm); /* viewer drops reference */
 84:     PetscObjectReference((PetscObject) v);  /* viewer drops reference */
 85:     PetscViewerVTKAddField(viewer, (PetscObject) dm, DMPlexVTKWriteAll, ft, (PetscObject) v);
 86:   } else if (ishdf5) {
 87: #if defined(PETSC_HAVE_HDF5)
 88:     VecView_Plex_Local_HDF5(v, viewer);
 89: #else
 90:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
 91: #endif
 92:   } else {
 93:     if (isseq) {VecView_Seq(v, viewer);}
 94:     else       {VecView_MPI(v, viewer);}
 95:   }
 96:   return(0);
 97: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

839:   Not collective

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

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

848:   Level: beginner

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

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

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

868:   Not collective

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

875:   Output Parameters:

877:   Level: beginner

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

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

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

898:   Not collective

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

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

907:   Level: beginner

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

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

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

928:   Not collective

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

935:   Output Parameter:

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

940:   Level: beginner

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

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

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

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

962:   Not collective

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

969:   Output Parameter:

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

974:   Level: beginner

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

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

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

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

998:   Not collective

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

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

1007:   Level: beginner

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

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

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

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

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

1036:   Not collective

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

1043:   Output Parameter:

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

1048:   Level: beginner

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

1061:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1062:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1064:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1065:   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);
1066:   for (c = 0; c < dof; ++c) {
1067:     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);
1068:     mesh->cones[off+c] = cone[c];
1069:   }
1070:   return(0);
1071: }

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

1078:   Not collective

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

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

1090:   Level: beginner

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

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

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

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

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

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

1126:   Not collective

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

1136:   Output Parameter:

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

1141:   Level: beginner

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

1154:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1155:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1157:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1158:   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);
1159:   for (c = 0; c < dof; ++c) {
1160:     PetscInt cdof, o = coneOrientation[c];

1162:     PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1163:     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);
1164:     mesh->coneOrientations[off+c] = o;
1165:   }
1166:   return(0);
1167: }

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

1180:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1181:   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);
1182:   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);
1183:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1184:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1185:   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);
1186:   mesh->cones[off+conePos] = conePoint;
1187:   return(0);
1188: }

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

1201:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1202:   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);
1203:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1204:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1205:   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);
1206:   mesh->coneOrientations[off+conePos] = coneOrientation;
1207:   return(0);
1208: }

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

1215:   Not collective

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

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

1224:   Level: beginner

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

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

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

1245:   Not collective

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

1252:   Output Parameter:

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

1257:   Level: beginner

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

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

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

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

1279:   Not collective

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

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

1288:   Level: beginner

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

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

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

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

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

1317:   Not collective

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

1324:   Output Parameter:

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

1329:   Level: beginner

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

1342:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1343:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1345:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1346:   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);
1347:   for (c = 0; c < dof; ++c) {
1348:     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);
1349:     mesh->supports[off+c] = support[c];
1350:   }
1351:   return(0);
1352: }

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

1365:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1366:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1367:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1368:   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);
1369:   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);
1370:   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);
1371:   mesh->supports[off+supportPos] = supportPoint;
1372:   return(0);
1373: }

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

1380:   Not collective

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

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

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

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

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

1401:   Level: beginner

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

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

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

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

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

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

1522: /*@C
1523:   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

1525:   Not collective

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

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

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

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

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

1547:   Level: beginner

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

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

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

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

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

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

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

1679:   Not collective

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

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

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

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

1697:   Level: beginner

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

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

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

1719:   Not collective

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

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

1728:   Level: beginner

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

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

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

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

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

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

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

1782:   Not collective

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

1787:   Output Parameter:

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

1792:   Level: beginner

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

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

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

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

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

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

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

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

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

1857:   Collective on dm

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

1862:   Output Parameter:

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

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

1872:   Level: beginner

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

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

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

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

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

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

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

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

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

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

1968:   Not Collective

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

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

1979:   Level: intermediate

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

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

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

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

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

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

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

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

2043:   Not Collective

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

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

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

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

2060:   Level: intermediate

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

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

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

2084:   Not Collective

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

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

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

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

2101:   Level: intermediate

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


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

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

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

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

2138:       DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
2139:       for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) {
2140:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2141:           offsets[p*(depth+2)+d+1] = i;
2142:           break;
2143:         }
2144:       }
2145:       if (i == closureSize) offsets[p*(depth+2)+d+1] = i;
2146:     }
2147:     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);
2148:   }
2149:   for (d = 0; d < depth+1; ++d) {
2150:     PetscInt dof;

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

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

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

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

2193:   Not Collective

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

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

2204:   Level: intermediate

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

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

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

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

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

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

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

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

2268:   Not Collective

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

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

2279:   Level: intermediate

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

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

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

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

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

2309:   Not Collective

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

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

2320:   Level: intermediate

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

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

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


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

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

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

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

2363:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
2364:       for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) {
2365:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2366:           offsets[p*(height+2)+h+1] = i;
2367:           break;
2368:         }
2369:       }
2370:       if (i == closureSize) offsets[p*(height+2)+h+1] = i;
2371:     }
2372:     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);
2373:   }
2374:   for (h = 0; h < height+1; ++h) {
2375:     PetscInt dof;

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

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

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

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

2418:   Not Collective

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

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

2427:   Level: intermediate

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

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


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

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

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

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

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

2556:   Not Collective

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

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

2564:   Level: developer

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

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

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

2586:   Not Collective

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

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

2594:   Level: developer

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

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

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

2619:   Not Collective

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

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

2629:   Level: developer

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

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

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

2662:   Not Collective

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

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

2672:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2951:   Not Collective

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

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

2968:   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
2969:   number of dof for field 0 on each edge.

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

2973:   Level: developer

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

2978: .keywords: mesh, elements
2979: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
2980: @*/
2981: 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)
2982: {
2983:   PetscSection   aSec;

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

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

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

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

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

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

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

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

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

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

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

3069: /******************************** FEM Support **********************************/

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3224:   Not collective

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

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

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

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

3242:   Level: intermediate

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

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

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

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

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

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

3334:   Not collective

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

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

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

3350:   Level: intermediate

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

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

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

3370: 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[])
3371: {
3372:   PetscInt        cdof;   /* The number of constraints on this point */
3373:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3374:   PetscScalar    *a;
3375:   PetscInt        off, cind = 0, k;
3376:   PetscErrorCode  ierr;

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

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

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

3446: 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[])
3447: {
3448:   PetscScalar    *a;
3449:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3450:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3451:   PetscInt        cind = 0, k, c;
3452:   PetscErrorCode  ierr;

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

3491: 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[])
3492: {
3493:   PetscScalar    *a;
3494:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3495:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3496:   PetscInt        cind = 0, k, c;
3497:   PetscErrorCode  ierr;

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

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

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

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

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

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

3593:   Not collective

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

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

3606:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3980: 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)
3981: {
3982:   Mat             cMat;
3983:   PetscSection    aSec, cSec;
3984:   IS              aIS;
3985:   PetscInt        aStart = -1, aEnd = -1;
3986:   const PetscInt  *anchors;
3987:   PetscInt        numFields, f, p, q, newP = 0;
3988:   PetscInt        newNumPoints = 0, newNumIndices = 0;
3989:   PetscInt        *newPoints, *indices, *newIndices;
3990:   PetscInt        maxAnchor, maxDof;
3991:   PetscInt        newOffsets[32];
3992:   PetscInt        *pointMatOffsets[32];
3993:   PetscInt        *newPointOffsets[32];
3994:   PetscScalar     *pointMat[32];
3995:   PetscScalar     *newValues=NULL,*tmpValues;
3996:   PetscBool       anyConstrained = PETSC_FALSE;
3997:   PetscErrorCode  ierr;

4002:   PetscSectionGetNumFields(section, &numFields);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4629:   Not collective

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

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

4643:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5034:   Level: developer

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

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

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

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

5066:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5244:   Level: developer

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

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

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

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

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

5332:   Level: developer

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

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

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

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

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

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

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

5392:   Level: developer

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

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

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

5428:         DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
5429:         for (cl = 0; cl < fclosureSize*2; cl += 2) {
5430:           const PetscInt p = fclosure[cl];
5431:           if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
5432:         }
5433:         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);
5434:         for (v = 0; v < fnumCorners; ++v) {
5435:           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]);
5436:         }
5437:         DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
5438:       }
5439:       DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
5440:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5441:     }
5442:   }
5443:   return(0);
5444: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5636:       PetscFVGetNumComponents(fv, &numComp[f]);
5637:       numDof[f*(dim+1)+dim] = numComp[f];
5638:     }
5639:   }
5640:   for (f = 0; f < numFields; ++f) {
5641:     PetscInt d;
5642:     for (d = 1; d < dim; ++d) {
5643:       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.");
5644:     }
5645:   }
5646:   DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section);
5647:   for (f = 0; f < numFields; ++f) {
5648:     PetscFE     fe;
5649:     const char *name;

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

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

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

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

5675:   Level: intermediate

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

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

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

5697:   Level: intermediate

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

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

5716:   not collective

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

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


5726:   Level: intermediate

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

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

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

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

5753:   collective on dm

5755:   Input Parameters:
5756: + dm - The DMPlex object
5757: . 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).
5758: - anchorIS - The list of all anchor points.  Must have a local communicator (PETSC_COMM_SELF or derivative).

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

5762:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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