Actual source code: plex.c

petsc-master 2016-06-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, mbs;
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:   MatGetBlockSize(*J, &mbs);
778:   if (mbs > 1) bs = mbs;
779:   PetscStrcmp(mtype, MATSHELL, &isShell);
780:   PetscStrcmp(mtype, MATBAIJ, &isBlock);
781:   PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
782:   PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
783:   PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
784:   PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
785:   PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
786:   if (!isShell) {
787:     PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
788:     PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin;
789:     PetscInt  pStart, pEnd, p, dof, cdof;

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

795:       PetscSectionGetDof(sectionGlobal, p, &dof);
796:       PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
797:       dof  = dof < 0 ? -(dof+1) : dof;
798:       bdof = cdof && (dof-cdof) ? 1 : dof;
799:       if (dof) {
800:         if (bs < 0)          {bs = bdof;}
801:         else if (bs != bdof) {bs = 1; break;}
802:       }
803:     }
804:     /* Must have same blocksize on all procs (some might have no points) */
805:     bsLocal = bs;
806:     MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
807:     bsLocal = bs < 0 ? bsMax : bs;
808:     MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
809:     if (bsMin != bsMax) {bs = 1;}
810:     else                {bs = bsMax;}
811:     PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
812:     DMPlexPreallocateOperator(dm, bs, dnz, onz, dnzu, onzu, *J, fillMatrix);
813:     PetscFree4(dnz, onz, dnzu, onzu);

815:     /* Set localtoglobalmapping on the matrix for MatSetValuesLocal() to work */
816:     DMGetLocalToGlobalMapping(dm,&ltog);
817:     MatSetLocalToGlobalMapping(*J,ltog,ltog);
818:   }
819:   return(0);
820: }

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

827:   Not collective

829:   Input Parameter:
830: . mesh - The DMPlex

832:   Output Parameters:
833: + pStart - The first mesh point
834: - pEnd   - The upper bound for mesh points

836:   Level: beginner

838: .seealso: DMPlexCreate(), DMPlexSetChart()
839: @*/
840: PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
841: {
842:   DM_Plex       *mesh = (DM_Plex*) dm->data;

847:   PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
848:   return(0);
849: }

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

856:   Not collective

858:   Input Parameters:
859: + mesh - The DMPlex
860: . pStart - The first mesh point
861: - pEnd   - The upper bound for mesh points

863:   Output Parameters:

865:   Level: beginner

867: .seealso: DMPlexCreate(), DMPlexGetChart()
868: @*/
869: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
870: {
871:   DM_Plex       *mesh = (DM_Plex*) dm->data;

876:   PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
877:   PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
878:   return(0);
879: }

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

886:   Not collective

888:   Input Parameters:
889: + mesh - The DMPlex
890: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

892:   Output Parameter:
893: . size - The cone size for point p

895:   Level: beginner

897: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
898: @*/
899: PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size)
900: {
901:   DM_Plex       *mesh = (DM_Plex*) dm->data;

907:   PetscSectionGetDof(mesh->coneSection, p, size);
908:   return(0);
909: }

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

916:   Not collective

918:   Input Parameters:
919: + mesh - The DMPlex
920: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
921: - size - The cone size for point p

923:   Output Parameter:

925:   Note:
926:   This should be called after DMPlexSetChart().

928:   Level: beginner

930: .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart()
931: @*/
932: PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size)
933: {
934:   DM_Plex       *mesh = (DM_Plex*) dm->data;

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

941:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
942:   return(0);
943: }

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

950:   Not collective

952:   Input Parameters:
953: + mesh - The DMPlex
954: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
955: - size - The additional cone size for point p

957:   Output Parameter:

959:   Note:
960:   This should be called after DMPlexSetChart().

962:   Level: beginner

964: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart()
965: @*/
966: PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size)
967: {
968:   DM_Plex       *mesh = (DM_Plex*) dm->data;
969:   PetscInt       csize;

974:   PetscSectionAddDof(mesh->coneSection, p, size);
975:   PetscSectionGetDof(mesh->coneSection, p, &csize);

977:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
978:   return(0);
979: }

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

986:   Not collective

988:   Input Parameters:
989: + mesh - The DMPlex
990: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

995:   Level: beginner

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

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

1003: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart()
1004: @*/
1005: PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
1006: {
1007:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1008:   PetscInt       off;

1014:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1015:   *cone = &mesh->cones[off];
1016:   return(0);
1017: }

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

1024:   Not collective

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

1031:   Output Parameter:

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

1036:   Level: beginner

1038: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1039: @*/
1040: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
1041: {
1042:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1043:   PetscInt       pStart, pEnd;
1044:   PetscInt       dof, off, c;

1049:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1050:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1052:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1053:   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);
1054:   for (c = 0; c < dof; ++c) {
1055:     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);
1056:     mesh->cones[off+c] = cone[c];
1057:   }
1058:   return(0);
1059: }

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

1066:   Not collective

1068:   Input Parameters:
1069: + mesh - The DMPlex
1070: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1078:   Level: beginner

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

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

1086: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
1087: @*/
1088: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
1089: {
1090:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1091:   PetscInt       off;

1096: #if defined(PETSC_USE_DEBUG)
1097:   {
1098:     PetscInt dof;
1099:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1101:   }
1102: #endif
1103:   PetscSectionGetOffset(mesh->coneSection, p, &off);

1105:   *coneOrientation = &mesh->coneOrientations[off];
1106:   return(0);
1107: }

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

1114:   Not collective

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

1124:   Output Parameter:

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

1129:   Level: beginner

1131: .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1132: @*/
1133: PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
1134: {
1135:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1136:   PetscInt       pStart, pEnd;
1137:   PetscInt       dof, off, c;

1142:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1143:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1145:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1146:   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);
1147:   for (c = 0; c < dof; ++c) {
1148:     PetscInt cdof, o = coneOrientation[c];

1150:     PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1151:     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);
1152:     mesh->coneOrientations[off+c] = o;
1153:   }
1154:   return(0);
1155: }

1159: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
1160: {
1161:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1162:   PetscInt       pStart, pEnd;
1163:   PetscInt       dof, off;

1168:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1169:   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);
1170:   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);
1171:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1172:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1173:   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);
1174:   mesh->cones[off+conePos] = conePoint;
1175:   return(0);
1176: }

1180: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1181: {
1182:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1183:   PetscInt       pStart, pEnd;
1184:   PetscInt       dof, off;

1189:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1190:   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);
1191:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1192:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1193:   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);
1194:   mesh->coneOrientations[off+conePos] = coneOrientation;
1195:   return(0);
1196: }

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

1203:   Not collective

1205:   Input Parameters:
1206: + mesh - The DMPlex
1207: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1209:   Output Parameter:
1210: . size - The support size for point p

1212:   Level: beginner

1214: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1215: @*/
1216: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1217: {
1218:   DM_Plex       *mesh = (DM_Plex*) dm->data;

1224:   PetscSectionGetDof(mesh->supportSection, p, size);
1225:   return(0);
1226: }

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

1233:   Not collective

1235:   Input Parameters:
1236: + mesh - The DMPlex
1237: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1238: - size - The support size for point p

1240:   Output Parameter:

1242:   Note:
1243:   This should be called after DMPlexSetChart().

1245:   Level: beginner

1247: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1248: @*/
1249: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1250: {
1251:   DM_Plex       *mesh = (DM_Plex*) dm->data;

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

1258:   mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1259:   return(0);
1260: }

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

1267:   Not collective

1269:   Input Parameters:
1270: + mesh - The DMPlex
1271: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1276:   Level: beginner

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

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

1284: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1285: @*/
1286: PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1287: {
1288:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1289:   PetscInt       off;

1295:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1296:   *support = &mesh->supports[off];
1297:   return(0);
1298: }

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

1305:   Not collective

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

1312:   Output Parameter:

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

1317:   Level: beginner

1319: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1320: @*/
1321: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1322: {
1323:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1324:   PetscInt       pStart, pEnd;
1325:   PetscInt       dof, off, c;

1330:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1331:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1333:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1334:   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);
1335:   for (c = 0; c < dof; ++c) {
1336:     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);
1337:     mesh->supports[off+c] = support[c];
1338:   }
1339:   return(0);
1340: }

1344: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1345: {
1346:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1347:   PetscInt       pStart, pEnd;
1348:   PetscInt       dof, off;

1353:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1354:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1355:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1356:   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);
1357:   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);
1358:   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);
1359:   mesh->supports[off+supportPos] = supportPoint;
1360:   return(0);
1361: }

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

1368:   Not collective

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

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

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

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

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

1389:   Level: beginner

1391: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1392: @*/
1393: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1394: {
1395:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1396:   PetscInt       *closure, *fifo;
1397:   const PetscInt *tmp = NULL, *tmpO = NULL;
1398:   PetscInt        tmpSize, t;
1399:   PetscInt        depth       = 0, maxSize;
1400:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1401:   PetscErrorCode  ierr;

1405:   DMPlexGetDepth(dm, &depth);
1406:   /* This is only 1-level */
1407:   if (useCone) {
1408:     DMPlexGetConeSize(dm, p, &tmpSize);
1409:     DMPlexGetCone(dm, p, &tmp);
1410:     DMPlexGetConeOrientation(dm, p, &tmpO);
1411:   } else {
1412:     DMPlexGetSupportSize(dm, p, &tmpSize);
1413:     DMPlexGetSupport(dm, p, &tmp);
1414:   }
1415:   if (depth == 1) {
1416:     if (*points) {
1417:       closure = *points;
1418:     } else {
1419:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1420:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1421:     }
1422:     closure[0] = p; closure[1] = 0;
1423:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1424:       closure[closureSize]   = tmp[t];
1425:       closure[closureSize+1] = tmpO ? tmpO[t] : 0;
1426:     }
1427:     if (numPoints) *numPoints = closureSize/2;
1428:     if (points)    *points    = closure;
1429:     return(0);
1430:   }
1431:   {
1432:     PetscInt c, coneSeries, s,supportSeries;

1434:     c = mesh->maxConeSize;
1435:     coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1436:     s = mesh->maxSupportSize;
1437:     supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1438:     maxSize = 2*PetscMax(coneSeries,supportSeries);
1439:   }
1440:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1441:   if (*points) {
1442:     closure = *points;
1443:   } else {
1444:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1445:   }
1446:   closure[0] = p; closure[1] = 0;
1447:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1448:     const PetscInt cp = tmp[t];
1449:     const PetscInt co = tmpO ? tmpO[t] : 0;

1451:     closure[closureSize]   = cp;
1452:     closure[closureSize+1] = co;
1453:     fifo[fifoSize]         = cp;
1454:     fifo[fifoSize+1]       = co;
1455:   }
1456:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1457:   while (fifoSize - fifoStart) {
1458:     const PetscInt q   = fifo[fifoStart];
1459:     const PetscInt o   = fifo[fifoStart+1];
1460:     const PetscInt rev = o >= 0 ? 0 : 1;
1461:     const PetscInt off = rev ? -(o+1) : o;

1463:     if (useCone) {
1464:       DMPlexGetConeSize(dm, q, &tmpSize);
1465:       DMPlexGetCone(dm, q, &tmp);
1466:       DMPlexGetConeOrientation(dm, q, &tmpO);
1467:     } else {
1468:       DMPlexGetSupportSize(dm, q, &tmpSize);
1469:       DMPlexGetSupport(dm, q, &tmp);
1470:       tmpO = NULL;
1471:     }
1472:     for (t = 0; t < tmpSize; ++t) {
1473:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1474:       const PetscInt cp = tmp[i];
1475:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1476:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1477:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1478:       PetscInt       co = tmpO ? tmpO[i] : 0;
1479:       PetscInt       c;

1481:       if (rev) {
1482:         PetscInt childSize, coff;
1483:         DMPlexGetConeSize(dm, cp, &childSize);
1484:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1485:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1486:       }
1487:       /* Check for duplicate */
1488:       for (c = 0; c < closureSize; c += 2) {
1489:         if (closure[c] == cp) break;
1490:       }
1491:       if (c == closureSize) {
1492:         closure[closureSize]   = cp;
1493:         closure[closureSize+1] = co;
1494:         fifo[fifoSize]         = cp;
1495:         fifo[fifoSize+1]       = co;
1496:         closureSize           += 2;
1497:         fifoSize              += 2;
1498:       }
1499:     }
1500:     fifoStart += 2;
1501:   }
1502:   if (numPoints) *numPoints = closureSize/2;
1503:   if (points)    *points    = closure;
1504:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1505:   return(0);
1506: }

1510: /*@C
1511:   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

1513:   Not collective

1515:   Input Parameters:
1516: + mesh - The DMPlex
1517: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1518: . orientation - The orientation of the point
1519: . useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
1520: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used

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

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

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

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

1535:   Level: beginner

1537: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1538: @*/
1539: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1540: {
1541:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1542:   PetscInt       *closure, *fifo;
1543:   const PetscInt *tmp = NULL, *tmpO = NULL;
1544:   PetscInt        tmpSize, t;
1545:   PetscInt        depth       = 0, maxSize;
1546:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1547:   PetscErrorCode  ierr;

1551:   DMPlexGetDepth(dm, &depth);
1552:   /* This is only 1-level */
1553:   if (useCone) {
1554:     DMPlexGetConeSize(dm, p, &tmpSize);
1555:     DMPlexGetCone(dm, p, &tmp);
1556:     DMPlexGetConeOrientation(dm, p, &tmpO);
1557:   } else {
1558:     DMPlexGetSupportSize(dm, p, &tmpSize);
1559:     DMPlexGetSupport(dm, p, &tmp);
1560:   }
1561:   if (depth == 1) {
1562:     if (*points) {
1563:       closure = *points;
1564:     } else {
1565:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1566:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1567:     }
1568:     closure[0] = p; closure[1] = ornt;
1569:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1570:       const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1571:       closure[closureSize]   = tmp[i];
1572:       closure[closureSize+1] = tmpO ? tmpO[i] : 0;
1573:     }
1574:     if (numPoints) *numPoints = closureSize/2;
1575:     if (points)    *points    = closure;
1576:     return(0);
1577:   }
1578:   {
1579:     PetscInt c, coneSeries, s,supportSeries;

1581:     c = mesh->maxConeSize;
1582:     coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1583:     s = mesh->maxSupportSize;
1584:     supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1585:     maxSize = 2*PetscMax(coneSeries,supportSeries);
1586:   }
1587:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1588:   if (*points) {
1589:     closure = *points;
1590:   } else {
1591:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1592:   }
1593:   closure[0] = p; closure[1] = ornt;
1594:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1595:     const PetscInt i  = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1596:     const PetscInt cp = tmp[i];
1597:     PetscInt       co = tmpO ? tmpO[i] : 0;

1599:     if (ornt < 0) {
1600:       PetscInt childSize, coff;
1601:       DMPlexGetConeSize(dm, cp, &childSize);
1602:       coff = co < 0 ? -(tmpO[i]+1) : tmpO[i];
1603:       co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1604:     }
1605:     closure[closureSize]   = cp;
1606:     closure[closureSize+1] = co;
1607:     fifo[fifoSize]         = cp;
1608:     fifo[fifoSize+1]       = co;
1609:   }
1610:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1611:   while (fifoSize - fifoStart) {
1612:     const PetscInt q   = fifo[fifoStart];
1613:     const PetscInt o   = fifo[fifoStart+1];
1614:     const PetscInt rev = o >= 0 ? 0 : 1;
1615:     const PetscInt off = rev ? -(o+1) : o;

1617:     if (useCone) {
1618:       DMPlexGetConeSize(dm, q, &tmpSize);
1619:       DMPlexGetCone(dm, q, &tmp);
1620:       DMPlexGetConeOrientation(dm, q, &tmpO);
1621:     } else {
1622:       DMPlexGetSupportSize(dm, q, &tmpSize);
1623:       DMPlexGetSupport(dm, q, &tmp);
1624:       tmpO = NULL;
1625:     }
1626:     for (t = 0; t < tmpSize; ++t) {
1627:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1628:       const PetscInt cp = tmp[i];
1629:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1630:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1631:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1632:       PetscInt       co = tmpO ? tmpO[i] : 0;
1633:       PetscInt       c;

1635:       if (rev) {
1636:         PetscInt childSize, coff;
1637:         DMPlexGetConeSize(dm, cp, &childSize);
1638:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1639:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1640:       }
1641:       /* Check for duplicate */
1642:       for (c = 0; c < closureSize; c += 2) {
1643:         if (closure[c] == cp) break;
1644:       }
1645:       if (c == closureSize) {
1646:         closure[closureSize]   = cp;
1647:         closure[closureSize+1] = co;
1648:         fifo[fifoSize]         = cp;
1649:         fifo[fifoSize+1]       = co;
1650:         closureSize           += 2;
1651:         fifoSize              += 2;
1652:       }
1653:     }
1654:     fifoStart += 2;
1655:   }
1656:   if (numPoints) *numPoints = closureSize/2;
1657:   if (points)    *points    = closure;
1658:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1659:   return(0);
1660: }

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

1667:   Not collective

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

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

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

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

1685:   Level: beginner

1687: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1688: @*/
1689: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1690: {

1697:   DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1698:   if (numPoints) *numPoints = 0;
1699:   return(0);
1700: }

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

1707:   Not collective

1709:   Input Parameter:
1710: . mesh - The DMPlex

1712:   Output Parameters:
1713: + maxConeSize - The maximum number of in-edges
1714: - maxSupportSize - The maximum number of out-edges

1716:   Level: beginner

1718: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1719: @*/
1720: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1721: {
1722:   DM_Plex *mesh = (DM_Plex*) dm->data;

1726:   if (maxConeSize)    *maxConeSize    = mesh->maxConeSize;
1727:   if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1728:   return(0);
1729: }

1733: PetscErrorCode DMSetUp_Plex(DM dm)
1734: {
1735:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1736:   PetscInt       size;

1741:   PetscSectionSetUp(mesh->coneSection);
1742:   PetscSectionGetStorageSize(mesh->coneSection, &size);
1743:   PetscMalloc1(size, &mesh->cones);
1744:   PetscCalloc1(size, &mesh->coneOrientations);
1745:   if (mesh->maxSupportSize) {
1746:     PetscSectionSetUp(mesh->supportSection);
1747:     PetscSectionGetStorageSize(mesh->supportSection, &size);
1748:     PetscMalloc1(size, &mesh->supports);
1749:   }
1750:   return(0);
1751: }

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

1760:   if (subdm) {DMClone(dm, subdm);}
1761:   DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1762:   return(0);
1763: }

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

1770:   Not collective

1772:   Input Parameter:
1773: . mesh - The DMPlex

1775:   Output Parameter:

1777:   Note:
1778:   This should be called after all calls to DMPlexSetCone()

1780:   Level: beginner

1782: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1783: @*/
1784: PetscErrorCode DMPlexSymmetrize(DM dm)
1785: {
1786:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1787:   PetscInt      *offsets;
1788:   PetscInt       supportSize;
1789:   PetscInt       pStart, pEnd, p;

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

1800:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1801:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1802:     for (c = off; c < off+dof; ++c) {
1803:       PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1804:     }
1805:   }
1806:   for (p = pStart; p < pEnd; ++p) {
1807:     PetscInt dof;

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

1811:     mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1812:   }
1813:   PetscSectionSetUp(mesh->supportSection);
1814:   /* Calculate supports */
1815:   PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1816:   PetscMalloc1(supportSize, &mesh->supports);
1817:   PetscCalloc1(pEnd - pStart, &offsets);
1818:   for (p = pStart; p < pEnd; ++p) {
1819:     PetscInt dof, off, c;

1821:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1822:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1823:     for (c = off; c < off+dof; ++c) {
1824:       const PetscInt q = mesh->cones[c];
1825:       PetscInt       offS;

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

1829:       mesh->supports[offS+offsets[q]] = p;
1830:       ++offsets[q];
1831:     }
1832:   }
1833:   PetscFree(offsets);
1834:   return(0);
1835: }

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

1845:   Collective on dm

1847:   Input Parameter:
1848: . mesh - The DMPlex

1850:   Output Parameter:

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

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

1860:   Level: beginner

1862: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1863: @*/
1864: PetscErrorCode DMPlexStratify(DM dm)
1865: {
1866:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1867:   DMLabel        label;
1868:   PetscInt       pStart, pEnd, p;
1869:   PetscInt       numRoots = 0, numLeaves = 0;

1874:   PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1875:   /* Calculate depth */
1876:   DMPlexGetChart(dm, &pStart, &pEnd);
1877:   DMCreateLabel(dm, "depth");
1878:   DMPlexGetDepthLabel(dm, &label);
1879:   /* Initialize roots and count leaves */
1880:   for (p = pStart; p < pEnd; ++p) {
1881:     PetscInt coneSize, supportSize;

1883:     DMPlexGetConeSize(dm, p, &coneSize);
1884:     DMPlexGetSupportSize(dm, p, &supportSize);
1885:     if (!coneSize && supportSize) {
1886:       ++numRoots;
1887:       DMLabelSetValue(label, p, 0);
1888:     } else if (!supportSize && coneSize) {
1889:       ++numLeaves;
1890:     } else if (!supportSize && !coneSize) {
1891:       /* Isolated points */
1892:       DMLabelSetValue(label, p, 0);
1893:     }
1894:   }
1895:   if (numRoots + numLeaves == (pEnd - pStart)) {
1896:     for (p = pStart; p < pEnd; ++p) {
1897:       PetscInt coneSize, supportSize;

1899:       DMPlexGetConeSize(dm, p, &coneSize);
1900:       DMPlexGetSupportSize(dm, p, &supportSize);
1901:       if (!supportSize && coneSize) {
1902:         DMLabelSetValue(label, p, 1);
1903:       }
1904:     }
1905:   } else {
1906:     IS       pointIS;
1907:     PetscInt numPoints = 0, level = 0;

1909:     DMLabelGetStratumIS(label, level, &pointIS);
1910:     if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1911:     while (numPoints) {
1912:       const PetscInt *points;
1913:       const PetscInt  newLevel = level+1;

1915:       ISGetIndices(pointIS, &points);
1916:       for (p = 0; p < numPoints; ++p) {
1917:         const PetscInt  point = points[p];
1918:         const PetscInt *support;
1919:         PetscInt        supportSize, s;

1921:         DMPlexGetSupportSize(dm, point, &supportSize);
1922:         DMPlexGetSupport(dm, point, &support);
1923:         for (s = 0; s < supportSize; ++s) {
1924:           DMLabelSetValue(label, support[s], newLevel);
1925:         }
1926:       }
1927:       ISRestoreIndices(pointIS, &points);
1928:       ++level;
1929:       ISDestroy(&pointIS);
1930:       DMLabelGetStratumIS(label, level, &pointIS);
1931:       if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1932:       else         {numPoints = 0;}
1933:     }
1934:     ISDestroy(&pointIS);
1935:   }
1936:   { /* just in case there is an empty process */
1937:     PetscInt numValues, maxValues = 0, v;

1939:     DMLabelGetNumValues(label,&numValues);
1940:     for (v = 0; v < numValues; v++) {
1941:       IS pointIS;

1943:       DMLabelGetStratumIS(label, v, &pointIS);
1944:       if (pointIS) {
1945:         PetscInt  min, max, numPoints;
1946:         PetscInt  start;
1947:         PetscBool contig;

1949:         ISGetLocalSize(pointIS, &numPoints);
1950:         ISGetMinMax(pointIS, &min, &max);
1951:         ISContiguousLocal(pointIS,min,max+1,&start,&contig);
1952:         if (start == 0 && contig) {
1953:           ISDestroy(&pointIS);
1954:           ISCreateStride(PETSC_COMM_SELF,numPoints,min,1,&pointIS);
1955:           DMLabelSetStratumIS(label, v, pointIS);
1956:         }
1957:       }
1958:       ISDestroy(&pointIS);
1959:     }
1960:     MPI_Allreduce(&numValues,&maxValues,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
1961:     for (v = numValues; v < maxValues; v++) {
1962:       DMLabelAddStratum(label,v);
1963:     }
1964:   }

1966:   DMLabelGetState(label, &mesh->depthState);
1967:   PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1968:   return(0);
1969: }

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

1976:   Not Collective

1978:   Input Parameters:
1979: + dm - The DMPlex object
1980: . numPoints - The number of input points for the join
1981: - points - The input points

1983:   Output Parameters:
1984: + numCoveredPoints - The number of points in the join
1985: - coveredPoints - The points in the join

1987:   Level: intermediate

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

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

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

1997: .keywords: mesh
1998: .seealso: DMPlexRestoreJoin(), DMPlexGetMeet()
1999: @*/
2000: PetscErrorCode DMPlexGetJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2001: {
2002:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2003:   PetscInt      *join[2];
2004:   PetscInt       joinSize, i = 0;
2005:   PetscInt       dof, off, p, c, m;

2013:   DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[0]);
2014:   DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1]);
2015:   /* Copy in support of first point */
2016:   PetscSectionGetDof(mesh->supportSection, points[0], &dof);
2017:   PetscSectionGetOffset(mesh->supportSection, points[0], &off);
2018:   for (joinSize = 0; joinSize < dof; ++joinSize) {
2019:     join[i][joinSize] = mesh->supports[off+joinSize];
2020:   }
2021:   /* Check each successive support */
2022:   for (p = 1; p < numPoints; ++p) {
2023:     PetscInt newJoinSize = 0;

2025:     PetscSectionGetDof(mesh->supportSection, points[p], &dof);
2026:     PetscSectionGetOffset(mesh->supportSection, points[p], &off);
2027:     for (c = 0; c < dof; ++c) {
2028:       const PetscInt point = mesh->supports[off+c];

2030:       for (m = 0; m < joinSize; ++m) {
2031:         if (point == join[i][m]) {
2032:           join[1-i][newJoinSize++] = point;
2033:           break;
2034:         }
2035:       }
2036:     }
2037:     joinSize = newJoinSize;
2038:     i        = 1-i;
2039:   }
2040:   *numCoveredPoints = joinSize;
2041:   *coveredPoints    = join[i];
2042:   DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2043:   return(0);
2044: }

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

2051:   Not Collective

2053:   Input Parameters:
2054: + dm - The DMPlex object
2055: . numPoints - The number of input points for the join
2056: - points - The input points

2058:   Output Parameters:
2059: + numCoveredPoints - The number of points in the join
2060: - coveredPoints - The points in the join

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

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

2068:   Level: intermediate

2070: .keywords: mesh
2071: .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet()
2072: @*/
2073: PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2074: {

2082:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2083:   if (numCoveredPoints) *numCoveredPoints = 0;
2084:   return(0);
2085: }

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

2092:   Not Collective

2094:   Input Parameters:
2095: + dm - The DMPlex object
2096: . numPoints - The number of input points for the join
2097: - points - The input points

2099:   Output Parameters:
2100: + numCoveredPoints - The number of points in the join
2101: - coveredPoints - The points in the join

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

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

2109:   Level: intermediate

2111: .keywords: mesh
2112: .seealso: DMPlexGetJoin(), DMPlexRestoreJoin(), DMPlexGetMeet()
2113: @*/
2114: PetscErrorCode DMPlexGetFullJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2115: {
2116:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2117:   PetscInt      *offsets, **closures;
2118:   PetscInt      *join[2];
2119:   PetscInt       depth = 0, maxSize, joinSize = 0, i = 0;
2120:   PetscInt       p, d, c, m, ms;


2129:   DMPlexGetDepth(dm, &depth);
2130:   PetscCalloc1(numPoints, &closures);
2131:   DMGetWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2132:   ms      = mesh->maxSupportSize;
2133:   maxSize = (ms > 1) ? ((PetscPowInt(ms,depth+1)-1)/(ms-1)) : depth + 1;
2134:   DMGetWorkArray(dm, maxSize, PETSC_INT, &join[0]);
2135:   DMGetWorkArray(dm, maxSize, PETSC_INT, &join[1]);

2137:   for (p = 0; p < numPoints; ++p) {
2138:     PetscInt closureSize;

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

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

2146:       DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
2147:       for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) {
2148:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2149:           offsets[p*(depth+2)+d+1] = i;
2150:           break;
2151:         }
2152:       }
2153:       if (i == closureSize) offsets[p*(depth+2)+d+1] = i;
2154:     }
2155:     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);
2156:   }
2157:   for (d = 0; d < depth+1; ++d) {
2158:     PetscInt dof;

2160:     /* Copy in support of first point */
2161:     dof = offsets[d+1] - offsets[d];
2162:     for (joinSize = 0; joinSize < dof; ++joinSize) {
2163:       join[i][joinSize] = closures[0][(offsets[d]+joinSize)*2];
2164:     }
2165:     /* Check each successive cone */
2166:     for (p = 1; p < numPoints && joinSize; ++p) {
2167:       PetscInt newJoinSize = 0;

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

2173:         for (m = 0; m < joinSize; ++m) {
2174:           if (point == join[i][m]) {
2175:             join[1-i][newJoinSize++] = point;
2176:             break;
2177:           }
2178:         }
2179:       }
2180:       joinSize = newJoinSize;
2181:       i        = 1-i;
2182:     }
2183:     if (joinSize) break;
2184:   }
2185:   *numCoveredPoints = joinSize;
2186:   *coveredPoints    = join[i];
2187:   for (p = 0; p < numPoints; ++p) {
2188:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, NULL, &closures[p]);
2189:   }
2190:   PetscFree(closures);
2191:   DMRestoreWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2192:   DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2193:   return(0);
2194: }

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

2201:   Not Collective

2203:   Input Parameters:
2204: + dm - The DMPlex object
2205: . numPoints - The number of input points for the meet
2206: - points - The input points

2208:   Output Parameters:
2209: + numCoveredPoints - The number of points in the meet
2210: - coveredPoints - The points in the meet

2212:   Level: intermediate

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

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

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

2222: .keywords: mesh
2223: .seealso: DMPlexRestoreMeet(), DMPlexGetJoin()
2224: @*/
2225: PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2226: {
2227:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2228:   PetscInt      *meet[2];
2229:   PetscInt       meetSize, i = 0;
2230:   PetscInt       dof, off, p, c, m;

2238:   DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[0]);
2239:   DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1]);
2240:   /* Copy in cone of first point */
2241:   PetscSectionGetDof(mesh->coneSection, points[0], &dof);
2242:   PetscSectionGetOffset(mesh->coneSection, points[0], &off);
2243:   for (meetSize = 0; meetSize < dof; ++meetSize) {
2244:     meet[i][meetSize] = mesh->cones[off+meetSize];
2245:   }
2246:   /* Check each successive cone */
2247:   for (p = 1; p < numPoints; ++p) {
2248:     PetscInt newMeetSize = 0;

2250:     PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2251:     PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2252:     for (c = 0; c < dof; ++c) {
2253:       const PetscInt point = mesh->cones[off+c];

2255:       for (m = 0; m < meetSize; ++m) {
2256:         if (point == meet[i][m]) {
2257:           meet[1-i][newMeetSize++] = point;
2258:           break;
2259:         }
2260:       }
2261:     }
2262:     meetSize = newMeetSize;
2263:     i        = 1-i;
2264:   }
2265:   *numCoveringPoints = meetSize;
2266:   *coveringPoints    = meet[i];
2267:   DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2268:   return(0);
2269: }

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

2276:   Not Collective

2278:   Input Parameters:
2279: + dm - The DMPlex object
2280: . numPoints - The number of input points for the meet
2281: - points - The input points

2283:   Output Parameters:
2284: + numCoveredPoints - The number of points in the meet
2285: - coveredPoints - The points in the meet

2287:   Level: intermediate

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

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

2295: .keywords: mesh
2296: .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin()
2297: @*/
2298: PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2299: {

2307:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2308:   if (numCoveredPoints) *numCoveredPoints = 0;
2309:   return(0);
2310: }

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

2317:   Not Collective

2319:   Input Parameters:
2320: + dm - The DMPlex object
2321: . numPoints - The number of input points for the meet
2322: - points - The input points

2324:   Output Parameters:
2325: + numCoveredPoints - The number of points in the meet
2326: - coveredPoints - The points in the meet

2328:   Level: intermediate

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

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

2336: .keywords: mesh
2337: .seealso: DMPlexGetMeet(), DMPlexRestoreMeet(), DMPlexGetJoin()
2338: @*/
2339: PetscErrorCode DMPlexGetFullMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2340: {
2341:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2342:   PetscInt      *offsets, **closures;
2343:   PetscInt      *meet[2];
2344:   PetscInt       height = 0, maxSize, meetSize = 0, i = 0;
2345:   PetscInt       p, h, c, m, mc;


2354:   DMPlexGetDepth(dm, &height);
2355:   PetscMalloc1(numPoints, &closures);
2356:   DMGetWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2357:   mc      = mesh->maxConeSize;
2358:   maxSize = (mc > 1) ? ((PetscPowInt(mc,height+1)-1)/(mc-1)) : height + 1;
2359:   DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[0]);
2360:   DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[1]);

2362:   for (p = 0; p < numPoints; ++p) {
2363:     PetscInt closureSize;

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

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

2371:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
2372:       for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) {
2373:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2374:           offsets[p*(height+2)+h+1] = i;
2375:           break;
2376:         }
2377:       }
2378:       if (i == closureSize) offsets[p*(height+2)+h+1] = i;
2379:     }
2380:     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);
2381:   }
2382:   for (h = 0; h < height+1; ++h) {
2383:     PetscInt dof;

2385:     /* Copy in cone of first point */
2386:     dof = offsets[h+1] - offsets[h];
2387:     for (meetSize = 0; meetSize < dof; ++meetSize) {
2388:       meet[i][meetSize] = closures[0][(offsets[h]+meetSize)*2];
2389:     }
2390:     /* Check each successive cone */
2391:     for (p = 1; p < numPoints && meetSize; ++p) {
2392:       PetscInt newMeetSize = 0;

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

2398:         for (m = 0; m < meetSize; ++m) {
2399:           if (point == meet[i][m]) {
2400:             meet[1-i][newMeetSize++] = point;
2401:             break;
2402:           }
2403:         }
2404:       }
2405:       meetSize = newMeetSize;
2406:       i        = 1-i;
2407:     }
2408:     if (meetSize) break;
2409:   }
2410:   *numCoveredPoints = meetSize;
2411:   *coveredPoints    = meet[i];
2412:   for (p = 0; p < numPoints; ++p) {
2413:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, NULL, &closures[p]);
2414:   }
2415:   PetscFree(closures);
2416:   DMRestoreWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2417:   DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2418:   return(0);
2419: }

2423: /*@C
2424:   DMPlexEqual - Determine if two DMs have the same topology

2426:   Not Collective

2428:   Input Parameters:
2429: + dmA - A DMPlex object
2430: - dmB - A DMPlex object

2432:   Output Parameters:
2433: . equal - PETSC_TRUE if the topologies are identical

2435:   Level: intermediate

2437:   Notes:
2438:   We are not solving graph isomorphism, so we do not permutation.

2440: .keywords: mesh
2441: .seealso: DMPlexGetCone()
2442: @*/
2443: PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal)
2444: {
2445:   PetscInt       depth, depthB, pStart, pEnd, pStartB, pEndB, p;


2453:   *equal = PETSC_FALSE;
2454:   DMPlexGetDepth(dmA, &depth);
2455:   DMPlexGetDepth(dmB, &depthB);
2456:   if (depth != depthB) return(0);
2457:   DMPlexGetChart(dmA, &pStart,  &pEnd);
2458:   DMPlexGetChart(dmB, &pStartB, &pEndB);
2459:   if ((pStart != pStartB) || (pEnd != pEndB)) return(0);
2460:   for (p = pStart; p < pEnd; ++p) {
2461:     const PetscInt *cone, *coneB, *ornt, *orntB, *support, *supportB;
2462:     PetscInt        coneSize, coneSizeB, c, supportSize, supportSizeB, s;

2464:     DMPlexGetConeSize(dmA, p, &coneSize);
2465:     DMPlexGetCone(dmA, p, &cone);
2466:     DMPlexGetConeOrientation(dmA, p, &ornt);
2467:     DMPlexGetConeSize(dmB, p, &coneSizeB);
2468:     DMPlexGetCone(dmB, p, &coneB);
2469:     DMPlexGetConeOrientation(dmB, p, &orntB);
2470:     if (coneSize != coneSizeB) return(0);
2471:     for (c = 0; c < coneSize; ++c) {
2472:       if (cone[c] != coneB[c]) return(0);
2473:       if (ornt[c] != orntB[c]) return(0);
2474:     }
2475:     DMPlexGetSupportSize(dmA, p, &supportSize);
2476:     DMPlexGetSupport(dmA, p, &support);
2477:     DMPlexGetSupportSize(dmB, p, &supportSizeB);
2478:     DMPlexGetSupport(dmB, p, &supportB);
2479:     if (supportSize != supportSizeB) return(0);
2480:     for (s = 0; s < supportSize; ++s) {
2481:       if (support[s] != supportB[s]) return(0);
2482:     }
2483:   }
2484:   *equal = PETSC_TRUE;
2485:   return(0);
2486: }

2490: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2491: {
2492:   MPI_Comm       comm;

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

2561: /*@
2562:   DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point

2564:   Not Collective

2566:   Input Parameter:
2567: . dm    - The DMPlex object

2569:   Output Parameter:
2570: . depthLabel - The DMLabel recording point depth

2572:   Level: developer

2574: .keywords: mesh, points
2575: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2576: @*/
2577: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
2578: {

2584:   if (!dm->depthLabel) {DMGetLabel(dm, "depth", &dm->depthLabel);}
2585:   *depthLabel = dm->depthLabel;
2586:   return(0);
2587: }

2591: /*@
2592:   DMPlexGetDepth - Get the depth of the DAG representing this mesh

2594:   Not Collective

2596:   Input Parameter:
2597: . dm    - The DMPlex object

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

2602:   Level: developer

2604: .keywords: mesh, points
2605: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2606: @*/
2607: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
2608: {
2609:   DMLabel        label;
2610:   PetscInt       d = 0;

2616:   DMPlexGetDepthLabel(dm, &label);
2617:   if (label) {DMLabelGetNumValues(label, &d);}
2618:   *depth = d-1;
2619:   return(0);
2620: }

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

2627:   Not Collective

2629:   Input Parameters:
2630: + dm           - The DMPlex object
2631: - stratumValue - The requested depth

2633:   Output Parameters:
2634: + start - The first point at this depth
2635: - end   - One beyond the last point at this depth

2637:   Level: developer

2639: .keywords: mesh, points
2640: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
2641: @*/
2642: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2643: {
2644:   DMLabel        label;
2645:   PetscInt       pStart, pEnd;

2652:   DMPlexGetChart(dm, &pStart, &pEnd);
2653:   if (pStart == pEnd) return(0);
2654:   if (stratumValue < 0) {
2655:     if (start) *start = pStart;
2656:     if (end)   *end   = pEnd;
2657:     return(0);
2658:   }
2659:   DMPlexGetDepthLabel(dm, &label);
2660:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2661:   DMLabelGetStratumBounds(label, stratumValue, start, end);
2662:   return(0);
2663: }

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

2670:   Not Collective

2672:   Input Parameters:
2673: + dm           - The DMPlex object
2674: - stratumValue - The requested height

2676:   Output Parameters:
2677: + start - The first point at this height
2678: - end   - One beyond the last point at this height

2680:   Level: developer

2682: .keywords: mesh, points
2683: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
2684: @*/
2685: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2686: {
2687:   DMLabel        label;
2688:   PetscInt       depth, pStart, pEnd;

2695:   DMPlexGetChart(dm, &pStart, &pEnd);
2696:   if (pStart == pEnd) return(0);
2697:   if (stratumValue < 0) {
2698:     if (start) *start = pStart;
2699:     if (end)   *end   = pEnd;
2700:     return(0);
2701:   }
2702:   DMPlexGetDepthLabel(dm, &label);
2703:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2704:   DMLabelGetNumValues(label, &depth);
2705:   DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
2706:   return(0);
2707: }

2711: /* Set the number of dof on each point and separate by fields */
2712: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
2713: {
2714:   PetscInt      *pMax;
2715:   PetscInt       depth, pStart = 0, pEnd = 0;
2716:   PetscInt       Nf, p, d, dep, f;
2717:   PetscBool     *isFE;

2721:   PetscMalloc1(numFields, &isFE);
2722:   DMGetNumFields(dm, &Nf);
2723:   for (f = 0; f < numFields; ++f) {
2724:     PetscObject  obj;
2725:     PetscClassId id;

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

2755:       for (f = 0; f < numFields; ++f) {
2756:         if (isFE[f] && p >= pMax[dep]) continue;
2757:         PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
2758:         tot += numDof[f*(dim+1)+d];
2759:       }
2760:       PetscSectionSetDof(*section, p, tot);
2761:     }
2762:   }
2763:   PetscFree(pMax);
2764:   PetscFree(isFE);
2765:   return(0);
2766: }

2770: /* Set the number of dof on each point and separate by fields
2771:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
2772: */
2773: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2774: {
2775:   PetscInt       numFields;
2776:   PetscInt       bc;
2777:   PetscSection   aSec;

2781:   PetscSectionGetNumFields(section, &numFields);
2782:   for (bc = 0; bc < numBC; ++bc) {
2783:     PetscInt        field = 0;
2784:     const PetscInt *comp;
2785:     const PetscInt *idx;
2786:     PetscInt        Nc = -1, n, i;

2788:     if (numFields) field = bcField[bc];
2789:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2790:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2791:     ISGetLocalSize(bcPoints[bc], &n);
2792:     ISGetIndices(bcPoints[bc], &idx);
2793:     for (i = 0; i < n; ++i) {
2794:       const PetscInt p = idx[i];
2795:       PetscInt       numConst;

2797:       if (numFields) {
2798:         PetscSectionGetFieldDof(section, p, field, &numConst);
2799:       } else {
2800:         PetscSectionGetDof(section, p, &numConst);
2801:       }
2802:       /* If Nc < 0, constrain every dof on the point */
2803:       if (Nc > 0) numConst = PetscMin(numConst, Nc);
2804:       if (numFields) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
2805:       PetscSectionAddConstraintDof(section, p, numConst);
2806:     }
2807:     ISRestoreIndices(bcPoints[bc], &idx);
2808:     if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2809:   }
2810:   DMPlexGetAnchors(dm, &aSec, NULL);
2811:   if (aSec) {
2812:     PetscInt aStart, aEnd, a;

2814:     PetscSectionGetChart(aSec, &aStart, &aEnd);
2815:     for (a = aStart; a < aEnd; a++) {
2816:       PetscInt dof, f;

2818:       PetscSectionGetDof(aSec, a, &dof);
2819:       if (dof) {
2820:         /* if there are point-to-point constraints, then all dofs are constrained */
2821:         PetscSectionGetDof(section, a, &dof);
2822:         PetscSectionSetConstraintDof(section, a, dof);
2823:         for (f = 0; f < numFields; f++) {
2824:           PetscSectionGetFieldDof(section, a, f, &dof);
2825:           PetscSectionSetFieldConstraintDof(section, a, f, dof);
2826:         }
2827:       }
2828:     }
2829:   }
2830:   return(0);
2831: }

2835: /* Set the constrained field indices on each point
2836:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
2837: */
2838: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2839: {
2840:   PetscSection   aSec;
2841:   PetscInt      *indices;
2842:   PetscInt       numFields, maxDof, pStart, pEnd, p, bc, f, d;

2846:   PetscSectionGetNumFields(section, &numFields);
2847:   if (!numFields) return(0);
2848:   /* Initialize all field indices to -1 */
2849:   PetscSectionGetChart(section, &pStart, &pEnd);
2850:   PetscSectionGetMaxDof(section, &maxDof);
2851:   PetscMalloc1(maxDof, &indices);
2852:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
2853:   for (p = pStart; p < pEnd; ++p) for (f = 0; f < numFields; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
2854:   /* Handle BC constraints */
2855:   for (bc = 0; bc < numBC; ++bc) {
2856:     const PetscInt  field = bcField[bc];
2857:     const PetscInt *comp, *idx;
2858:     PetscInt        Nc = -1, n, i;

2860:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2861:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2862:     ISGetLocalSize(bcPoints[bc], &n);
2863:     ISGetIndices(bcPoints[bc], &idx);
2864:     for (i = 0; i < n; ++i) {
2865:       const PetscInt  p = idx[i];
2866:       const PetscInt *find;
2867:       PetscInt        fcdof, c;

2869:       PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
2870:       if (Nc < 0) {
2871:         for (d = 0; d < fcdof; ++d) indices[d] = d;
2872:       } else {
2873:         PetscSectionGetFieldConstraintIndices(section, p, field, &find);
2874:         for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
2875:         for (c = 0; c < Nc; ++c) indices[d+c] = comp[c];
2876:         PetscSortInt(d+Nc, indices);
2877:         for (c = d+Nc; c < fcdof; ++c) indices[c] = -1;
2878:       }
2879:       PetscSectionSetFieldConstraintIndices(section, p, field, indices);
2880:     }
2881:     if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2882:     ISRestoreIndices(bcPoints[bc], &idx);
2883:   }
2884:   /* Handle anchors */
2885:   DMPlexGetAnchors(dm, &aSec, NULL);
2886:   if (aSec) {
2887:     PetscInt aStart, aEnd, a;

2889:     for (d = 0; d < maxDof; ++d) indices[d] = d;
2890:     PetscSectionGetChart(aSec, &aStart, &aEnd);
2891:     for (a = aStart; a < aEnd; a++) {
2892:       PetscInt dof, fdof, f;

2894:       PetscSectionGetDof(aSec, a, &dof);
2895:       if (dof) {
2896:         /* if there are point-to-point constraints, then all dofs are constrained */
2897:         for (f = 0; f < numFields; f++) {
2898:           PetscSectionGetFieldDof(section, a, f, &fdof);
2899:           PetscSectionSetFieldConstraintIndices(section, a, f, indices);
2900:         }
2901:       }
2902:     }
2903:   }
2904:   PetscFree(indices);
2905:   return(0);
2906: }

2910: /* Set the constrained indices on each point */
2911: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
2912: {
2913:   PetscInt      *indices;
2914:   PetscInt       numFields, maxDof, pStart, pEnd, p, f, d;

2918:   PetscSectionGetNumFields(section, &numFields);
2919:   PetscSectionGetMaxDof(section, &maxDof);
2920:   PetscSectionGetChart(section, &pStart, &pEnd);
2921:   PetscMalloc1(maxDof, &indices);
2922:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
2923:   for (p = pStart; p < pEnd; ++p) {
2924:     PetscInt cdof, d;

2926:     PetscSectionGetConstraintDof(section, p, &cdof);
2927:     if (cdof) {
2928:       if (numFields) {
2929:         PetscInt numConst = 0, foff = 0;

2931:         for (f = 0; f < numFields; ++f) {
2932:           const PetscInt *find;
2933:           PetscInt        fcdof, fdof;

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

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

2959:   Not Collective

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

2973:   Output Parameter:
2974: . section - The PetscSection object

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

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

2981:   Level: developer

2983:   Fortran Notes:
2984:   A Fortran 90 version is available as DMPlexCreateSectionF90()

2986: .keywords: mesh, elements
2987: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
2988: @*/
2989: 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)
2990: {
2991:   PetscSection   aSec;

2995:   DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
2996:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
2997:   if (perm) {PetscSectionSetPermutation(*section, perm);}
2998:   PetscSectionSetUp(*section);
2999:   DMPlexGetAnchors(dm,&aSec,NULL);
3000:   if (numBC || aSec) {
3001:     DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
3002:     DMPlexCreateSectionBCIndices(dm, *section);
3003:   }
3004:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
3005:   return(0);
3006: }

3010: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
3011: {
3012:   PetscSection   section, s;
3013:   Mat            m;
3014:   PetscInt       maxHeight;

3018:   DMClone(dm, cdm);
3019:   DMPlexGetMaxProjectionHeight(dm, &maxHeight);
3020:   DMPlexSetMaxProjectionHeight(*cdm, maxHeight);
3021:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
3022:   DMSetDefaultSection(*cdm, section);
3023:   PetscSectionDestroy(&section);
3024:   PetscSectionCreate(PETSC_COMM_SELF, &s);
3025:   MatCreate(PETSC_COMM_SELF, &m);
3026:   DMSetDefaultConstraints(*cdm, s, m);
3027:   PetscSectionDestroy(&s);
3028:   MatDestroy(&m);
3029:   return(0);
3030: }

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

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

3046: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
3047: {
3048:   DM_Plex *mesh = (DM_Plex*) dm->data;

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

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

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

3070: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
3071: {
3072:   DM_Plex *mesh = (DM_Plex*) dm->data;

3076:   if (coneOrientations) *coneOrientations = mesh->coneOrientations;
3077:   return(0);
3078: }

3080: /******************************** FEM Support **********************************/

3084: PetscErrorCode DMPlexCreateSpectralClosurePermutation(DM dm, PetscSection section)
3085: {
3086:   PetscInt      *perm;
3087:   PetscInt       dim, eStart, k, Nf, f, Nc, c, i, size = 0, offset = 0, foffset = 0;

3091:   if (!section) {DMGetDefaultSection(dm, &section);}
3092:   DMGetDimension(dm, &dim);
3093:   PetscSectionGetNumFields(section, &Nf);
3094:   if (dim != 2) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Can only do quads now");
3095:   for (f = 0; f < Nf; ++f) {
3096:     /* An order k SEM disc has k-1 dofs on an edge */
3097:     DMPlexGetDepthStratum(dm, 1, &eStart, NULL);
3098:     PetscSectionGetFieldDof(section, eStart, f, &k);
3099:     PetscSectionGetFieldComponents(section, f, &Nc);
3100:     k = k/Nc + 1;
3101:     size += PetscSqr(k+1)*Nc;
3102:   }
3103:   PetscMalloc1(size, &perm);
3104:   for (f = 0; f < Nf; ++f) {
3105:     /* The original quad closure is oriented clockwise, {f, e_b, e_r, e_t, e_l, v_lb, v_rb, v_tr, v_tl} */
3106:     DMPlexGetDepthStratum(dm, 1, &eStart, NULL);
3107:     PetscSectionGetFieldDof(section, eStart, f, &k);
3108:     PetscSectionGetFieldComponents(section, f, &Nc);
3109:     k = k/Nc + 1;
3110:     /* The SEM order is

3112:        v_lb, {e_b}, v_rb,
3113:        e^{(k-1)-i}_l, {f^{((k-1)-i)*(k-1)}}, e^i_r,
3114:        v_lt, reverse {e_t}, v_rt
3115:     */
3116:     {
3117:       const PetscInt of   = 0;
3118:       const PetscInt oeb  = of   + PetscSqr(k-1);
3119:       const PetscInt oer  = oeb  + (k-1);
3120:       const PetscInt oet  = oer  + (k-1);
3121:       const PetscInt oel  = oet  + (k-1);
3122:       const PetscInt ovlb = oel  + (k-1);
3123:       const PetscInt ovrb = ovlb + 1;
3124:       const PetscInt ovrt = ovrb + 1;
3125:       const PetscInt ovlt = ovrt + 1;
3126:       PetscInt       o;

3128:       /* bottom */
3129:       for (c = 0; c < Nc; ++c, ++offset) perm[offset] = ovlb*Nc + c + foffset;
3130:       for (o = oeb; o < oer; ++o) for (c = 0; c < Nc; ++c, ++offset) perm[offset] = o*Nc + c + foffset;
3131:       for (c = 0; c < Nc; ++c, ++offset) perm[offset] = ovrb*Nc + c + foffset;
3132:       /* middle */
3133:       for (i = 0; i < k-1; ++i) {
3134:         for (c = 0; c < Nc; ++c, ++offset) perm[offset] = (oel+(k-2)-i)*Nc + c + foffset;
3135:         for (o = of+(k-1)*i; o < of+(k-1)*(i+1); ++o) for (c = 0; c < Nc; ++c, ++offset) perm[offset] = o*Nc + c + foffset;
3136:         for (c = 0; c < Nc; ++c, ++offset) perm[offset] = (oer+i)*Nc + c + foffset;
3137:       }
3138:       /* top */
3139:       for (c = 0; c < Nc; ++c, ++offset) perm[offset] = ovlt*Nc + c + foffset;
3140:       for (o = oel-1; o >= oet; --o) for (c = 0; c < Nc; ++c, ++offset) perm[offset] = o*Nc + c + foffset;
3141:       for (c = 0; c < Nc; ++c, ++offset) perm[offset] = ovrt*Nc + c + foffset;
3142:       foffset = offset;
3143:     }
3144:   }
3145:   /* Check permutation */
3146:   {
3147:     PetscInt *check;

3149:     PetscMalloc1(size, &check);
3150:     for (i = 0; i < size; ++i) {check[i] = -1; if (perm[i] < 0 || perm[i] >= size) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid permutation index p[%D] = %D", i, perm[i]);}
3151:     for (i = 0; i < size; ++i) check[perm[i]] = i;
3152:     for (i = 0; i < size; ++i) {if (check[i] < 0) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Missing permutation index %D", i);}
3153:     PetscFree(check);
3154:   }
3155:   PetscSectionSetClosurePermutation_Internal(section, (PetscObject) dm, size, PETSC_OWN_POINTER, perm);
3156:   return(0);
3157: }

3161: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3162: {
3163:   PetscScalar    *array, *vArray;
3164:   const PetscInt *cone, *coneO;
3165:   PetscInt        pStart, pEnd, p, numPoints, size = 0, offset = 0;
3166:   PetscErrorCode  ierr;

3169:   PetscSectionGetChart(section, &pStart, &pEnd);
3170:   DMPlexGetConeSize(dm, point, &numPoints);
3171:   DMPlexGetCone(dm, point, &cone);
3172:   DMPlexGetConeOrientation(dm, point, &coneO);
3173:   if (!values || !*values) {
3174:     if ((point >= pStart) && (point < pEnd)) {
3175:       PetscInt dof;

3177:       PetscSectionGetDof(section, point, &dof);
3178:       size += dof;
3179:     }
3180:     for (p = 0; p < numPoints; ++p) {
3181:       const PetscInt cp = cone[p];
3182:       PetscInt       dof;

3184:       if ((cp < pStart) || (cp >= pEnd)) continue;
3185:       PetscSectionGetDof(section, cp, &dof);
3186:       size += dof;
3187:     }
3188:     if (!values) {
3189:       if (csize) *csize = size;
3190:       return(0);
3191:     }
3192:     DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
3193:   } else {
3194:     array = *values;
3195:   }
3196:   size = 0;
3197:   VecGetArray(v, &vArray);
3198:   if ((point >= pStart) && (point < pEnd)) {
3199:     PetscInt     dof, off, d;
3200:     PetscScalar *varr;

3202:     PetscSectionGetDof(section, point, &dof);
3203:     PetscSectionGetOffset(section, point, &off);
3204:     varr = &vArray[off];
3205:     for (d = 0; d < dof; ++d, ++offset) {
3206:       array[offset] = varr[d];
3207:     }
3208:     size += dof;
3209:   }
3210:   for (p = 0; p < numPoints; ++p) {
3211:     const PetscInt cp = cone[p];
3212:     PetscInt       o  = coneO[p];
3213:     PetscInt       dof, off, d;
3214:     PetscScalar   *varr;

3216:     if ((cp < pStart) || (cp >= pEnd)) continue;
3217:     PetscSectionGetDof(section, cp, &dof);
3218:     PetscSectionGetOffset(section, cp, &off);
3219:     varr = &vArray[off];
3220:     if (o >= 0) {
3221:       for (d = 0; d < dof; ++d, ++offset) {
3222:         array[offset] = varr[d];
3223:       }
3224:     } else {
3225:       for (d = dof-1; d >= 0; --d, ++offset) {
3226:         array[offset] = varr[d];
3227:       }
3228:     }
3229:     size += dof;
3230:   }
3231:   VecRestoreArray(v, &vArray);
3232:   if (!*values) {
3233:     if (csize) *csize = size;
3234:     *values = array;
3235:   } else {
3236:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3237:     *csize = size;
3238:   }
3239:   return(0);
3240: }

3244: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], const PetscInt perm[], const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3245: {
3246:   PetscInt       offset = 0, p;

3250:   *size = 0;
3251:   for (p = 0; p < numPoints*2; p += 2) {
3252:     const PetscInt point = points[p];
3253:     const PetscInt o     = points[p+1];
3254:     PetscInt       dof, off, d;
3255:     const PetscScalar *varr;

3257:     PetscSectionGetDof(section, point, &dof);
3258:     PetscSectionGetOffset(section, point, &off);
3259:     varr = &vArray[off];
3260:     if (o >= 0) {
3261:       for (d = 0; d < dof; ++d, ++offset)    array[perm ? perm[offset] : offset] = varr[d];
3262:     } else {
3263:       for (d = dof-1; d >= 0; --d, ++offset) array[perm ? perm[offset] : offset] = varr[d];
3264:     }
3265:   }
3266:   *size = offset;
3267:   return(0);
3268: }

3272: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Fields_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], PetscInt numFields, const PetscInt perm[], const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3273: {
3274:   PetscInt       offset = 0, f;

3278:   *size = 0;
3279:   for (f = 0; f < numFields; ++f) {
3280:     PetscInt fcomp, p;

3282:     PetscSectionGetFieldComponents(section, f, &fcomp);
3283:     for (p = 0; p < numPoints*2; p += 2) {
3284:       const PetscInt point = points[p];
3285:       const PetscInt o     = points[p+1];
3286:       PetscInt       fdof, foff, d, c;
3287:       const PetscScalar *varr;

3289:       PetscSectionGetFieldDof(section, point, f, &fdof);
3290:       PetscSectionGetFieldOffset(section, point, f, &foff);
3291:       varr = &vArray[foff];
3292:       if (o >= 0) {
3293:         for (d = 0; d < fdof; ++d, ++offset) array[perm ? perm[offset] : offset] = varr[d];
3294:       } else {
3295:         for (d = fdof/fcomp-1; d >= 0; --d) {
3296:           for (c = 0; c < fcomp; ++c, ++offset) {
3297:             array[perm ? perm[offset] : offset] = varr[d*fcomp+c];
3298:           }
3299:         }
3300:       }
3301:     }
3302:   }
3303:   *size = offset;
3304:   return(0);
3305: }

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

3312:   Not collective

3314:   Input Parameters:
3315: + dm - The DM
3316: . section - The section describing the layout in v, or NULL to use the default section
3317: . v - The local vector
3318: - point - The sieve point in the DM

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

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

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

3330:   Level: intermediate

3332: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3333: @*/
3334: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3335: {
3336:   PetscSection    clSection;
3337:   IS              clPoints;
3338:   PetscScalar    *array, *vArray;
3339:   PetscInt       *points = NULL;
3340:   const PetscInt *clp, *perm;
3341:   PetscInt        depth, numFields, numPoints, size;
3342:   PetscErrorCode  ierr;

3346:   if (!section) {DMGetDefaultSection(dm, &section);}
3349:   DMPlexGetDepth(dm, &depth);
3350:   PetscSectionGetNumFields(section, &numFields);
3351:   if (depth == 1 && numFields < 2) {
3352:     DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
3353:     return(0);
3354:   }
3355:   /* Get points */
3356:   PetscSectionGetClosureInversePermutation_Internal(section, (PetscObject) dm, NULL, &perm);
3357:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3358:   if (!clPoints) {
3359:     PetscInt pStart, pEnd, p, q;

3361:     PetscSectionGetChart(section, &pStart, &pEnd);
3362:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3363:     /* Compress out points not in the section */
3364:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3365:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3366:         points[q*2]   = points[p];
3367:         points[q*2+1] = points[p+1];
3368:         ++q;
3369:       }
3370:     }
3371:     numPoints = q;
3372:   } else {
3373:     PetscInt dof, off;

3375:     PetscSectionGetDof(clSection, point, &dof);
3376:     PetscSectionGetOffset(clSection, point, &off);
3377:     ISGetIndices(clPoints, &clp);
3378:     numPoints = dof/2;
3379:     points    = (PetscInt *) &clp[off];
3380:   }
3381:   /* Get array */
3382:   if (!values || !*values) {
3383:     PetscInt asize = 0, dof, p;

3385:     for (p = 0; p < numPoints*2; p += 2) {
3386:       PetscSectionGetDof(section, points[p], &dof);
3387:       asize += dof;
3388:     }
3389:     if (!values) {
3390:       if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3391:       else           {ISRestoreIndices(clPoints, &clp);}
3392:       if (csize) *csize = asize;
3393:       return(0);
3394:     }
3395:     DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
3396:   } else {
3397:     array = *values;
3398:   }
3399:   VecGetArray(v, &vArray);
3400:   /* Get values */
3401:   if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, perm, vArray, &size, array);}
3402:   else               {DMPlexVecGetClosure_Static(section, numPoints, points, perm, vArray, &size, array);}
3403:   /* Cleanup points */
3404:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3405:   else           {ISRestoreIndices(clPoints, &clp);}
3406:   /* Cleanup array */
3407:   VecRestoreArray(v, &vArray);
3408:   if (!*values) {
3409:     if (csize) *csize = size;
3410:     *values = array;
3411:   } else {
3412:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %D < actual size %D", *csize, size);
3413:     *csize = size;
3414:   }
3415:   return(0);
3416: }

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

3423:   Not collective

3425:   Input Parameters:
3426: + dm - The DM
3427: . section - The section describing the layout in v, or NULL to use the default section
3428: . v - The local vector
3429: . point - The sieve point in the DM
3430: . csize - The number of values in the closure, or NULL
3431: - values - The array of values, which is a borrowed array and should not be freed

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

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

3439:   Level: intermediate

3441: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3442: @*/
3443: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3444: {
3445:   PetscInt       size = 0;

3449:   /* Should work without recalculating size */
3450:   DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
3451:   return(0);
3452: }

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

3459: PETSC_STATIC_INLINE PetscErrorCode updatePoint_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, PetscInt orientation, const PetscInt perm[], const PetscScalar values[], PetscInt offset, PetscScalar array[])
3460: {
3461:   PetscInt        cdof;   /* The number of constraints on this point */
3462:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3463:   PetscScalar    *a;
3464:   PetscInt        off, cind = 0, k;
3465:   PetscErrorCode  ierr;

3468:   PetscSectionGetConstraintDof(section, point, &cdof);
3469:   PetscSectionGetOffset(section, point, &off);
3470:   a    = &array[off];
3471:   if (!cdof || setBC) {
3472:     if (orientation >= 0) {
3473:       for (k = 0; k < dof; ++k) {
3474:         fuse(&a[k], values[perm ? perm[offset+k] : offset+k]);
3475:       }
3476:     } else {
3477:       for (k = 0; k < dof; ++k) {
3478:         fuse(&a[k], values[perm ? perm[offset+dof-k-1] : offset+dof-k-1]);
3479:       }
3480:     }
3481:   } else {
3482:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3483:     if (orientation >= 0) {
3484:       for (k = 0; k < dof; ++k) {
3485:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3486:         fuse(&a[k], values[perm ? perm[offset+k] : offset+k]);
3487:       }
3488:     } else {
3489:       for (k = 0; k < dof; ++k) {
3490:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3491:         fuse(&a[k], values[perm ? perm[offset+dof-k-1] : offset+dof-k-1]);
3492:       }
3493:     }
3494:   }
3495:   return(0);
3496: }

3500: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscInt perm[], const PetscScalar values[], PetscInt offset, PetscScalar array[])
3501: {
3502:   PetscInt        cdof;   /* The number of constraints on this point */
3503:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3504:   PetscScalar    *a;
3505:   PetscInt        off, cind = 0, k;
3506:   PetscErrorCode  ierr;

3509:   PetscSectionGetConstraintDof(section, point, &cdof);
3510:   PetscSectionGetOffset(section, point, &off);
3511:   a    = &array[off];
3512:   if (cdof) {
3513:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3514:     if (orientation >= 0) {
3515:       for (k = 0; k < dof; ++k) {
3516:         if ((cind < cdof) && (k == cdofs[cind])) {
3517:           fuse(&a[k], values[perm ? perm[offset+k] : offset+k]);
3518:           ++cind;
3519:         }
3520:       }
3521:     } else {
3522:       for (k = 0; k < dof; ++k) {
3523:         if ((cind < cdof) && (k == cdofs[cind])) {
3524:           fuse(&a[k], values[perm ? perm[offset+dof-k-1] : offset+dof-k-1]);
3525:           ++cind;
3526:         }
3527:       }
3528:     }
3529:   }
3530:   return(0);
3531: }

3535: PETSC_STATIC_INLINE PetscErrorCode updatePointFields_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, const PetscInt perm[], const PetscScalar values[], PetscInt *offset, PetscScalar array[])
3536: {
3537:   PetscScalar    *a;
3538:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3539:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3540:   PetscInt        cind = 0, k, c;
3541:   PetscErrorCode  ierr;

3544:   PetscSectionGetFieldDof(section, point, f, &fdof);
3545:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3546:   PetscSectionGetFieldOffset(section, point, f, &foff);
3547:   a    = &array[foff];
3548:   if (!fcdof || setBC) {
3549:     if (o >= 0) {
3550:       for (k = 0; k < fdof; ++k) fuse(&a[k], values[perm ? perm[foffset+k] : foffset+k]);
3551:     } else {
3552:       for (k = fdof/fcomp-1; k >= 0; --k) {
3553:         for (c = 0; c < fcomp; ++c) {
3554:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[perm ? perm[foffset+k*fcomp+c] : foffset+k*fcomp+c]);
3555:         }
3556:       }
3557:     }
3558:   } else {
3559:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3560:     if (o >= 0) {
3561:       for (k = 0; k < fdof; ++k) {
3562:         if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
3563:         fuse(&a[k], values[perm ? perm[foffset+k] : foffset+k]);
3564:       }
3565:     } else {
3566:       for (k = fdof/fcomp-1; k >= 0; --k) {
3567:         for (c = 0; c < fcomp; ++c) {
3568:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
3569:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[perm ? perm[foffset+k*fcomp+c] : foffset+k*fcomp+c]);
3570:         }
3571:       }
3572:     }
3573:   }
3574:   *offset += fdof;
3575:   return(0);
3576: }

3580: PETSC_STATIC_INLINE PetscErrorCode updatePointFieldsBC_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), const PetscInt perm[], const PetscScalar values[], PetscInt *offset, PetscScalar array[])
3581: {
3582:   PetscScalar    *a;
3583:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3584:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3585:   PetscInt        cind = 0, k, c;
3586:   PetscErrorCode  ierr;

3589:   PetscSectionGetFieldDof(section, point, f, &fdof);
3590:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3591:   PetscSectionGetFieldOffset(section, point, f, &foff);
3592:   a    = &array[foff];
3593:   if (fcdof) {
3594:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3595:     if (o >= 0) {
3596:       for (k = 0; k < fdof; ++k) {
3597:         if ((cind < fcdof) && (k == fcdofs[cind])) {
3598:           fuse(&a[k], values[perm ? perm[foffset+k] : foffset+k]);
3599:           ++cind;
3600:         }
3601:       }
3602:     } else {
3603:       for (k = fdof/fcomp-1; k >= 0; --k) {
3604:         for (c = 0; c < fcomp; ++c) {
3605:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {
3606:             fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[perm ? perm[foffset+k*fcomp+c] : foffset+k*fcomp+c]);
3607:             ++cind;
3608:           }
3609:         }
3610:       }
3611:     }
3612:   }
3613:   *offset += fdof;
3614:   return(0);
3615: }

3619: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3620: {
3621:   PetscScalar    *array;
3622:   const PetscInt *cone, *coneO;
3623:   PetscInt        pStart, pEnd, p, numPoints, off, dof;
3624:   PetscErrorCode  ierr;

3627:   PetscSectionGetChart(section, &pStart, &pEnd);
3628:   DMPlexGetConeSize(dm, point, &numPoints);
3629:   DMPlexGetCone(dm, point, &cone);
3630:   DMPlexGetConeOrientation(dm, point, &coneO);
3631:   VecGetArray(v, &array);
3632:   for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
3633:     const PetscInt cp = !p ? point : cone[p-1];
3634:     const PetscInt o  = !p ? 0     : coneO[p-1];

3636:     if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
3637:     PetscSectionGetDof(section, cp, &dof);
3638:     /* ADD_VALUES */
3639:     {
3640:       const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3641:       PetscScalar    *a;
3642:       PetscInt        cdof, coff, cind = 0, k;

3644:       PetscSectionGetConstraintDof(section, cp, &cdof);
3645:       PetscSectionGetOffset(section, cp, &coff);
3646:       a    = &array[coff];
3647:       if (!cdof) {
3648:         if (o >= 0) {
3649:           for (k = 0; k < dof; ++k) {
3650:             a[k] += values[off+k];
3651:           }
3652:         } else {
3653:           for (k = 0; k < dof; ++k) {
3654:             a[k] += values[off+dof-k-1];
3655:           }
3656:         }
3657:       } else {
3658:         PetscSectionGetConstraintIndices(section, cp, &cdofs);
3659:         if (o >= 0) {
3660:           for (k = 0; k < dof; ++k) {
3661:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3662:             a[k] += values[off+k];
3663:           }
3664:         } else {
3665:           for (k = 0; k < dof; ++k) {
3666:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3667:             a[k] += values[off+dof-k-1];
3668:           }
3669:         }
3670:       }
3671:     }
3672:   }
3673:   VecRestoreArray(v, &array);
3674:   return(0);
3675: }

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

3682:   Not collective

3684:   Input Parameters:
3685: + dm - The DM
3686: . section - The section describing the layout in v, or NULL to use the default section
3687: . v - The local vector
3688: . point - The sieve point in the DM
3689: . values - The array of values
3690: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

3695:   Level: intermediate

3697: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
3698: @*/
3699: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3700: {
3701:   PetscSection    clSection;
3702:   IS              clPoints;
3703:   PetscScalar    *array;
3704:   PetscInt       *points = NULL;
3705:   const PetscInt *clp, *perm;
3706:   PetscInt        depth, numFields, numPoints, p;
3707:   PetscErrorCode  ierr;

3711:   if (!section) {DMGetDefaultSection(dm, &section);}
3714:   DMPlexGetDepth(dm, &depth);
3715:   PetscSectionGetNumFields(section, &numFields);
3716:   if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
3717:     DMPlexVecSetClosure_Depth1_Static(dm, section, v, point, values, mode);
3718:     return(0);
3719:   }
3720:   /* Get points */
3721:   PetscSectionGetClosureInversePermutation_Internal(section, (PetscObject) dm, NULL, &perm);
3722:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3723:   if (!clPoints) {
3724:     PetscInt pStart, pEnd, q;

3726:     PetscSectionGetChart(section, &pStart, &pEnd);
3727:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3728:     /* Compress out points not in the section */
3729:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3730:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3731:         points[q*2]   = points[p];
3732:         points[q*2+1] = points[p+1];
3733:         ++q;
3734:       }
3735:     }
3736:     numPoints = q;
3737:   } else {
3738:     PetscInt dof, off;

3740:     PetscSectionGetDof(clSection, point, &dof);
3741:     PetscSectionGetOffset(clSection, point, &off);
3742:     ISGetIndices(clPoints, &clp);
3743:     numPoints = dof/2;
3744:     points    = (PetscInt *) &clp[off];
3745:   }
3746:   /* Get array */
3747:   VecGetArray(v, &array);
3748:   /* Get values */
3749:   if (numFields > 0) {
3750:     PetscInt offset = 0, fcomp, f;
3751:     for (f = 0; f < numFields; ++f) {
3752:       PetscSectionGetFieldComponents(section, f, &fcomp);
3753:       switch (mode) {
3754:       case INSERT_VALUES:
3755:         for (p = 0; p < numPoints*2; p += 2) {
3756:           const PetscInt point = points[p];
3757:           const PetscInt o     = points[p+1];
3758:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, perm, values, &offset, array);
3759:         } break;
3760:       case INSERT_ALL_VALUES:
3761:         for (p = 0; p < numPoints*2; p += 2) {
3762:           const PetscInt point = points[p];
3763:           const PetscInt o     = points[p+1];
3764:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, perm, values, &offset, array);
3765:         } break;
3766:       case INSERT_BC_VALUES:
3767:         for (p = 0; p < numPoints*2; p += 2) {
3768:           const PetscInt point = points[p];
3769:           const PetscInt o     = points[p+1];
3770:           updatePointFieldsBC_private(section, point, o, f, fcomp, insert, perm, values, &offset, array);
3771:         } break;
3772:       case ADD_VALUES:
3773:         for (p = 0; p < numPoints*2; p += 2) {
3774:           const PetscInt point = points[p];
3775:           const PetscInt o     = points[p+1];
3776:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, perm, values, &offset, array);
3777:         } break;
3778:       case ADD_ALL_VALUES:
3779:         for (p = 0; p < numPoints*2; p += 2) {
3780:           const PetscInt point = points[p];
3781:           const PetscInt o     = points[p+1];
3782:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, perm, values, &offset, array);
3783:         } break;
3784:       case ADD_BC_VALUES:
3785:         for (p = 0; p < numPoints*2; p += 2) {
3786:           const PetscInt point = points[p];
3787:           const PetscInt o     = points[p+1];
3788:           updatePointFieldsBC_private(section, point, o, f, fcomp, add, perm, values, &offset, array);
3789:         } break;
3790:       default:
3791:         SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3792:       }
3793:     }
3794:   } else {
3795:     PetscInt dof, off;

3797:     switch (mode) {
3798:     case INSERT_VALUES:
3799:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3800:         PetscInt o = points[p+1];
3801:         PetscSectionGetDof(section, points[p], &dof);
3802:         updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, perm, values, off, array);
3803:       } break;
3804:     case INSERT_ALL_VALUES:
3805:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3806:         PetscInt o = points[p+1];
3807:         PetscSectionGetDof(section, points[p], &dof);
3808:         updatePoint_private(section, points[p], dof, insert, PETSC_TRUE,  o, perm, values, off, array);
3809:       } break;
3810:     case INSERT_BC_VALUES:
3811:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3812:         PetscInt o = points[p+1];
3813:         PetscSectionGetDof(section, points[p], &dof);
3814:         updatePointBC_private(section, points[p], dof, insert,  o, perm, values, off, array);
3815:       } break;
3816:     case ADD_VALUES:
3817:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3818:         PetscInt o = points[p+1];
3819:         PetscSectionGetDof(section, points[p], &dof);
3820:         updatePoint_private(section, points[p], dof, add,    PETSC_FALSE, o, perm, values, off, array);
3821:       } break;
3822:     case ADD_ALL_VALUES:
3823:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3824:         PetscInt o = points[p+1];
3825:         PetscSectionGetDof(section, points[p], &dof);
3826:         updatePoint_private(section, points[p], dof, add,    PETSC_TRUE,  o, perm, values, off, array);
3827:       } break;
3828:     case ADD_BC_VALUES:
3829:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3830:         PetscInt o = points[p+1];
3831:         PetscSectionGetDof(section, points[p], &dof);
3832:         updatePointBC_private(section, points[p], dof, add,  o, perm, values, off, array);
3833:       } break;
3834:     default:
3835:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3836:     }
3837:   }
3838:   /* Cleanup points */
3839:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3840:   else           {ISRestoreIndices(clPoints, &clp);}
3841:   /* Cleanup array */
3842:   VecRestoreArray(v, &array);
3843:   return(0);
3844: }

3848: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
3849: {
3850:   PetscSection    clSection;
3851:   IS              clPoints;
3852:   PetscScalar    *array;
3853:   PetscInt       *points = NULL;
3854:   const PetscInt *clp, *perm;
3855:   PetscInt        numFields, numPoints, p;
3856:   PetscInt        offset = 0, fcomp, f;
3857:   PetscErrorCode  ierr;

3861:   if (!section) {DMGetDefaultSection(dm, &section);}
3864:   PetscSectionGetNumFields(section, &numFields);
3865:   /* Get points */
3866:   PetscSectionGetClosureInversePermutation_Internal(section, (PetscObject) dm, NULL, &perm);
3867:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3868:   if (!clPoints) {
3869:     PetscInt pStart, pEnd, q;

3871:     PetscSectionGetChart(section, &pStart, &pEnd);
3872:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3873:     /* Compress out points not in the section */
3874:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3875:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3876:         points[q*2]   = points[p];
3877:         points[q*2+1] = points[p+1];
3878:         ++q;
3879:       }
3880:     }
3881:     numPoints = q;
3882:   } else {
3883:     PetscInt dof, off;

3885:     PetscSectionGetDof(clSection, point, &dof);
3886:     PetscSectionGetOffset(clSection, point, &off);
3887:     ISGetIndices(clPoints, &clp);
3888:     numPoints = dof/2;
3889:     points    = (PetscInt *) &clp[off];
3890:   }
3891:   /* Get array */
3892:   VecGetArray(v, &array);
3893:   /* Get values */
3894:   for (f = 0; f < numFields; ++f) {
3895:     PetscSectionGetFieldComponents(section, f, &fcomp);
3896:     if (!fieldActive[f]) {
3897:       for (p = 0; p < numPoints*2; p += 2) {
3898:         PetscInt fdof;
3899:         PetscSectionGetFieldDof(section, points[p], f, &fdof);
3900:         offset += fdof;
3901:       }
3902:       continue;
3903:     }
3904:     switch (mode) {
3905:     case INSERT_VALUES:
3906:       for (p = 0; p < numPoints*2; p += 2) {
3907:         const PetscInt point = points[p];
3908:         const PetscInt o     = points[p+1];
3909:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, perm, values, &offset, array);
3910:       } break;
3911:     case INSERT_ALL_VALUES:
3912:       for (p = 0; p < numPoints*2; p += 2) {
3913:         const PetscInt point = points[p];
3914:         const PetscInt o     = points[p+1];
3915:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, perm, values, &offset, array);
3916:         } break;
3917:     case INSERT_BC_VALUES:
3918:       for (p = 0; p < numPoints*2; p += 2) {
3919:         const PetscInt point = points[p];
3920:         const PetscInt o     = points[p+1];
3921:         updatePointFieldsBC_private(section, point, o, f, fcomp, insert, perm, values, &offset, array);
3922:       } break;
3923:     case ADD_VALUES:
3924:       for (p = 0; p < numPoints*2; p += 2) {
3925:         const PetscInt point = points[p];
3926:         const PetscInt o     = points[p+1];
3927:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, perm, values, &offset, array);
3928:       } break;
3929:     case ADD_ALL_VALUES:
3930:       for (p = 0; p < numPoints*2; p += 2) {
3931:         const PetscInt point = points[p];
3932:         const PetscInt o     = points[p+1];
3933:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, perm, values, &offset, array);
3934:       } break;
3935:     default:
3936:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3937:     }
3938:   }
3939:   /* Cleanup points */
3940:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3941:   else           {ISRestoreIndices(clPoints, &clp);}
3942:   /* Cleanup array */
3943:   VecRestoreArray(v, &array);
3944:   return(0);
3945: }

3949: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
3950: {
3951:   PetscMPIInt    rank;
3952:   PetscInt       i, j;

3956:   MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
3957:   PetscViewerASCIIPrintf(viewer, "[%d]mat for sieve point %D\n", rank, point);
3958:   for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
3959:   for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
3960:   numCIndices = numCIndices ? numCIndices : numRIndices;
3961:   for (i = 0; i < numRIndices; i++) {
3962:     PetscViewerASCIIPrintf(viewer, "[%d]", rank);
3963:     for (j = 0; j < numCIndices; j++) {
3964: #if defined(PETSC_USE_COMPLEX)
3965:       PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
3966: #else
3967:       PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
3968: #endif
3969:     }
3970:     PetscViewerASCIIPrintf(viewer, "\n");
3971:   }
3972:   return(0);
3973: }

3977: /* . off - The global offset of this point */
3978: PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
3979: {
3980:   PetscInt        dof;    /* The number of unknowns on this point */
3981:   PetscInt        cdof;   /* The number of constraints on this point */
3982:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3983:   PetscInt        cind = 0, k;
3984:   PetscErrorCode  ierr;

3987:   PetscSectionGetDof(section, point, &dof);
3988:   PetscSectionGetConstraintDof(section, point, &cdof);
3989:   if (!cdof || setBC) {
3990:     if (orientation >= 0) {
3991:       for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
3992:     } else {
3993:       for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
3994:     }
3995:   } else {
3996:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3997:     if (orientation >= 0) {
3998:       for (k = 0; k < dof; ++k) {
3999:         if ((cind < cdof) && (k == cdofs[cind])) {
4000:           /* Insert check for returning constrained indices */
4001:           indices[*loff+k] = -(off+k+1);
4002:           ++cind;
4003:         } else {
4004:           indices[*loff+k] = off+k-cind;
4005:         }
4006:       }
4007:     } else {
4008:       for (k = 0; k < dof; ++k) {
4009:         if ((cind < cdof) && (k == cdofs[cind])) {
4010:           /* Insert check for returning constrained indices */
4011:           indices[*loff+dof-k-1] = -(off+k+1);
4012:           ++cind;
4013:         } else {
4014:           indices[*loff+dof-k-1] = off+k-cind;
4015:         }
4016:       }
4017:     }
4018:   }
4019:   *loff += dof;
4020:   return(0);
4021: }

4025: /* . off - The global offset of this point */
4026: PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection section, PetscInt point, PetscInt off, PetscInt foffs[], PetscBool setBC, PetscInt orientation, PetscInt indices[])
4027: {
4028:   PetscInt       numFields, foff, f;

4032:   PetscSectionGetNumFields(section, &numFields);
4033:   for (f = 0, foff = 0; f < numFields; ++f) {
4034:     PetscInt        fdof, fcomp, cfdof;
4035:     const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4036:     PetscInt        cind = 0, k, c;

4038:     PetscSectionGetFieldComponents(section, f, &fcomp);
4039:     PetscSectionGetFieldDof(section, point, f, &fdof);
4040:     PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
4041:     if (!cfdof || setBC) {
4042:       if (orientation >= 0) {
4043:         for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
4044:       } else {
4045:         for (k = fdof/fcomp-1; k >= 0; --k) {
4046:           for (c = 0; c < fcomp; ++c) {
4047:             indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
4048:           }
4049:         }
4050:       }
4051:     } else {
4052:       PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
4053:       if (orientation >= 0) {
4054:         for (k = 0; k < fdof; ++k) {
4055:           if ((cind < cfdof) && (k == fcdofs[cind])) {
4056:             indices[foffs[f]+k] = -(off+foff+k+1);
4057:             ++cind;
4058:           } else {
4059:             indices[foffs[f]+k] = off+foff+k-cind;
4060:           }
4061:         }
4062:       } else {
4063:         for (k = fdof/fcomp-1; k >= 0; --k) {
4064:           for (c = 0; c < fcomp; ++c) {
4065:             if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
4066:               indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
4067:               ++cind;
4068:             } else {
4069:               indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
4070:             }
4071:           }
4072:         }
4073:       }
4074:     }
4075:     foff     += (setBC ? fdof : (fdof - cfdof));
4076:     foffs[f] += fdof;
4077:   }
4078:   return(0);
4079: }

4083: 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)
4084: {
4085:   Mat             cMat;
4086:   PetscSection    aSec, cSec;
4087:   IS              aIS;
4088:   PetscInt        aStart = -1, aEnd = -1;
4089:   const PetscInt  *anchors;
4090:   PetscInt        numFields, f, p, q, newP = 0;
4091:   PetscInt        newNumPoints = 0, newNumIndices = 0;
4092:   PetscInt        *newPoints, *indices, *newIndices;
4093:   PetscInt        maxAnchor, maxDof;
4094:   PetscInt        newOffsets[32];
4095:   PetscInt        *pointMatOffsets[32];
4096:   PetscInt        *newPointOffsets[32];
4097:   PetscScalar     *pointMat[32];
4098:   PetscScalar     *newValues=NULL,*tmpValues;
4099:   PetscBool       anyConstrained = PETSC_FALSE;
4100:   PetscErrorCode  ierr;

4105:   PetscSectionGetNumFields(section, &numFields);

4107:   DMPlexGetAnchors(dm,&aSec,&aIS);
4108:   /* if there are point-to-point constraints */
4109:   if (aSec) {
4110:     PetscMemzero(newOffsets, 32 * sizeof(PetscInt));
4111:     ISGetIndices(aIS,&anchors);
4112:     PetscSectionGetChart(aSec,&aStart,&aEnd);
4113:     /* figure out how many points are going to be in the new element matrix
4114:      * (we allow double counting, because it's all just going to be summed
4115:      * into the global matrix anyway) */
4116:     for (p = 0; p < 2*numPoints; p+=2) {
4117:       PetscInt b    = points[p];
4118:       PetscInt bDof = 0, bSecDof;

4120:       PetscSectionGetDof(section,b,&bSecDof);
4121:       if (!bSecDof) {
4122:         continue;
4123:       }
4124:       if (b >= aStart && b < aEnd) {
4125:         PetscSectionGetDof(aSec,b,&bDof);
4126:       }
4127:       if (bDof) {
4128:         /* this point is constrained */
4129:         /* it is going to be replaced by its anchors */
4130:         PetscInt bOff, q;

4132:         anyConstrained = PETSC_TRUE;
4133:         newNumPoints  += bDof;
4134:         PetscSectionGetOffset(aSec,b,&bOff);
4135:         for (q = 0; q < bDof; q++) {
4136:           PetscInt a = anchors[bOff + q];
4137:           PetscInt aDof;

4139:           PetscSectionGetDof(section,a,&aDof);
4140:           newNumIndices += aDof;
4141:           for (f = 0; f < numFields; ++f) {
4142:             PetscInt fDof;

4144:             PetscSectionGetFieldDof(section, a, f, &fDof);
4145:             newOffsets[f+1] += fDof;
4146:           }
4147:         }
4148:       }
4149:       else {
4150:         /* this point is not constrained */
4151:         newNumPoints++;
4152:         newNumIndices += bSecDof;
4153:         for (f = 0; f < numFields; ++f) {
4154:           PetscInt fDof;

4156:           PetscSectionGetFieldDof(section, b, f, &fDof);
4157:           newOffsets[f+1] += fDof;
4158:         }
4159:       }
4160:     }
4161:   }
4162:   if (!anyConstrained) {
4163:     if (outNumPoints)  *outNumPoints  = 0;
4164:     if (outNumIndices) *outNumIndices = 0;
4165:     if (outPoints)     *outPoints     = NULL;
4166:     if (outValues)     *outValues     = NULL;
4167:     if (aSec) {ISRestoreIndices(aIS,&anchors);}
4168:     return(0);
4169:   }

4171:   if (outNumPoints)  *outNumPoints  = newNumPoints;
4172:   if (outNumIndices) *outNumIndices = newNumIndices;

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

4176:   if (!outPoints && !outValues) {
4177:     if (offsets) {
4178:       for (f = 0; f <= numFields; f++) {
4179:         offsets[f] = newOffsets[f];
4180:       }
4181:     }
4182:     if (aSec) {ISRestoreIndices(aIS,&anchors);}
4183:     return(0);
4184:   }

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

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

4190:   /* workspaces */
4191:   if (numFields) {
4192:     for (f = 0; f < numFields; f++) {
4193:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4194:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4195:     }
4196:   }
4197:   else {
4198:     DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4199:     DMGetWorkArray(dm,numPoints,PETSC_INT,&newPointOffsets[0]);
4200:   }

4202:   /* get workspaces for the point-to-point matrices */
4203:   if (numFields) {
4204:     PetscInt totalOffset, totalMatOffset;

4206:     for (p = 0; p < numPoints; p++) {
4207:       PetscInt b    = points[2*p];
4208:       PetscInt bDof = 0, bSecDof;

4210:       PetscSectionGetDof(section,b,&bSecDof);
4211:       if (!bSecDof) {
4212:         for (f = 0; f < numFields; f++) {
4213:           newPointOffsets[f][p + 1] = 0;
4214:           pointMatOffsets[f][p + 1] = 0;
4215:         }
4216:         continue;
4217:       }
4218:       if (b >= aStart && b < aEnd) {
4219:         PetscSectionGetDof(aSec, b, &bDof);
4220:       }
4221:       if (bDof) {
4222:         for (f = 0; f < numFields; f++) {
4223:           PetscInt fDof, q, bOff, allFDof = 0;

4225:           PetscSectionGetFieldDof(section, b, f, &fDof);
4226:           PetscSectionGetOffset(aSec, b, &bOff);
4227:           for (q = 0; q < bDof; q++) {
4228:             PetscInt a = anchors[bOff + q];
4229:             PetscInt aFDof;

4231:             PetscSectionGetFieldDof(section, a, f, &aFDof);
4232:             allFDof += aFDof;
4233:           }
4234:           newPointOffsets[f][p+1] = allFDof;
4235:           pointMatOffsets[f][p+1] = fDof * allFDof;
4236:         }
4237:       }
4238:       else {
4239:         for (f = 0; f < numFields; f++) {
4240:           PetscInt fDof;

4242:           PetscSectionGetFieldDof(section, b, f, &fDof);
4243:           newPointOffsets[f][p+1] = fDof;
4244:           pointMatOffsets[f][p+1] = 0;
4245:         }
4246:       }
4247:     }
4248:     for (f = 0, totalOffset = 0, totalMatOffset = 0; f < numFields; f++) {
4249:       newPointOffsets[f][0] = totalOffset;
4250:       pointMatOffsets[f][0] = totalMatOffset;
4251:       for (p = 0; p < numPoints; p++) {
4252:         newPointOffsets[f][p+1] += newPointOffsets[f][p];
4253:         pointMatOffsets[f][p+1] += pointMatOffsets[f][p];
4254:       }
4255:       totalOffset    = newPointOffsets[f][numPoints];
4256:       totalMatOffset = pointMatOffsets[f][numPoints];
4257:       DMGetWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4258:     }
4259:   }
4260:   else {
4261:     for (p = 0; p < numPoints; p++) {
4262:       PetscInt b    = points[2*p];
4263:       PetscInt bDof = 0, bSecDof;

4265:       PetscSectionGetDof(section,b,&bSecDof);
4266:       if (!bSecDof) {
4267:         newPointOffsets[0][p + 1] = 0;
4268:         pointMatOffsets[0][p + 1] = 0;
4269:         continue;
4270:       }
4271:       if (b >= aStart && b < aEnd) {
4272:         PetscSectionGetDof(aSec, b, &bDof);
4273:       }
4274:       if (bDof) {
4275:         PetscInt bOff, q, allDof = 0;

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

4281:           PetscSectionGetDof(section, a, &aDof);
4282:           allDof += aDof;
4283:         }
4284:         newPointOffsets[0][p+1] = allDof;
4285:         pointMatOffsets[0][p+1] = bSecDof * allDof;
4286:       }
4287:       else {
4288:         newPointOffsets[0][p+1] = bSecDof;
4289:         pointMatOffsets[0][p+1] = 0;
4290:       }
4291:     }
4292:     newPointOffsets[0][0] = 0;
4293:     pointMatOffsets[0][0] = 0;
4294:     for (p = 0; p < numPoints; p++) {
4295:       newPointOffsets[0][p+1] += newPointOffsets[0][p];
4296:       pointMatOffsets[0][p+1] += pointMatOffsets[0][p];
4297:     }
4298:     DMGetWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4299:   }

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

4304:   /* get the point-to-point matrices; construct newPoints */
4305:   PetscSectionGetMaxDof(aSec, &maxAnchor);
4306:   PetscSectionGetMaxDof(section, &maxDof);
4307:   DMGetWorkArray(dm,maxDof,PETSC_INT,&indices);
4308:   DMGetWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4309:   if (numFields) {
4310:     for (p = 0, newP = 0; p < numPoints; p++) {
4311:       PetscInt b    = points[2*p];
4312:       PetscInt o    = points[2*p+1];
4313:       PetscInt bDof = 0, bSecDof;

4315:       PetscSectionGetDof(section, b, &bSecDof);
4316:       if (!bSecDof) {
4317:         continue;
4318:       }
4319:       if (b >= aStart && b < aEnd) {
4320:         PetscSectionGetDof(aSec, b, &bDof);
4321:       }
4322:       if (bDof) {
4323:         PetscInt fStart[32], fEnd[32], fAnchorStart[32], fAnchorEnd[32], bOff, q;

4325:         fStart[0] = 0;
4326:         fEnd[0]   = 0;
4327:         for (f = 0; f < numFields; f++) {
4328:           PetscInt fDof;

4330:           PetscSectionGetFieldDof(cSec, b, f, &fDof);
4331:           fStart[f+1] = fStart[f] + fDof;
4332:           fEnd[f+1]   = fStart[f+1];
4333:         }
4334:         PetscSectionGetOffset(cSec, b, &bOff);
4335:         DMPlexGetIndicesPointFields_Internal(cSec, b, bOff, fEnd, PETSC_TRUE, o, indices);

4337:         fAnchorStart[0] = 0;
4338:         fAnchorEnd[0]   = 0;
4339:         for (f = 0; f < numFields; f++) {
4340:           PetscInt fDof = newPointOffsets[f][p + 1] - newPointOffsets[f][p];

4342:           fAnchorStart[f+1] = fAnchorStart[f] + fDof;
4343:           fAnchorEnd[f+1]   = fAnchorStart[f + 1];
4344:         }
4345:         PetscSectionGetOffset(aSec, b, &bOff);
4346:         for (q = 0; q < bDof; q++) {
4347:           PetscInt a = anchors[bOff + q], aOff;

4349:           /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4350:           newPoints[2*(newP + q)]     = a;
4351:           newPoints[2*(newP + q) + 1] = 0;
4352:           PetscSectionGetOffset(section, a, &aOff);
4353:           DMPlexGetIndicesPointFields_Internal(section, a, aOff, fAnchorEnd, PETSC_TRUE, 0, newIndices);
4354:         }
4355:         newP += bDof;

4357:         if (outValues) {
4358:           /* get the point-to-point submatrix */
4359:           for (f = 0; f < numFields; f++) {
4360:             MatGetValues(cMat,fEnd[f]-fStart[f],indices + fStart[f],fAnchorEnd[f] - fAnchorStart[f],newIndices + fAnchorStart[f],pointMat[f] + pointMatOffsets[f][p]);
4361:           }
4362:         }
4363:       }
4364:       else {
4365:         newPoints[2 * newP]     = b;
4366:         newPoints[2 * newP + 1] = o;
4367:         newP++;
4368:       }
4369:     }
4370:   } else {
4371:     for (p = 0; p < numPoints; p++) {
4372:       PetscInt b    = points[2*p];
4373:       PetscInt o    = points[2*p+1];
4374:       PetscInt bDof = 0, bSecDof;

4376:       PetscSectionGetDof(section, b, &bSecDof);
4377:       if (!bSecDof) {
4378:         continue;
4379:       }
4380:       if (b >= aStart && b < aEnd) {
4381:         PetscSectionGetDof(aSec, b, &bDof);
4382:       }
4383:       if (bDof) {
4384:         PetscInt bEnd = 0, bAnchorEnd = 0, bOff;

4386:         PetscSectionGetOffset(cSec, b, &bOff);
4387:         DMPlexGetIndicesPoint_Internal(cSec, b, bOff, &bEnd, PETSC_TRUE, o, indices);

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

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

4395:           newPoints[2*(newP + q)]     = a;
4396:           newPoints[2*(newP + q) + 1] = 0;
4397:           PetscSectionGetOffset(section, a, &aOff);
4398:           DMPlexGetIndicesPoint_Internal(section, a, aOff, &bAnchorEnd, PETSC_TRUE, 0, newIndices);
4399:         }
4400:         newP += bDof;

4402:         /* get the point-to-point submatrix */
4403:         if (outValues) {
4404:           MatGetValues(cMat,bEnd,indices,bAnchorEnd,newIndices,pointMat[0] + pointMatOffsets[0][p]);
4405:         }
4406:       }
4407:       else {
4408:         newPoints[2 * newP]     = b;
4409:         newPoints[2 * newP + 1] = o;
4410:         newP++;
4411:       }
4412:     }
4413:   }

4415:   if (outValues) {
4416:     DMGetWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4417:     PetscMemzero(tmpValues,newNumIndices*numIndices*sizeof(*tmpValues));
4418:     /* multiply constraints on the right */
4419:     if (numFields) {
4420:       for (f = 0; f < numFields; f++) {
4421:         PetscInt oldOff = offsets[f];

4423:         for (p = 0; p < numPoints; p++) {
4424:           PetscInt cStart = newPointOffsets[f][p];
4425:           PetscInt b      = points[2 * p];
4426:           PetscInt c, r, k;
4427:           PetscInt dof;

4429:           PetscSectionGetFieldDof(section,b,f,&dof);
4430:           if (!dof) {
4431:             continue;
4432:           }
4433:           if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4434:             PetscInt nCols         = newPointOffsets[f][p+1]-cStart;
4435:             const PetscScalar *mat = pointMat[f] + pointMatOffsets[f][p];

4437:             for (r = 0; r < numIndices; r++) {
4438:               for (c = 0; c < nCols; c++) {
4439:                 for (k = 0; k < dof; k++) {
4440:                   tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4441:                 }
4442:               }
4443:             }
4444:           }
4445:           else {
4446:             /* copy this column as is */
4447:             for (r = 0; r < numIndices; r++) {
4448:               for (c = 0; c < dof; c++) {
4449:                 tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4450:               }
4451:             }
4452:           }
4453:           oldOff += dof;
4454:         }
4455:       }
4456:     }
4457:     else {
4458:       PetscInt oldOff = 0;
4459:       for (p = 0; p < numPoints; p++) {
4460:         PetscInt cStart = newPointOffsets[0][p];
4461:         PetscInt b      = points[2 * p];
4462:         PetscInt c, r, k;
4463:         PetscInt dof;

4465:         PetscSectionGetDof(section,b,&dof);
4466:         if (!dof) {
4467:           continue;
4468:         }
4469:         if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4470:           PetscInt nCols         = newPointOffsets[0][p+1]-cStart;
4471:           const PetscScalar *mat = pointMat[0] + pointMatOffsets[0][p];

4473:           for (r = 0; r < numIndices; r++) {
4474:             for (c = 0; c < nCols; c++) {
4475:               for (k = 0; k < dof; k++) {
4476:                 tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4477:               }
4478:             }
4479:           }
4480:         }
4481:         else {
4482:           /* copy this column as is */
4483:           for (r = 0; r < numIndices; r++) {
4484:             for (c = 0; c < dof; c++) {
4485:               tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4486:             }
4487:           }
4488:         }
4489:         oldOff += dof;
4490:       }
4491:     }

4493:     if (multiplyLeft) {
4494:       DMGetWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4495:       PetscMemzero(newValues,newNumIndices*newNumIndices*sizeof(*newValues));
4496:       /* multiply constraints transpose on the left */
4497:       if (numFields) {
4498:         for (f = 0; f < numFields; f++) {
4499:           PetscInt oldOff = offsets[f];

4501:           for (p = 0; p < numPoints; p++) {
4502:             PetscInt rStart = newPointOffsets[f][p];
4503:             PetscInt b      = points[2 * p];
4504:             PetscInt c, r, k;
4505:             PetscInt dof;

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

4512:               for (r = 0; r < nRows; r++) {
4513:                 for (c = 0; c < newNumIndices; c++) {
4514:                   for (k = 0; k < dof; k++) {
4515:                     newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4516:                   }
4517:                 }
4518:               }
4519:             }
4520:             else {
4521:               /* copy this row as is */
4522:               for (r = 0; r < dof; r++) {
4523:                 for (c = 0; c < newNumIndices; c++) {
4524:                   newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4525:                 }
4526:               }
4527:             }
4528:             oldOff += dof;
4529:           }
4530:         }
4531:       }
4532:       else {
4533:         PetscInt oldOff = 0;

4535:         for (p = 0; p < numPoints; p++) {
4536:           PetscInt rStart = newPointOffsets[0][p];
4537:           PetscInt b      = points[2 * p];
4538:           PetscInt c, r, k;
4539:           PetscInt dof;

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

4546:             for (r = 0; r < nRows; r++) {
4547:               for (c = 0; c < newNumIndices; c++) {
4548:                 for (k = 0; k < dof; k++) {
4549:                   newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4550:                 }
4551:               }
4552:             }
4553:           }
4554:           else {
4555:             /* copy this row as is */
4556:             for (r = 0; r < dof; r++) {
4557:               for (c = 0; c < newNumIndices; c++) {
4558:                 newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4559:               }
4560:             }
4561:           }
4562:           oldOff += dof;
4563:         }
4564:       }

4566:       DMRestoreWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4567:     }
4568:     else {
4569:       newValues = tmpValues;
4570:     }
4571:   }

4573:   /* clean up */
4574:   DMRestoreWorkArray(dm,maxDof,PETSC_INT,&indices);
4575:   DMRestoreWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);

4577:   if (numFields) {
4578:     for (f = 0; f < numFields; f++) {
4579:       DMRestoreWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4580:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4581:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4582:     }
4583:   }
4584:   else {
4585:     DMRestoreWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4586:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4587:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[0]);
4588:   }
4589:   ISRestoreIndices(aIS,&anchors);

4591:   /* output */
4592:   if (outPoints) {
4593:     *outPoints = newPoints;
4594:   }
4595:   else {
4596:     DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4597:   }
4598:   if (outValues) {
4599:     *outValues = newValues;
4600:   }
4601:   for (f = 0; f <= numFields; f++) {
4602:     offsets[f] = newOffsets[f];
4603:   }
4604:   return(0);
4605: }

4609: PetscErrorCode DMPlexGetClosureIndices(DM dm, PetscSection section, PetscSection globalSection, PetscInt point, PetscInt *numIndices, PetscInt **indices, PetscInt *outOffsets)
4610: {
4611:   PetscSection    clSection;
4612:   IS              clPoints;
4613:   const PetscInt *clp;
4614:   PetscInt       *points = NULL, *pointsNew;
4615:   PetscInt        numPoints, numPointsNew;
4616:   PetscInt        offsets[32];
4617:   PetscInt        Nf, Nind, NindNew, off, globalOff, f, p;
4618:   PetscErrorCode  ierr;

4626:   PetscSectionGetNumFields(section, &Nf);
4627:   if (Nf > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", Nf);
4628:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
4629:   /* Get points in closure */
4630:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4631:   if (!clPoints) {
4632:     PetscInt pStart, pEnd, q;

4634:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4635:     /* Compress out points not in the section */
4636:     PetscSectionGetChart(section, &pStart, &pEnd);
4637:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4638:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4639:         points[q*2]   = points[p];
4640:         points[q*2+1] = points[p+1];
4641:         ++q;
4642:       }
4643:     }
4644:     numPoints = q;
4645:   } else {
4646:     PetscInt dof, off;

4648:     PetscSectionGetDof(clSection, point, &dof);
4649:     numPoints = dof/2;
4650:     PetscSectionGetOffset(clSection, point, &off);
4651:     ISGetIndices(clPoints, &clp);
4652:     points = (PetscInt *) &clp[off];
4653:   }
4654:   /* Get number of indices and indices per field */
4655:   for (p = 0, Nind = 0; p < numPoints*2; p += 2) {
4656:     PetscInt dof, fdof;

4658:     PetscSectionGetDof(section, points[p], &dof);
4659:     for (f = 0; f < Nf; ++f) {
4660:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
4661:       offsets[f+1] += fdof;
4662:     }
4663:     Nind += dof;
4664:   }
4665:   for (f = 1; f < Nf; ++f) offsets[f+1] += offsets[f];
4666:   if (Nf && offsets[Nf] != Nind) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[Nf], Nind);
4667:   /* Correct for hanging node constraints */
4668:   {
4669:     DMPlexAnchorsModifyMat(dm, section, numPoints, Nind, points, NULL, &numPointsNew, &NindNew, &pointsNew, NULL, offsets, PETSC_TRUE);
4670:     if (numPointsNew) {
4671:       if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4672:       else           {ISRestoreIndices(clPoints, &clp);}
4673:       numPoints = numPointsNew;
4674:       Nind      = NindNew;
4675:       points    = pointsNew;
4676:     }
4677:   }
4678:   /* Calculate indices */
4679:   DMGetWorkArray(dm, Nind, PETSC_INT, indices);
4680:   if (Nf) {
4681:     if (outOffsets) {
4682:       PetscInt f;

4684:       for (f = 0; f <= Nf; f++) {
4685:         outOffsets[f] = offsets[f];
4686:       }
4687:     }
4688:     for (p = 0; p < numPoints*2; p += 2) {
4689:       PetscInt o = points[p+1];
4690:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4691:       DMPlexGetIndicesPointFields_Internal(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, *indices);
4692:     }
4693:   } else {
4694:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
4695:       PetscInt o = points[p+1];
4696:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4697:       DMPlexGetIndicesPoint_Internal(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, *indices);
4698:     }
4699:   }
4700:   /* Cleanup points */
4701:   if (numPointsNew) {
4702:     DMRestoreWorkArray(dm, 2*numPointsNew, PETSC_INT, &pointsNew);
4703:   } else {
4704:     if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4705:     else           {ISRestoreIndices(clPoints, &clp);}
4706:   }
4707:   if (numIndices) *numIndices = Nind;
4708:   return(0);
4709: }

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

4720:   DMRestoreWorkArray(dm, 0, PETSC_INT, indices);
4721:   return(0);
4722: }

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

4729:   Not collective

4731:   Input Parameters:
4732: + dm - The DM
4733: . section - The section describing the layout in v, or NULL to use the default section
4734: . globalSection - The section describing the layout in v, or NULL to use the default global section
4735: . A - The matrix
4736: . point - The sieve point in the DM
4737: . values - The array of values
4738: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

4743:   Level: intermediate

4745: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
4746: @*/
4747: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4748: {
4749:   DM_Plex        *mesh   = (DM_Plex*) dm->data;
4750:   PetscSection    clSection;
4751:   IS              clPoints;
4752:   PetscInt       *points = NULL, *newPoints;
4753:   const PetscInt *clp;
4754:   PetscInt       *indices;
4755:   PetscInt        offsets[32];
4756:   PetscInt        numFields, numPoints, newNumPoints, numIndices, newNumIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
4757:   PetscScalar    *newValues;
4758:   PetscErrorCode  ierr;

4762:   if (!section) {DMGetDefaultSection(dm, &section);}
4764:   if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
4767:   PetscSectionGetNumFields(section, &numFields);
4768:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4769:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
4770:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4771:   if (!clPoints) {
4772:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4773:     /* Compress out points not in the section */
4774:     PetscSectionGetChart(section, &pStart, &pEnd);
4775:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4776:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4777:         points[q*2]   = points[p];
4778:         points[q*2+1] = points[p+1];
4779:         ++q;
4780:       }
4781:     }
4782:     numPoints = q;
4783:   } else {
4784:     PetscInt dof, off;

4786:     PetscSectionGetDof(clSection, point, &dof);
4787:     numPoints = dof/2;
4788:     PetscSectionGetOffset(clSection, point, &off);
4789:     ISGetIndices(clPoints, &clp);
4790:     points = (PetscInt *) &clp[off];
4791:   }
4792:   for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
4793:     PetscInt fdof;

4795:     PetscSectionGetDof(section, points[p], &dof);
4796:     for (f = 0; f < numFields; ++f) {
4797:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
4798:       offsets[f+1] += fdof;
4799:     }
4800:     numIndices += dof;
4801:   }
4802:   for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];

4804:   if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
4805:   DMPlexAnchorsModifyMat(dm,section,numPoints,numIndices,points,values,&newNumPoints,&newNumIndices,&newPoints,&newValues,offsets,PETSC_TRUE);
4806:   if (newNumPoints) {
4807:     if (!clPoints) {
4808:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4809:     } else {
4810:       ISRestoreIndices(clPoints, &clp);
4811:     }
4812:     numPoints  = newNumPoints;
4813:     numIndices = newNumIndices;
4814:     points     = newPoints;
4815:     values     = newValues;
4816:   }
4817:   DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
4818:   if (numFields) {
4819:     for (p = 0; p < numPoints*2; p += 2) {
4820:       PetscInt o = points[p+1];
4821:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4822:       DMPlexGetIndicesPointFields_Internal(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
4823:     }
4824:   } else {
4825:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
4826:       PetscInt o = points[p+1];
4827:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4828:       DMPlexGetIndicesPoint_Internal(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
4829:     }
4830:   }
4831:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
4832:   MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
4833:   if (mesh->printFEM > 1) {
4834:     PetscInt i;
4835:     PetscPrintf(PETSC_COMM_SELF, "  Indices:");
4836:     for (i = 0; i < numIndices; ++i) {PetscPrintf(PETSC_COMM_SELF, " %d", indices[i]);}
4837:     PetscPrintf(PETSC_COMM_SELF, "\n");
4838:   }
4839:   if (ierr) {
4840:     PetscMPIInt    rank;
4841:     PetscErrorCode ierr2;

4843:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4844:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4845:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
4846:     ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
4847: 
4848:   }
4849:   if (newNumPoints) {
4850:     DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4851:     DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4852:   }
4853:   else {
4854:     if (!clPoints) {
4855:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4856:     } else {
4857:       ISRestoreIndices(clPoints, &clp);
4858:     }
4859:   }
4860:   DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
4861:   return(0);
4862: }

4866: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4867: {
4868:   DM_Plex        *mesh   = (DM_Plex*) dmf->data;
4869:   PetscInt       *fpoints = NULL, *ftotpoints = NULL;
4870:   PetscInt       *cpoints = NULL;
4871:   PetscInt       *findices, *cindices;
4872:   PetscInt        foffsets[32], coffsets[32];
4873:   CellRefiner     cellRefiner;
4874:   PetscInt        numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4875:   PetscErrorCode  ierr;

4880:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4882:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
4884:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4886:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4889:   PetscSectionGetNumFields(fsection, &numFields);
4890:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4891:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4892:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4893:   /* Column indices */
4894:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4895:   maxFPoints = numCPoints;
4896:   /* Compress out points not in the section */
4897:   /*   TODO: Squeeze out points with 0 dof as well */
4898:   PetscSectionGetChart(csection, &pStart, &pEnd);
4899:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4900:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4901:       cpoints[q*2]   = cpoints[p];
4902:       cpoints[q*2+1] = cpoints[p+1];
4903:       ++q;
4904:     }
4905:   }
4906:   numCPoints = q;
4907:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4908:     PetscInt fdof;

4910:     PetscSectionGetDof(csection, cpoints[p], &dof);
4911:     if (!dof) continue;
4912:     for (f = 0; f < numFields; ++f) {
4913:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4914:       coffsets[f+1] += fdof;
4915:     }
4916:     numCIndices += dof;
4917:   }
4918:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4919:   /* Row indices */
4920:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4921:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4922:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4923:   for (r = 0, q = 0; r < numSubcells; ++r) {
4924:     /* TODO Map from coarse to fine cells */
4925:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4926:     /* Compress out points not in the section */
4927:     PetscSectionGetChart(fsection, &pStart, &pEnd);
4928:     for (p = 0; p < numFPoints*2; p += 2) {
4929:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4930:         PetscSectionGetDof(fsection, fpoints[p], &dof);
4931:         if (!dof) continue;
4932:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4933:         if (s < q) continue;
4934:         ftotpoints[q*2]   = fpoints[p];
4935:         ftotpoints[q*2+1] = fpoints[p+1];
4936:         ++q;
4937:       }
4938:     }
4939:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4940:   }
4941:   numFPoints = q;
4942:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4943:     PetscInt fdof;

4945:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4946:     if (!dof) continue;
4947:     for (f = 0; f < numFields; ++f) {
4948:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4949:       foffsets[f+1] += fdof;
4950:     }
4951:     numFIndices += dof;
4952:   }
4953:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

4955:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4956:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4957:   DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4958:   DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4959:   if (numFields) {
4960:     for (p = 0; p < numFPoints*2; p += 2) {
4961:       PetscInt o = ftotpoints[p+1];
4962:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4963:       DMPlexGetIndicesPointFields_Internal(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4964:     }
4965:     for (p = 0; p < numCPoints*2; p += 2) {
4966:       PetscInt o = cpoints[p+1];
4967:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4968:       DMPlexGetIndicesPointFields_Internal(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4969:     }
4970:   } else {
4971:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4972:       PetscInt o = ftotpoints[p+1];
4973:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4974:       DMPlexGetIndicesPoint_Internal(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4975:     }
4976:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4977:       PetscInt o = cpoints[p+1];
4978:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4979:       DMPlexGetIndicesPoint_Internal(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4980:     }
4981:   }
4982:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
4983:   MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
4984:   if (ierr) {
4985:     PetscMPIInt    rank;
4986:     PetscErrorCode ierr2;

4988:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4989:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4990:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
4991:     ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
4992:     ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
4993: 
4994:   }
4995:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4996:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4997:   DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4998:   DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4999:   return(0);
5000: }

5004: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
5005: {
5006:   PetscInt      *fpoints = NULL, *ftotpoints = NULL;
5007:   PetscInt      *cpoints = NULL;
5008:   PetscInt       foffsets[32], coffsets[32];
5009:   CellRefiner    cellRefiner;
5010:   PetscInt       numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;

5016:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
5018:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
5020:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
5022:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
5024:   PetscSectionGetNumFields(fsection, &numFields);
5025:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
5026:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
5027:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
5028:   /* Column indices */
5029:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5030:   maxFPoints = numCPoints;
5031:   /* Compress out points not in the section */
5032:   /*   TODO: Squeeze out points with 0 dof as well */
5033:   PetscSectionGetChart(csection, &pStart, &pEnd);
5034:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
5035:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
5036:       cpoints[q*2]   = cpoints[p];
5037:       cpoints[q*2+1] = cpoints[p+1];
5038:       ++q;
5039:     }
5040:   }
5041:   numCPoints = q;
5042:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
5043:     PetscInt fdof;

5045:     PetscSectionGetDof(csection, cpoints[p], &dof);
5046:     if (!dof) continue;
5047:     for (f = 0; f < numFields; ++f) {
5048:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
5049:       coffsets[f+1] += fdof;
5050:     }
5051:     numCIndices += dof;
5052:   }
5053:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
5054:   /* Row indices */
5055:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
5056:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
5057:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
5058:   for (r = 0, q = 0; r < numSubcells; ++r) {
5059:     /* TODO Map from coarse to fine cells */
5060:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
5061:     /* Compress out points not in the section */
5062:     PetscSectionGetChart(fsection, &pStart, &pEnd);
5063:     for (p = 0; p < numFPoints*2; p += 2) {
5064:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
5065:         PetscSectionGetDof(fsection, fpoints[p], &dof);
5066:         if (!dof) continue;
5067:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
5068:         if (s < q) continue;
5069:         ftotpoints[q*2]   = fpoints[p];
5070:         ftotpoints[q*2+1] = fpoints[p+1];
5071:         ++q;
5072:       }
5073:     }
5074:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
5075:   }
5076:   numFPoints = q;
5077:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
5078:     PetscInt fdof;

5080:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
5081:     if (!dof) continue;
5082:     for (f = 0; f < numFields; ++f) {
5083:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
5084:       foffsets[f+1] += fdof;
5085:     }
5086:     numFIndices += dof;
5087:   }
5088:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

5090:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
5091:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
5092:   if (numFields) {
5093:     for (p = 0; p < numFPoints*2; p += 2) {
5094:       PetscInt o = ftotpoints[p+1];
5095:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5096:       DMPlexGetIndicesPointFields_Internal(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
5097:     }
5098:     for (p = 0; p < numCPoints*2; p += 2) {
5099:       PetscInt o = cpoints[p+1];
5100:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5101:       DMPlexGetIndicesPointFields_Internal(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5102:     }
5103:   } else {
5104:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5105:       PetscInt o = ftotpoints[p+1];
5106:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5107:       DMPlexGetIndicesPoint_Internal(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5108:     }
5109:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5110:       PetscInt o = cpoints[p+1];
5111:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5112:       DMPlexGetIndicesPoint_Internal(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5113:     }
5114:   }
5115:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5116:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5117:   return(0);
5118: }

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

5125:   Input Parameter:
5126: . dm - The DMPlex object

5128:   Output Parameters:
5129: + cMax - The first hybrid cell
5130: . fMax - The first hybrid face
5131: . eMax - The first hybrid edge
5132: - vMax - The first hybrid vertex

5134:   Level: developer

5136: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
5137: @*/
5138: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
5139: {
5140:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5141:   PetscInt       dim;

5146:   DMGetDimension(dm, &dim);
5147:   if (cMax) *cMax = mesh->hybridPointMax[dim];
5148:   if (fMax) *fMax = mesh->hybridPointMax[dim-1];
5149:   if (eMax) *eMax = mesh->hybridPointMax[1];
5150:   if (vMax) *vMax = mesh->hybridPointMax[0];
5151:   return(0);
5152: }

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

5159:   Input Parameters:
5160: . dm   - The DMPlex object
5161: . cMax - The first hybrid cell
5162: . fMax - The first hybrid face
5163: . eMax - The first hybrid edge
5164: - vMax - The first hybrid vertex

5166:   Level: developer

5168: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5169: @*/
5170: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5171: {
5172:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5173:   PetscInt       dim;

5178:   DMGetDimension(dm, &dim);
5179:   if (cMax >= 0) mesh->hybridPointMax[dim]   = cMax;
5180:   if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5181:   if (eMax >= 0) mesh->hybridPointMax[1]     = eMax;
5182:   if (vMax >= 0) mesh->hybridPointMax[0]     = vMax;
5183:   return(0);
5184: }

5188: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5189: {
5190:   DM_Plex *mesh = (DM_Plex*) dm->data;

5195:   *cellHeight = mesh->vtkCellHeight;
5196:   return(0);
5197: }

5201: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5202: {
5203:   DM_Plex *mesh = (DM_Plex*) dm->data;

5207:   mesh->vtkCellHeight = cellHeight;
5208:   return(0);
5209: }

5213: /* We can easily have a form that takes an IS instead */
5214: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5215: {
5216:   PetscSection   section, globalSection;
5217:   PetscInt      *numbers, p;

5221:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
5222:   PetscSectionSetChart(section, pStart, pEnd);
5223:   for (p = pStart; p < pEnd; ++p) {
5224:     PetscSectionSetDof(section, p, 1);
5225:   }
5226:   PetscSectionSetUp(section);
5227:   PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
5228:   PetscMalloc1(pEnd - pStart, &numbers);
5229:   for (p = pStart; p < pEnd; ++p) {
5230:     PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5231:     if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5232:     else                       numbers[p-pStart] += shift;
5233:   }
5234:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5235:   if (globalSize) {
5236:     PetscLayout layout;
5237:     PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5238:     PetscLayoutGetSize(layout, globalSize);
5239:     PetscLayoutDestroy(&layout);
5240:   }
5241:   PetscSectionDestroy(&section);
5242:   PetscSectionDestroy(&globalSection);
5243:   return(0);
5244: }

5248: PetscErrorCode DMPlexCreateCellNumbering_Internal(DM dm, PetscBool includeHybrid, IS *globalCellNumbers)
5249: {
5250:   PetscInt       cellHeight, cStart, cEnd, cMax;

5254:   DMPlexGetVTKCellHeight(dm, &cellHeight);
5255:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5256:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5257:   if (cMax >= 0 && !includeHybrid) cEnd = PetscMin(cEnd, cMax);
5258:   DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, globalCellNumbers);
5259:   return(0);
5260: }

5264: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5265: {
5266:   DM_Plex       *mesh = (DM_Plex*) dm->data;

5271:   if (!mesh->globalCellNumbers) {DMPlexCreateCellNumbering_Internal(dm, PETSC_FALSE, &mesh->globalCellNumbers);}
5272:   *globalCellNumbers = mesh->globalCellNumbers;
5273:   return(0);
5274: }

5278: PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM dm, PetscBool includeHybrid, IS *globalVertexNumbers)
5279: {
5280:   PetscInt       vStart, vEnd, vMax;

5285:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5286:   DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5287:   if (vMax >= 0 && !includeHybrid) vEnd = PetscMin(vEnd, vMax);
5288:   DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, globalVertexNumbers);
5289:   return(0);
5290: }

5294: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5295: {
5296:   DM_Plex       *mesh = (DM_Plex*) dm->data;

5301:   if (!mesh->globalVertexNumbers) {DMPlexCreateVertexNumbering_Internal(dm, PETSC_FALSE, &mesh->globalVertexNumbers);}
5302:   *globalVertexNumbers = mesh->globalVertexNumbers;
5303:   return(0);
5304: }

5308: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5309: {
5310:   IS             nums[4];
5311:   PetscInt       depths[4];
5312:   PetscInt       depth, d, shift = 0;

5317:   DMPlexGetDepth(dm, &depth);
5318:   /* For unstratified meshes use dim instead of depth */
5319:   if (depth < 0) {DMGetDimension(dm, &depth);}
5320:   depths[0] = depth; depths[1] = 0;
5321:   for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5322:   for (d = 0; d <= depth; ++d) {
5323:     PetscInt pStart, pEnd, gsize;

5325:     DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5326:     DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5327:     shift += gsize;
5328:   }
5329:   ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5330:   for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5331:   return(0);
5332: }

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

5339:   Input Parameters:
5340:   + dm - The DMPlex object

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

5344:   Level: developer

5346: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5347: @*/
5348: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5349: {
5350:   PetscSection    coneSection, supportSection;
5351:   const PetscInt *cone, *support;
5352:   PetscInt        coneSize, c, supportSize, s;
5353:   PetscInt        pStart, pEnd, p, csize, ssize;
5354:   PetscErrorCode  ierr;

5358:   DMPlexGetConeSection(dm, &coneSection);
5359:   DMPlexGetSupportSection(dm, &supportSection);
5360:   /* Check that point p is found in the support of its cone points, and vice versa */
5361:   DMPlexGetChart(dm, &pStart, &pEnd);
5362:   for (p = pStart; p < pEnd; ++p) {
5363:     DMPlexGetConeSize(dm, p, &coneSize);
5364:     DMPlexGetCone(dm, p, &cone);
5365:     for (c = 0; c < coneSize; ++c) {
5366:       PetscBool dup = PETSC_FALSE;
5367:       PetscInt  d;
5368:       for (d = c-1; d >= 0; --d) {
5369:         if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5370:       }
5371:       DMPlexGetSupportSize(dm, cone[c], &supportSize);
5372:       DMPlexGetSupport(dm, cone[c], &support);
5373:       for (s = 0; s < supportSize; ++s) {
5374:         if (support[s] == p) break;
5375:       }
5376:       if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5377:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5378:         for (s = 0; s < coneSize; ++s) {
5379:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
5380:         }
5381:         PetscPrintf(PETSC_COMM_SELF, "\n");
5382:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
5383:         for (s = 0; s < supportSize; ++s) {
5384:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
5385:         }
5386:         PetscPrintf(PETSC_COMM_SELF, "\n");
5387:         if (dup) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
5388:         else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
5389:       }
5390:     }
5391:     DMPlexGetSupportSize(dm, p, &supportSize);
5392:     DMPlexGetSupport(dm, p, &support);
5393:     for (s = 0; s < supportSize; ++s) {
5394:       DMPlexGetConeSize(dm, support[s], &coneSize);
5395:       DMPlexGetCone(dm, support[s], &cone);
5396:       for (c = 0; c < coneSize; ++c) {
5397:         if (cone[c] == p) break;
5398:       }
5399:       if (c >= coneSize) {
5400:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
5401:         for (c = 0; c < supportSize; ++c) {
5402:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
5403:         }
5404:         PetscPrintf(PETSC_COMM_SELF, "\n");
5405:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
5406:         for (c = 0; c < coneSize; ++c) {
5407:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
5408:         }
5409:         PetscPrintf(PETSC_COMM_SELF, "\n");
5410:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
5411:       }
5412:     }
5413:   }
5414:   PetscSectionGetStorageSize(coneSection, &csize);
5415:   PetscSectionGetStorageSize(supportSection, &ssize);
5416:   if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
5417:   return(0);
5418: }

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

5425:   Input Parameters:
5426: + dm - The DMPlex object
5427: . isSimplex - Are the cells simplices or tensor products
5428: - cellHeight - Normally 0

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

5432:   Level: developer

5434: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
5435: @*/
5436: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5437: {
5438:   PetscInt       dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;

5443:   DMGetDimension(dm, &dim);
5444:   switch (dim) {
5445:   case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
5446:   case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
5447:   case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
5448:   default:
5449:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
5450:   }
5451:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5452:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5453:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5454:   cMax = cMax >= 0 ? cMax : cEnd;
5455:   for (c = cStart; c < cMax; ++c) {
5456:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5458:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5459:     for (cl = 0; cl < closureSize*2; cl += 2) {
5460:       const PetscInt p = closure[cl];
5461:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5462:     }
5463:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5464:     if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has  %d vertices != %d", c, coneSize, numCorners);
5465:   }
5466:   for (c = cMax; c < cEnd; ++c) {
5467:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5469:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5470:     for (cl = 0; cl < closureSize*2; cl += 2) {
5471:       const PetscInt p = closure[cl];
5472:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5473:     }
5474:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5475:     if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has  %d vertices > %d", c, coneSize, numHybridCorners);
5476:   }
5477:   return(0);
5478: }

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

5485:   Input Parameters:
5486: + dm - The DMPlex object
5487: . isSimplex - Are the cells simplices or tensor products
5488: - cellHeight - Normally 0

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

5492:   Level: developer

5494: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
5495: @*/
5496: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5497: {
5498:   PetscInt       pMax[4];
5499:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, h;

5504:   DMGetDimension(dm, &dim);
5505:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5506:   DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
5507:   for (h = cellHeight; h < dim; ++h) {
5508:     DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
5509:     for (c = cStart; c < cEnd; ++c) {
5510:       const PetscInt *cone, *ornt, *faces;
5511:       PetscInt        numFaces, faceSize, coneSize,f;
5512:       PetscInt       *closure = NULL, closureSize, cl, numCorners = 0;

5514:       if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
5515:       DMPlexGetConeSize(dm, c, &coneSize);
5516:       DMPlexGetCone(dm, c, &cone);
5517:       DMPlexGetConeOrientation(dm, c, &ornt);
5518:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5519:       for (cl = 0; cl < closureSize*2; cl += 2) {
5520:         const PetscInt p = closure[cl];
5521:         if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
5522:       }
5523:       DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
5524:       if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
5525:       for (f = 0; f < numFaces; ++f) {
5526:         PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;

5528:         DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
5529:         for (cl = 0; cl < fclosureSize*2; cl += 2) {
5530:           const PetscInt p = fclosure[cl];
5531:           if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
5532:         }
5533:         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);
5534:         for (v = 0; v < fnumCorners; ++v) {
5535:           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]);
5536:         }
5537:         DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
5538:       }
5539:       DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
5540:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5541:     }
5542:   }
5543:   return(0);
5544: }

5548: /* Pointwise interpolation
5549:      Just code FEM for now
5550:      u^f = I u^c
5551:      sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
5552:      u^f_i = sum_j psi^f_i I phi^c_j u^c_j
5553:      I_{ij} = psi^f_i phi^c_j
5554: */
5555: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
5556: {
5557:   PetscSection   gsc, gsf;
5558:   PetscInt       m, n;
5559:   void          *ctx;
5560:   DM             cdm;
5561:   PetscBool      regular;

5565:   DMGetDefaultGlobalSection(dmFine, &gsf);
5566:   PetscSectionGetConstrainedStorageSize(gsf, &m);
5567:   DMGetDefaultGlobalSection(dmCoarse, &gsc);
5568:   PetscSectionGetConstrainedStorageSize(gsc, &n);

5570:   MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
5571:   MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
5572:   MatSetType(*interpolation, dmCoarse->mattype);
5573:   DMGetApplicationContext(dmFine, &ctx);

5575:   DMGetCoarseDM(dmFine, &cdm);
5576:   DMPlexGetRegularRefinement(dmFine, &regular);
5577:   if (regular && cdm == dmCoarse) {DMPlexComputeInterpolatorNested(dmCoarse, dmFine, *interpolation, ctx);}
5578:   else                            {DMPlexComputeInterpolatorGeneral(dmCoarse, dmFine, *interpolation, ctx);}
5579:   MatViewFromOptions(*interpolation, NULL, "-interp_mat_view");
5580:   /* Use naive scaling */
5581:   DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
5582:   return(0);
5583: }

5587: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat)
5588: {
5590:   VecScatter     ctx;

5593:   DMPlexComputeInjectorFEM(dmCoarse, dmFine, &ctx, NULL);
5594:   MatCreateScatter(PetscObjectComm((PetscObject)ctx), ctx, mat);
5595:   VecScatterDestroy(&ctx);
5596:   return(0);
5597: }

5601: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
5602: {
5603:   PetscSection   section;
5604:   IS            *bcPoints, *bcComps;
5605:   PetscBool     *isFE;
5606:   PetscInt      *bcFields, *numComp, *numDof;
5607:   PetscInt       depth, dim, numBd, numBC = 0, numFields, bd, bc = 0, f;
5608:   PetscInt       cStart, cEnd, cEndInterior;

5612:   DMGetNumFields(dm, &numFields);
5613:   if (!numFields) return(0);
5614:   /* FE and FV boundary conditions are handled slightly differently */
5615:   PetscMalloc1(numFields, &isFE);
5616:   for (f = 0; f < numFields; ++f) {
5617:     PetscObject  obj;
5618:     PetscClassId id;

5620:     DMGetField(dm, f, &obj);
5621:     PetscObjectGetClassId(obj, &id);
5622:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
5623:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
5624:     else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
5625:   }
5626:   /* Allocate boundary point storage for FEM boundaries */
5627:   DMPlexGetDepth(dm, &depth);
5628:   DMGetDimension(dm, &dim);
5629:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
5630:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
5631:   PetscDSGetNumBoundary(dm->prob, &numBd);
5632:   for (bd = 0; bd < numBd; ++bd) {
5633:     PetscInt    field;
5634:     PetscBool   isEssential;
5635:     const char  *labelName;
5636:     DMLabel     label;

5638:     PetscDSGetBoundary(dm->prob, bd, &isEssential, NULL, &labelName, &field, NULL, NULL, NULL, NULL, NULL, NULL);
5639:     DMGetLabel(dm,labelName,&label);
5640:     if (label && isFE[field] && isEssential) ++numBC;
5641:   }
5642:   /* Add ghost cell boundaries for FVM */
5643:   for (f = 0; f < numFields; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
5644:   PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);
5645:   /* Constrain ghost cells for FV */
5646:   for (f = 0; f < numFields; ++f) {
5647:     PetscInt *newidx, c;

5649:     if (isFE[f] || cEndInterior < 0) continue;
5650:     PetscMalloc1(cEnd-cEndInterior,&newidx);
5651:     for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
5652:     bcFields[bc] = f;
5653:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5654:   }
5655:   /* Handle FEM Dirichlet boundaries */
5656:   for (bd = 0; bd < numBd; ++bd) {
5657:     const char     *bdLabel;
5658:     DMLabel         label;
5659:     const PetscInt *comps;
5660:     const PetscInt *values;
5661:     PetscInt        bd2, field, numComps, numValues;
5662:     PetscBool       isEssential, duplicate = PETSC_FALSE;

5664:     PetscDSGetBoundary(dm->prob, bd, &isEssential, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
5665:     DMGetLabel(dm, bdLabel, &label);
5666:     if (!isFE[field] || !label) continue;
5667:     /* Only want to modify label once */
5668:     for (bd2 = 0; bd2 < bd; ++bd2) {
5669:       const char *bdname;
5670:       PetscDSGetBoundary(dm->prob, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
5671:       PetscStrcmp(bdname, bdLabel, &duplicate);
5672:       if (duplicate) break;
5673:     }
5674:     if (!duplicate && (isFE[field])) {
5675:       /* don't complete cells, which are just present to give orientation to the boundary */
5676:       DMPlexLabelComplete(dm, label);
5677:     }
5678:     /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
5679:     if (isEssential) {
5680:       PetscInt       *newidx;
5681:       PetscInt        n, newn = 0, p, v;

5683:       bcFields[bc] = field;
5684:       if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
5685:       for (v = 0; v < numValues; ++v) {
5686:         IS              tmp;
5687:         const PetscInt *idx;

5689:         DMGetStratumIS(dm, bdLabel, values[v], &tmp);
5690:         if (!tmp) continue;
5691:         ISGetLocalSize(tmp, &n);
5692:         ISGetIndices(tmp, &idx);
5693:         if (isFE[field]) {
5694:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
5695:         } else {
5696:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
5697:         }
5698:         ISRestoreIndices(tmp, &idx);
5699:         ISDestroy(&tmp);
5700:       }
5701:       PetscMalloc1(newn,&newidx);
5702:       newn = 0;
5703:       for (v = 0; v < numValues; ++v) {
5704:         IS              tmp;
5705:         const PetscInt *idx;

5707:         DMGetStratumIS(dm, bdLabel, values[v], &tmp);
5708:         if (!tmp) continue;
5709:         ISGetLocalSize(tmp, &n);
5710:         ISGetIndices(tmp, &idx);
5711:         if (isFE[field]) {
5712:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
5713:         } else {
5714:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
5715:         }
5716:         ISRestoreIndices(tmp, &idx);
5717:         ISDestroy(&tmp);
5718:       }
5719:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5720:     }
5721:   }
5722:   /* Handle discretization */
5723:   PetscCalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
5724:   for (f = 0; f < numFields; ++f) {
5725:     PetscObject obj;

5727:     DMGetField(dm, f, &obj);
5728:     if (isFE[f]) {
5729:       PetscFE         fe = (PetscFE) obj;
5730:       const PetscInt *numFieldDof;
5731:       PetscInt        d;

5733:       PetscFEGetNumComponents(fe, &numComp[f]);
5734:       PetscFEGetNumDof(fe, &numFieldDof);
5735:       for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
5736:     } else {
5737:       PetscFV fv = (PetscFV) obj;

5739:       PetscFVGetNumComponents(fv, &numComp[f]);
5740:       numDof[f*(dim+1)+dim] = numComp[f];
5741:     }
5742:   }
5743:   for (f = 0; f < numFields; ++f) {
5744:     PetscInt d;
5745:     for (d = 1; d < dim; ++d) {
5746:       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.");
5747:     }
5748:   }
5749:   DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section);
5750:   for (f = 0; f < numFields; ++f) {
5751:     PetscFE     fe;
5752:     const char *name;

5754:     DMGetField(dm, f, (PetscObject *) &fe);
5755:     PetscObjectGetName((PetscObject) fe, &name);
5756:     PetscSectionSetFieldName(section, f, name);
5757:   }
5758:   DMSetDefaultSection(dm, section);
5759:   PetscSectionDestroy(&section);
5760:   for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
5761:   PetscFree3(bcFields,bcPoints,bcComps);
5762:   PetscFree2(numComp,numDof);
5763:   PetscFree(isFE);
5764:   return(0);
5765: }

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

5772:   Input Parameter:
5773: . dm - The DMPlex object

5775:   Output Parameter:
5776: . regular - The flag

5778:   Level: intermediate

5780: .seealso: DMPlexSetRegularRefinement()
5781: @*/
5782: PetscErrorCode DMPlexGetRegularRefinement(DM dm, PetscBool *regular)
5783: {
5787:   *regular = ((DM_Plex *) dm->data)->regularRefinement;
5788:   return(0);
5789: }

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

5796:   Input Parameters:
5797: + dm - The DMPlex object
5798: - regular - The flag

5800:   Level: intermediate

5802: .seealso: DMPlexGetRegularRefinement()
5803: @*/
5804: PetscErrorCode DMPlexSetRegularRefinement(DM dm, PetscBool regular)
5805: {
5808:   ((DM_Plex *) dm->data)->regularRefinement = regular;
5809:   return(0);
5810: }

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

5819:   not collective

5821:   Input Parameters:
5822: . dm - The DMPlex object

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


5829:   Level: intermediate

5831: .seealso: DMPlexSetAnchors(), DMGetConstraints(), DMSetConstraints()
5832: @*/
5833: PetscErrorCode DMPlexGetAnchors(DM dm, PetscSection *anchorSection, IS *anchorIS)
5834: {
5835:   DM_Plex *plex = (DM_Plex *)dm->data;

5840:   if (!plex->anchorSection && !plex->anchorIS && plex->createanchors) {(*plex->createanchors)(dm);}
5841:   if (anchorSection) *anchorSection = plex->anchorSection;
5842:   if (anchorIS) *anchorIS = plex->anchorIS;
5843:   return(0);
5844: }

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

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

5856:   collective on dm

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

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

5865:   Level: intermediate

5867: .seealso: DMPlexGetAnchors(), DMGetConstraints(), DMSetConstraints()
5868: @*/
5869: PetscErrorCode DMPlexSetAnchors(DM dm, PetscSection anchorSection, IS anchorIS)
5870: {
5871:   DM_Plex        *plex = (DM_Plex *)dm->data;
5872:   PetscMPIInt    result;

5877:   if (anchorSection) {
5879:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorSection),&result);
5880:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor section must have local communicator");
5881:   }
5882:   if (anchorIS) {
5884:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorIS),&result);
5885:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor IS must have local communicator");
5886:   }

5888:   PetscObjectReference((PetscObject)anchorSection);
5889:   PetscSectionDestroy(&plex->anchorSection);
5890:   plex->anchorSection = anchorSection;

5892:   PetscObjectReference((PetscObject)anchorIS);
5893:   ISDestroy(&plex->anchorIS);
5894:   plex->anchorIS = anchorIS;

5896: #if defined(PETSC_USE_DEBUG)
5897:   if (anchorIS && anchorSection) {
5898:     PetscInt size, a, pStart, pEnd;
5899:     const PetscInt *anchors;

5901:     PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5902:     ISGetLocalSize(anchorIS,&size);
5903:     ISGetIndices(anchorIS,&anchors);
5904:     for (a = 0; a < size; a++) {
5905:       PetscInt p;

5907:       p = anchors[a];
5908:       if (p >= pStart && p < pEnd) {
5909:         PetscInt dof;

5911:         PetscSectionGetDof(anchorSection,p,&dof);
5912:         if (dof) {
5913:           PetscErrorCode ierr2;

5915:           ierr2 = ISRestoreIndices(anchorIS,&anchors);CHKERRQ(ierr2);
5916:           SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point %d cannot be constrained and an anchor",p);
5917:         }
5918:       }
5919:     }
5920:     ISRestoreIndices(anchorIS,&anchors);
5921:   }
5922: #endif
5923:   /* reset the generic constraints */
5924:   DMSetDefaultConstraints(dm,NULL,NULL);
5925:   return(0);
5926: }

5930: static PetscErrorCode DMPlexCreateConstraintSection_Anchors(DM dm, PetscSection section, PetscSection *cSec)
5931: {
5932:   PetscSection anchorSection;
5933:   PetscInt pStart, pEnd, sStart, sEnd, p, dof, numFields, f;

5938:   DMPlexGetAnchors(dm,&anchorSection,NULL);
5939:   PetscSectionCreate(PETSC_COMM_SELF,cSec);
5940:   PetscSectionGetNumFields(section,&numFields);
5941:   if (numFields) {
5942:     PetscInt f;
5943:     PetscSectionSetNumFields(*cSec,numFields);

5945:     for (f = 0; f < numFields; f++) {
5946:       PetscInt numComp;

5948:       PetscSectionGetFieldComponents(section,f,&numComp);
5949:       PetscSectionSetFieldComponents(*cSec,f,numComp);
5950:     }
5951:   }
5952:   PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5953:   PetscSectionGetChart(section,&sStart,&sEnd);
5954:   pStart = PetscMax(pStart,sStart);
5955:   pEnd   = PetscMin(pEnd,sEnd);
5956:   pEnd   = PetscMax(pStart,pEnd);
5957:   PetscSectionSetChart(*cSec,pStart,pEnd);
5958:   for (p = pStart; p < pEnd; p++) {
5959:     PetscSectionGetDof(anchorSection,p,&dof);
5960:     if (dof) {
5961:       PetscSectionGetDof(section,p,&dof);
5962:       PetscSectionSetDof(*cSec,p,dof);
5963:       for (f = 0; f < numFields; f++) {
5964:         PetscSectionGetFieldDof(section,p,f,&dof);
5965:         PetscSectionSetFieldDof(*cSec,p,f,dof);
5966:       }
5967:     }
5968:   }
5969:   PetscSectionSetUp(*cSec);
5970:   return(0);
5971: }

5975: static PetscErrorCode DMPlexCreateConstraintMatrix_Anchors(DM dm, PetscSection section, PetscSection cSec, Mat *cMat)
5976: {
5977:   PetscSection aSec;
5978:   PetscInt pStart, pEnd, p, dof, aDof, aOff, off, nnz, annz, m, n, q, a, offset, *i, *j;
5979:   const PetscInt *anchors;
5980:   PetscInt numFields, f;
5981:   IS aIS;

5986:   PetscSectionGetStorageSize(cSec, &m);
5987:   PetscSectionGetStorageSize(section, &n);
5988:   MatCreate(PETSC_COMM_SELF,cMat);
5989:   MatSetSizes(*cMat,m,n,m,n);
5990:   MatSetType(*cMat,MATSEQAIJ);
5991:   DMPlexGetAnchors(dm,&aSec,&aIS);
5992:   ISGetIndices(aIS,&anchors);
5993:   /* cSec will be a subset of aSec and section */
5994:   PetscSectionGetChart(cSec,&pStart,&pEnd);
5995:   PetscMalloc1(m+1,&i);
5996:   i[0] = 0;
5997:   PetscSectionGetNumFields(section,&numFields);
5998:   for (p = pStart; p < pEnd; p++) {
5999:     PetscInt rDof, rOff, r;

6001:     PetscSectionGetDof(aSec,p,&rDof);
6002:     if (!rDof) continue;
6003:     PetscSectionGetOffset(aSec,p,&rOff);
6004:     if (numFields) {
6005:       for (f = 0; f < numFields; f++) {
6006:         annz = 0;
6007:         for (r = 0; r < rDof; r++) {
6008:           a = anchors[rOff + r];
6009:           PetscSectionGetFieldDof(section,a,f,&aDof);
6010:           annz += aDof;
6011:         }
6012:         PetscSectionGetFieldDof(cSec,p,f,&dof);
6013:         PetscSectionGetFieldOffset(cSec,p,f,&off);
6014:         for (q = 0; q < dof; q++) {
6015:           i[off + q + 1] = i[off + q] + annz;
6016:         }
6017:       }
6018:     }
6019:     else {
6020:       annz = 0;
6021:       for (q = 0; q < dof; q++) {
6022:         a = anchors[off + q];
6023:         PetscSectionGetDof(section,a,&aDof);
6024:         annz += aDof;
6025:       }
6026:       PetscSectionGetDof(cSec,p,&dof);
6027:       PetscSectionGetOffset(cSec,p,&off);
6028:       for (q = 0; q < dof; q++) {
6029:         i[off + q + 1] = i[off + q] + annz;
6030:       }
6031:     }
6032:   }
6033:   nnz = i[m];
6034:   PetscMalloc1(nnz,&j);
6035:   offset = 0;
6036:   for (p = pStart; p < pEnd; p++) {
6037:     if (numFields) {
6038:       for (f = 0; f < numFields; f++) {
6039:         PetscSectionGetFieldDof(cSec,p,f,&dof);
6040:         for (q = 0; q < dof; q++) {
6041:           PetscInt rDof, rOff, r;
6042:           PetscSectionGetDof(aSec,p,&rDof);
6043:           PetscSectionGetOffset(aSec,p,&rOff);
6044:           for (r = 0; r < rDof; r++) {
6045:             PetscInt s;

6047:             a = anchors[rOff + r];
6048:             PetscSectionGetFieldDof(section,a,f,&aDof);
6049:             PetscSectionGetFieldOffset(section,a,f,&aOff);
6050:             for (s = 0; s < aDof; s++) {
6051:               j[offset++] = aOff + s;
6052:             }
6053:           }
6054:         }
6055:       }
6056:     }
6057:     else {
6058:       PetscSectionGetDof(cSec,p,&dof);
6059:       for (q = 0; q < dof; q++) {
6060:         PetscInt rDof, rOff, r;
6061:         PetscSectionGetDof(aSec,p,&rDof);
6062:         PetscSectionGetOffset(aSec,p,&rOff);
6063:         for (r = 0; r < rDof; r++) {
6064:           PetscInt s;

6066:           a = anchors[rOff + r];
6067:           PetscSectionGetDof(section,a,&aDof);
6068:           PetscSectionGetOffset(section,a,&aOff);
6069:           for (s = 0; s < aDof; s++) {
6070:             j[offset++] = aOff + s;
6071:           }
6072:         }
6073:       }
6074:     }
6075:   }
6076:   MatSeqAIJSetPreallocationCSR(*cMat,i,j,NULL);
6077:   PetscFree(i);
6078:   PetscFree(j);
6079:   ISRestoreIndices(aIS,&anchors);
6080:   return(0);
6081: }

6085: PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm)
6086: {
6087:   DM_Plex        *plex = (DM_Plex *)dm->data;
6088:   PetscSection   anchorSection, section, cSec;
6089:   Mat            cMat;

6094:   DMPlexGetAnchors(dm,&anchorSection,NULL);
6095:   if (anchorSection) {
6096:     PetscDS  ds;
6097:     PetscInt nf;

6099:     DMGetDefaultSection(dm,&section);
6100:     DMPlexCreateConstraintSection_Anchors(dm,section,&cSec);
6101:     DMPlexCreateConstraintMatrix_Anchors(dm,section,cSec,&cMat);
6102:     DMGetDS(dm,&ds);
6103:     PetscDSGetNumFields(ds,&nf);
6104:     if (nf && plex->computeanchormatrix) {(*plex->computeanchormatrix)(dm,section,cSec,cMat);}
6105:     DMSetDefaultConstraints(dm,cSec,cMat);
6106:     PetscSectionDestroy(&cSec);
6107:     MatDestroy(&cMat);
6108:   }
6109:   return(0);
6110: }