Actual source code: plex.c

petsc-master 2015-08-04
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <petsc/private/isimpl.h>
  3: #include <petscsf.h>
  4: #include <petscds.h>

  6: /* Logging support */
  7: 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;

  9: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
 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, 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:     PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
371:     PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
372:     PetscViewerASCIISynchronizedPrintf(viewer, "Max sizes cone: %D support: %D\n", maxConeSize, maxSupportSize);
373:     PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
374:     PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
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:     PetscSectionGetChart(coordSection, &pStart, NULL);
397:     if (pStart >= 0) {PetscSectionVecView(coordSection, coordinates, viewer);}
398:     DMPlexGetLabel(dm, "marker", &markers);
399:     DMLabelView(markers,viewer);
400:     if (size > 1) {
401:       PetscSF sf;

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

419:     DMGetDimension(dm, &dim);
420:     DMPlexGetDepth(dm, &depth);
421:     DMPlexGetNumLabels(dm, &numLabels);
422:     numLabels  = PetscMax(numLabels, 10);
423:     numColors  = 10;
424:     numLColors = 10;
425:     PetscCalloc3(numLabels, &names, numColors, &colors, numLColors, &lcolors);
426:     PetscOptionsGetReal(((PetscObject) viewer)->prefix, "-dm_plex_view_scale", &scale, NULL);
427:     PetscOptionsGetBool(((PetscObject) viewer)->prefix, "-dm_plex_view_numbers", &useNumbers, NULL);
428:     PetscOptionsGetStringArray(((PetscObject) viewer)->prefix, "-dm_plex_view_labels", names, &numLabels, &useLabels);
429:     if (!useLabels) numLabels = 0;
430:     PetscOptionsGetStringArray(((PetscObject) viewer)->prefix, "-dm_plex_view_colors", colors, &numColors, &useColors);
431:     if (!useColors) {
432:       numColors = 3;
433:       for (c = 0; c < numColors; ++c) {PetscStrallocpy(defcolors[c], &colors[c]);}
434:     }
435:     PetscOptionsGetStringArray(((PetscObject) viewer)->prefix, "-dm_plex_view_lcolors", lcolors, &numLColors, &useColors);
436:     if (!useColors) {
437:       numLColors = 4;
438:       for (c = 0; c < numLColors; ++c) {PetscStrallocpy(deflcolors[c], &lcolors[c]);}
439:     }
440:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
441:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
442:     PetscObjectGetName((PetscObject) dm, &name);
443:     PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
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:     for (v = vStart; v < vEnd; ++v) {
467:       PetscInt  off, dof, d;
468:       PetscBool isLabeled = PETSC_FALSE;

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

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

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

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

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

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

635:       DMPlexGetLabelName(dm, l, &name);
636:       DMPlexGetLabel(dm, name, &label);
637:       DMLabelGetNumValues(label, &numValues);
638:       PetscViewerASCIIPrintf(viewer, "  %s: %d strata of sizes (", name, numValues);
639:       DMLabelGetValueIS(label, &valueIS);
640:       ISGetIndices(valueIS, &values);
641:       for (v = 0; v < numValues; ++v) {
642:         PetscInt size;

644:         DMLabelGetStratumSize(label, values[v], &size);
645:         if (v > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
646:         PetscViewerASCIIPrintf(viewer, "%d", size);
647:       }
648:       PetscViewerASCIIPrintf(viewer, ")\n");
649:       ISRestoreIndices(valueIS, &values);
650:       ISDestroy(&valueIS);
651:     }
652:   }
653:   return(0);
654: }

658: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
659: {
660:   PetscBool      iascii, ishdf5, isvtk;

666:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
667:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,   &isvtk);
668:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,  &ishdf5);
669:   if (iascii) {
670:     DMPlexView_Ascii(dm, viewer);
671:   } else if (ishdf5) {
672: #if defined(PETSC_HAVE_HDF5)
673:     PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
674:     DMPlexView_HDF5(dm, viewer);
675:     PetscViewerPopFormat(viewer);
676: #else
677:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
678: #endif
679:   }
680:   else if (isvtk) {
681:     DMPlexVTKWriteAll((PetscObject) dm,viewer);
682:   }
683:   return(0);
684: }

688: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
689: {
690:   PetscBool      isbinary, ishdf5;

696:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
697:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,   &ishdf5);
698:   if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");}
699:   else if (ishdf5) {
700: #if defined(PETSC_HAVE_HDF5)
701:     DMPlexLoad_HDF5(dm, viewer);
702: #else
703:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
704: #endif
705:   }
706:   return(0);
707: }

711: static PetscErrorCode BoundaryDestroy(DMBoundary *boundary)
712: {
713:   DMBoundary     b, next;

717:   if (!boundary) return(0);
718:   b = *boundary;
719:   *boundary = NULL;
720:   for (; b; b = next) {
721:     next = b->next;
722:     PetscFree(b->comps);
723:     PetscFree(b->ids);
724:     PetscFree(b->name);
725:     PetscFree(b->labelname);
726:     PetscFree(b);
727:   }
728:   return(0);
729: }

733: PetscErrorCode DMDestroy_Plex(DM dm)
734: {
735:   DM_Plex       *mesh = (DM_Plex*) dm->data;
736:   PlexLabel      next = mesh->labels;

740:   if (--mesh->refct > 0) return(0);
741:   PetscSectionDestroy(&mesh->coneSection);
742:   PetscFree(mesh->cones);
743:   PetscFree(mesh->coneOrientations);
744:   PetscSectionDestroy(&mesh->supportSection);
745:   PetscFree(mesh->supports);
746:   PetscFree(mesh->facesTmp);
747:   PetscFree(mesh->tetgenOpts);
748:   PetscFree(mesh->triangleOpts);
749:   PetscPartitionerDestroy(&mesh->partitioner);
750:   while (next) {
751:     PlexLabel tmp = next->next;

753:     DMLabelDestroy(&next->label);
754:     PetscFree(next);
755:     next = tmp;
756:   }
757:   DMDestroy(&mesh->coarseMesh);
758:   DMLabelDestroy(&mesh->subpointMap);
759:   ISDestroy(&mesh->globalVertexNumbers);
760:   ISDestroy(&mesh->globalCellNumbers);
761:   BoundaryDestroy(&mesh->boundary);
762:   PetscSectionDestroy(&mesh->anchorSection);
763:   ISDestroy(&mesh->anchorIS);
764:   PetscSectionDestroy(&mesh->parentSection);
765:   PetscFree(mesh->parents);
766:   PetscFree(mesh->childIDs);
767:   PetscSectionDestroy(&mesh->childSection);
768:   PetscFree(mesh->children);
769:   DMDestroy(&mesh->referenceTree);
770:   /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
771:   PetscFree(mesh);
772:   return(0);
773: }

777: PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
778: {
779:   PetscSection   sectionGlobal;
780:   PetscInt       bs = -1;
781:   PetscInt       localSize;
782:   PetscBool      isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock;
784:   MatType        mtype;
785:   ISLocalToGlobalMapping ltog;

788:   MatInitializePackage();
789:   mtype = dm->mattype;
790:   DMGetDefaultGlobalSection(dm, &sectionGlobal);
791:   /* PetscSectionGetStorageSize(sectionGlobal, &localSize); */
792:   PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
793:   MatCreate(PetscObjectComm((PetscObject)dm), J);
794:   MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
795:   MatSetType(*J, mtype);
796:   MatSetFromOptions(*J);
797:   PetscStrcmp(mtype, MATSHELL, &isShell);
798:   PetscStrcmp(mtype, MATBAIJ, &isBlock);
799:   PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
800:   PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
801:   PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
802:   PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
803:   PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
804:   if (!isShell) {
805:     PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
806:     PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin;

808:     if (bs < 0) {
809:       if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
810:         PetscInt pStart, pEnd, p, dof, cdof;

812:         PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
813:         for (p = pStart; p < pEnd; ++p) {
814:           PetscSectionGetDof(sectionGlobal, p, &dof);
815:           PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
816:           if (dof-cdof) {
817:             if (bs < 0) {
818:               bs = dof-cdof;
819:             } else if (bs != dof-cdof) {
820:               /* Layout does not admit a pointwise block size */
821:               bs = 1;
822:               break;
823:             }
824:           }
825:         }
826:         /* Must have same blocksize on all procs (some might have no points) */
827:         bsLocal = bs;
828:         MPI_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
829:         bsLocal = bs < 0 ? bsMax : bs;
830:         MPI_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
831:         if (bsMin != bsMax) {
832:           bs = 1;
833:         } else {
834:           bs = bsMax;
835:         }
836:       } else {
837:         bs = 1;
838:       }
839:     }
840:     PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
841:     DMPlexPreallocateOperator(dm, bs, dnz, onz, dnzu, onzu, *J, fillMatrix);
842:     PetscFree4(dnz, onz, dnzu, onzu);

844:     /* Set localtoglobalmapping on the matrix for MatSetValuesLocal() to work */
845:     DMGetLocalToGlobalMapping(dm,&ltog);
846:     MatSetLocalToGlobalMapping(*J,ltog,ltog);
847:   }
848:   return(0);
849: }

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

856:   Not collective

858:   Input Parameter:
859: . mesh - The DMPlex

861:   Output Parameters:
862: + pStart - The first mesh point
863: - pEnd   - The upper bound for mesh points

865:   Level: beginner

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

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

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

885:   Not collective

887:   Input Parameters:
888: + mesh - The DMPlex
889: . pStart - The first mesh point
890: - pEnd   - The upper bound for mesh points

892:   Output Parameters:

894:   Level: beginner

896: .seealso: DMPlexCreate(), DMPlexGetChart()
897: @*/
898: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
899: {
900:   DM_Plex       *mesh = (DM_Plex*) dm->data;

905:   PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
906:   PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
907:   return(0);
908: }

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

915:   Not collective

917:   Input Parameters:
918: + mesh - The DMPlex
919: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

921:   Output Parameter:
922: . size - The cone size for point p

924:   Level: beginner

926: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
927: @*/
928: PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size)
929: {
930:   DM_Plex       *mesh = (DM_Plex*) dm->data;

936:   PetscSectionGetDof(mesh->coneSection, p, size);
937:   return(0);
938: }

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

945:   Not collective

947:   Input Parameters:
948: + mesh - The DMPlex
949: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
950: - size - The cone size for point p

952:   Output Parameter:

954:   Note:
955:   This should be called after DMPlexSetChart().

957:   Level: beginner

959: .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart()
960: @*/
961: PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size)
962: {
963:   DM_Plex       *mesh = (DM_Plex*) dm->data;

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

970:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
971:   return(0);
972: }

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

979:   Not collective

981:   Input Parameters:
982: + mesh - The DMPlex
983: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
984: - size - The additional cone size for point p

986:   Output Parameter:

988:   Note:
989:   This should be called after DMPlexSetChart().

991:   Level: beginner

993: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart()
994: @*/
995: PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size)
996: {
997:   DM_Plex       *mesh = (DM_Plex*) dm->data;
998:   PetscInt       csize;

1003:   PetscSectionAddDof(mesh->coneSection, p, size);
1004:   PetscSectionGetDof(mesh->coneSection, p, &csize);

1006:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
1007:   return(0);
1008: }

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

1015:   Not collective

1017:   Input Parameters:
1018: + mesh - The DMPlex
1019: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1024:   Level: beginner

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

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

1032: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart()
1033: @*/
1034: PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
1035: {
1036:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1037:   PetscInt       off;

1043:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1044:   *cone = &mesh->cones[off];
1045:   return(0);
1046: }

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

1053:   Not collective

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

1060:   Output Parameter:

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

1065:   Level: beginner

1067: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1068: @*/
1069: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
1070: {
1071:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1072:   PetscInt       pStart, pEnd;
1073:   PetscInt       dof, off, c;

1078:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1079:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1081:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1082:   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);
1083:   for (c = 0; c < dof; ++c) {
1084:     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);
1085:     mesh->cones[off+c] = cone[c];
1086:   }
1087:   return(0);
1088: }

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

1095:   Not collective

1097:   Input Parameters:
1098: + mesh - The DMPlex
1099: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1107:   Level: beginner

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

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

1115: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
1116: @*/
1117: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
1118: {
1119:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1120:   PetscInt       off;

1125: #if defined(PETSC_USE_DEBUG)
1126:   {
1127:     PetscInt dof;
1128:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1130:   }
1131: #endif
1132:   PetscSectionGetOffset(mesh->coneSection, p, &off);

1134:   *coneOrientation = &mesh->coneOrientations[off];
1135:   return(0);
1136: }

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

1143:   Not collective

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

1153:   Output Parameter:

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

1158:   Level: beginner

1160: .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1161: @*/
1162: PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
1163: {
1164:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1165:   PetscInt       pStart, pEnd;
1166:   PetscInt       dof, off, c;

1171:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1172:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1174:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1175:   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);
1176:   for (c = 0; c < dof; ++c) {
1177:     PetscInt cdof, o = coneOrientation[c];

1179:     PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1180:     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);
1181:     mesh->coneOrientations[off+c] = o;
1182:   }
1183:   return(0);
1184: }

1188: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
1189: {
1190:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1191:   PetscInt       pStart, pEnd;
1192:   PetscInt       dof, off;

1197:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1198:   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);
1199:   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);
1200:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1201:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1202:   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);
1203:   mesh->cones[off+conePos] = conePoint;
1204:   return(0);
1205: }

1209: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1210: {
1211:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1212:   PetscInt       pStart, pEnd;
1213:   PetscInt       dof, off;

1218:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1219:   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);
1220:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1221:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1222:   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);
1223:   mesh->coneOrientations[off+conePos] = coneOrientation;
1224:   return(0);
1225: }

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

1232:   Not collective

1234:   Input Parameters:
1235: + mesh - The DMPlex
1236: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1238:   Output Parameter:
1239: . size - The support size for point p

1241:   Level: beginner

1243: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1244: @*/
1245: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1246: {
1247:   DM_Plex       *mesh = (DM_Plex*) dm->data;

1253:   PetscSectionGetDof(mesh->supportSection, p, size);
1254:   return(0);
1255: }

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

1262:   Not collective

1264:   Input Parameters:
1265: + mesh - The DMPlex
1266: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1267: - size - The support size for point p

1269:   Output Parameter:

1271:   Note:
1272:   This should be called after DMPlexSetChart().

1274:   Level: beginner

1276: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1277: @*/
1278: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1279: {
1280:   DM_Plex       *mesh = (DM_Plex*) dm->data;

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

1287:   mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1288:   return(0);
1289: }

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

1296:   Not collective

1298:   Input Parameters:
1299: + mesh - The DMPlex
1300: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1305:   Level: beginner

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

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

1313: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1314: @*/
1315: PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1316: {
1317:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1318:   PetscInt       off;

1324:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1325:   *support = &mesh->supports[off];
1326:   return(0);
1327: }

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

1334:   Not collective

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

1341:   Output Parameter:

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

1346:   Level: beginner

1348: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1349: @*/
1350: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1351: {
1352:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1353:   PetscInt       pStart, pEnd;
1354:   PetscInt       dof, off, c;

1359:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1360:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1362:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1363:   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);
1364:   for (c = 0; c < dof; ++c) {
1365:     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);
1366:     mesh->supports[off+c] = support[c];
1367:   }
1368:   return(0);
1369: }

1373: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1374: {
1375:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1376:   PetscInt       pStart, pEnd;
1377:   PetscInt       dof, off;

1382:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1383:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1384:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1385:   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);
1386:   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);
1387:   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);
1388:   mesh->supports[off+supportPos] = supportPoint;
1389:   return(0);
1390: }

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

1397:   Not collective

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

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

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

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

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

1418:   Level: beginner

1420: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1421: @*/
1422: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1423: {
1424:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1425:   PetscInt       *closure, *fifo;
1426:   const PetscInt *tmp = NULL, *tmpO = NULL;
1427:   PetscInt        tmpSize, t;
1428:   PetscInt        depth       = 0, maxSize;
1429:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1430:   PetscErrorCode  ierr;

1434:   DMPlexGetDepth(dm, &depth);
1435:   /* This is only 1-level */
1436:   if (useCone) {
1437:     DMPlexGetConeSize(dm, p, &tmpSize);
1438:     DMPlexGetCone(dm, p, &tmp);
1439:     DMPlexGetConeOrientation(dm, p, &tmpO);
1440:   } else {
1441:     DMPlexGetSupportSize(dm, p, &tmpSize);
1442:     DMPlexGetSupport(dm, p, &tmp);
1443:   }
1444:   if (depth == 1) {
1445:     if (*points) {
1446:       closure = *points;
1447:     } else {
1448:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1449:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1450:     }
1451:     closure[0] = p; closure[1] = 0;
1452:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1453:       closure[closureSize]   = tmp[t];
1454:       closure[closureSize+1] = tmpO ? tmpO[t] : 0;
1455:     }
1456:     if (numPoints) *numPoints = closureSize/2;
1457:     if (points)    *points    = closure;
1458:     return(0);
1459:   }
1460:   {
1461:     PetscInt c, coneSeries, s,supportSeries;

1463:     c = mesh->maxConeSize;
1464:     coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1465:     s = mesh->maxSupportSize;
1466:     supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1467:     maxSize = 2*PetscMax(coneSeries,supportSeries);
1468:   }
1469:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1470:   if (*points) {
1471:     closure = *points;
1472:   } else {
1473:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1474:   }
1475:   closure[0] = p; closure[1] = 0;
1476:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1477:     const PetscInt cp = tmp[t];
1478:     const PetscInt co = tmpO ? tmpO[t] : 0;

1480:     closure[closureSize]   = cp;
1481:     closure[closureSize+1] = co;
1482:     fifo[fifoSize]         = cp;
1483:     fifo[fifoSize+1]       = co;
1484:   }
1485:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1486:   while (fifoSize - fifoStart) {
1487:     const PetscInt q   = fifo[fifoStart];
1488:     const PetscInt o   = fifo[fifoStart+1];
1489:     const PetscInt rev = o >= 0 ? 0 : 1;
1490:     const PetscInt off = rev ? -(o+1) : o;

1492:     if (useCone) {
1493:       DMPlexGetConeSize(dm, q, &tmpSize);
1494:       DMPlexGetCone(dm, q, &tmp);
1495:       DMPlexGetConeOrientation(dm, q, &tmpO);
1496:     } else {
1497:       DMPlexGetSupportSize(dm, q, &tmpSize);
1498:       DMPlexGetSupport(dm, q, &tmp);
1499:       tmpO = NULL;
1500:     }
1501:     for (t = 0; t < tmpSize; ++t) {
1502:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1503:       const PetscInt cp = tmp[i];
1504:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1505:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1506:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1507:       PetscInt       co = tmpO ? tmpO[i] : 0;
1508:       PetscInt       c;

1510:       if (rev) {
1511:         PetscInt childSize, coff;
1512:         DMPlexGetConeSize(dm, cp, &childSize);
1513:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1514:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1515:       }
1516:       /* Check for duplicate */
1517:       for (c = 0; c < closureSize; c += 2) {
1518:         if (closure[c] == cp) break;
1519:       }
1520:       if (c == closureSize) {
1521:         closure[closureSize]   = cp;
1522:         closure[closureSize+1] = co;
1523:         fifo[fifoSize]         = cp;
1524:         fifo[fifoSize+1]       = co;
1525:         closureSize           += 2;
1526:         fifoSize              += 2;
1527:       }
1528:     }
1529:     fifoStart += 2;
1530:   }
1531:   if (numPoints) *numPoints = closureSize/2;
1532:   if (points)    *points    = closure;
1533:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1534:   return(0);
1535: }

1539: /*@C
1540:   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

1542:   Not collective

1544:   Input Parameters:
1545: + mesh - The DMPlex
1546: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1547: . orientation - The orientation of the point
1548: . useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
1549: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used

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

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

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

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

1564:   Level: beginner

1566: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1567: @*/
1568: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1569: {
1570:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1571:   PetscInt       *closure, *fifo;
1572:   const PetscInt *tmp = NULL, *tmpO = NULL;
1573:   PetscInt        tmpSize, t;
1574:   PetscInt        depth       = 0, maxSize;
1575:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1576:   PetscErrorCode  ierr;

1580:   DMPlexGetDepth(dm, &depth);
1581:   /* This is only 1-level */
1582:   if (useCone) {
1583:     DMPlexGetConeSize(dm, p, &tmpSize);
1584:     DMPlexGetCone(dm, p, &tmp);
1585:     DMPlexGetConeOrientation(dm, p, &tmpO);
1586:   } else {
1587:     DMPlexGetSupportSize(dm, p, &tmpSize);
1588:     DMPlexGetSupport(dm, p, &tmp);
1589:   }
1590:   if (depth == 1) {
1591:     if (*points) {
1592:       closure = *points;
1593:     } else {
1594:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1595:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1596:     }
1597:     closure[0] = p; closure[1] = ornt;
1598:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1599:       const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1600:       closure[closureSize]   = tmp[i];
1601:       closure[closureSize+1] = tmpO ? tmpO[i] : 0;
1602:     }
1603:     if (numPoints) *numPoints = closureSize/2;
1604:     if (points)    *points    = closure;
1605:     return(0);
1606:   }
1607:   {
1608:     PetscInt c, coneSeries, s,supportSeries;

1610:     c = mesh->maxConeSize;
1611:     coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1612:     s = mesh->maxSupportSize;
1613:     supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1614:     maxSize = 2*PetscMax(coneSeries,supportSeries);
1615:   }
1616:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1617:   if (*points) {
1618:     closure = *points;
1619:   } else {
1620:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1621:   }
1622:   closure[0] = p; closure[1] = ornt;
1623:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1624:     const PetscInt i  = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1625:     const PetscInt cp = tmp[i];
1626:     PetscInt       co = tmpO ? tmpO[i] : 0;

1628:     if (ornt < 0) {
1629:       PetscInt childSize, coff;
1630:       DMPlexGetConeSize(dm, cp, &childSize);
1631:       coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1632:       co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1633:     }
1634:     closure[closureSize]   = cp;
1635:     closure[closureSize+1] = co;
1636:     fifo[fifoSize]         = cp;
1637:     fifo[fifoSize+1]       = co;
1638:   }
1639:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1640:   while (fifoSize - fifoStart) {
1641:     const PetscInt q   = fifo[fifoStart];
1642:     const PetscInt o   = fifo[fifoStart+1];
1643:     const PetscInt rev = o >= 0 ? 0 : 1;
1644:     const PetscInt off = rev ? -(o+1) : o;

1646:     if (useCone) {
1647:       DMPlexGetConeSize(dm, q, &tmpSize);
1648:       DMPlexGetCone(dm, q, &tmp);
1649:       DMPlexGetConeOrientation(dm, q, &tmpO);
1650:     } else {
1651:       DMPlexGetSupportSize(dm, q, &tmpSize);
1652:       DMPlexGetSupport(dm, q, &tmp);
1653:       tmpO = NULL;
1654:     }
1655:     for (t = 0; t < tmpSize; ++t) {
1656:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1657:       const PetscInt cp = tmp[i];
1658:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1659:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1660:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1661:       PetscInt       co = tmpO ? tmpO[i] : 0;
1662:       PetscInt       c;

1664:       if (rev) {
1665:         PetscInt childSize, coff;
1666:         DMPlexGetConeSize(dm, cp, &childSize);
1667:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1668:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1669:       }
1670:       /* Check for duplicate */
1671:       for (c = 0; c < closureSize; c += 2) {
1672:         if (closure[c] == cp) break;
1673:       }
1674:       if (c == closureSize) {
1675:         closure[closureSize]   = cp;
1676:         closure[closureSize+1] = co;
1677:         fifo[fifoSize]         = cp;
1678:         fifo[fifoSize+1]       = co;
1679:         closureSize           += 2;
1680:         fifoSize              += 2;
1681:       }
1682:     }
1683:     fifoStart += 2;
1684:   }
1685:   if (numPoints) *numPoints = closureSize/2;
1686:   if (points)    *points    = closure;
1687:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1688:   return(0);
1689: }

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

1696:   Not collective

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

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

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

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

1714:   Level: beginner

1716: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1717: @*/
1718: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1719: {

1726:   DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1727:   if (numPoints) *numPoints = 0;
1728:   return(0);
1729: }

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

1736:   Not collective

1738:   Input Parameter:
1739: . mesh - The DMPlex

1741:   Output Parameters:
1742: + maxConeSize - The maximum number of in-edges
1743: - maxSupportSize - The maximum number of out-edges

1745:   Level: beginner

1747: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1748: @*/
1749: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1750: {
1751:   DM_Plex *mesh = (DM_Plex*) dm->data;

1755:   if (maxConeSize)    *maxConeSize    = mesh->maxConeSize;
1756:   if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1757:   return(0);
1758: }

1762: PetscErrorCode DMSetUp_Plex(DM dm)
1763: {
1764:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1765:   PetscInt       size;

1770:   PetscSectionSetUp(mesh->coneSection);
1771:   PetscSectionGetStorageSize(mesh->coneSection, &size);
1772:   PetscMalloc1(size, &mesh->cones);
1773:   PetscCalloc1(size, &mesh->coneOrientations);
1774:   if (mesh->maxSupportSize) {
1775:     PetscSectionSetUp(mesh->supportSection);
1776:     PetscSectionGetStorageSize(mesh->supportSection, &size);
1777:     PetscMalloc1(size, &mesh->supports);
1778:   }
1779:   return(0);
1780: }

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

1789:   if (subdm) {DMClone(dm, subdm);}
1790:   DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1791:   return(0);
1792: }

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

1799:   Not collective

1801:   Input Parameter:
1802: . mesh - The DMPlex

1804:   Output Parameter:

1806:   Note:
1807:   This should be called after all calls to DMPlexSetCone()

1809:   Level: beginner

1811: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1812: @*/
1813: PetscErrorCode DMPlexSymmetrize(DM dm)
1814: {
1815:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1816:   PetscInt      *offsets;
1817:   PetscInt       supportSize;
1818:   PetscInt       pStart, pEnd, p;

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

1829:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1830:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1831:     for (c = off; c < off+dof; ++c) {
1832:       PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1833:     }
1834:   }
1835:   for (p = pStart; p < pEnd; ++p) {
1836:     PetscInt dof;

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

1840:     mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1841:   }
1842:   PetscSectionSetUp(mesh->supportSection);
1843:   /* Calculate supports */
1844:   PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1845:   PetscMalloc1(supportSize, &mesh->supports);
1846:   PetscCalloc1(pEnd - pStart, &offsets);
1847:   for (p = pStart; p < pEnd; ++p) {
1848:     PetscInt dof, off, c;

1850:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1851:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1852:     for (c = off; c < off+dof; ++c) {
1853:       const PetscInt q = mesh->cones[c];
1854:       PetscInt       offS;

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

1858:       mesh->supports[offS+offsets[q]] = p;
1859:       ++offsets[q];
1860:     }
1861:   }
1862:   PetscFree(offsets);
1863:   return(0);
1864: }

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

1874:   Not collective

1876:   Input Parameter:
1877: . mesh - The DMPlex

1879:   Output Parameter:

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

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

1889:   Level: beginner

1891: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1892: @*/
1893: PetscErrorCode DMPlexStratify(DM dm)
1894: {
1895:   DMLabel        label;
1896:   PetscInt       pStart, pEnd, p;
1897:   PetscInt       numRoots = 0, numLeaves = 0;

1902:   PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1903:   /* Calculate depth */
1904:   DMPlexGetChart(dm, &pStart, &pEnd);
1905:   DMPlexCreateLabel(dm, "depth");
1906:   DMPlexGetDepthLabel(dm, &label);
1907:   /* Initialize roots and count leaves */
1908:   for (p = pStart; p < pEnd; ++p) {
1909:     PetscInt coneSize, supportSize;

1911:     DMPlexGetConeSize(dm, p, &coneSize);
1912:     DMPlexGetSupportSize(dm, p, &supportSize);
1913:     if (!coneSize && supportSize) {
1914:       ++numRoots;
1915:       DMLabelSetValue(label, p, 0);
1916:     } else if (!supportSize && coneSize) {
1917:       ++numLeaves;
1918:     } else if (!supportSize && !coneSize) {
1919:       /* Isolated points */
1920:       DMLabelSetValue(label, p, 0);
1921:     }
1922:   }
1923:   if (numRoots + numLeaves == (pEnd - pStart)) {
1924:     for (p = pStart; p < pEnd; ++p) {
1925:       PetscInt coneSize, supportSize;

1927:       DMPlexGetConeSize(dm, p, &coneSize);
1928:       DMPlexGetSupportSize(dm, p, &supportSize);
1929:       if (!supportSize && coneSize) {
1930:         DMLabelSetValue(label, p, 1);
1931:       }
1932:     }
1933:   } else {
1934:     IS       pointIS;
1935:     PetscInt numPoints = 0, level = 0;

1937:     DMLabelGetStratumIS(label, level, &pointIS);
1938:     if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1939:     while (numPoints) {
1940:       const PetscInt *points;
1941:       const PetscInt  newLevel = level+1;

1943:       ISGetIndices(pointIS, &points);
1944:       for (p = 0; p < numPoints; ++p) {
1945:         const PetscInt  point = points[p];
1946:         const PetscInt *support;
1947:         PetscInt        supportSize, s;

1949:         DMPlexGetSupportSize(dm, point, &supportSize);
1950:         DMPlexGetSupport(dm, point, &support);
1951:         for (s = 0; s < supportSize; ++s) {
1952:           DMLabelSetValue(label, support[s], newLevel);
1953:         }
1954:       }
1955:       ISRestoreIndices(pointIS, &points);
1956:       ++level;
1957:       ISDestroy(&pointIS);
1958:       DMLabelGetStratumIS(label, level, &pointIS);
1959:       if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1960:       else         {numPoints = 0;}
1961:     }
1962:     ISDestroy(&pointIS);
1963:   }
1964:   PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1965:   return(0);
1966: }

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

1973:   Not Collective

1975:   Input Parameters:
1976: + dm - The DMPlex object
1977: . numPoints - The number of input points for the join
1978: - points - The input points

1980:   Output Parameters:
1981: + numCoveredPoints - The number of points in the join
1982: - coveredPoints - The points in the join

1984:   Level: intermediate

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

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

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

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

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

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

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

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

2048:   Not Collective

2050:   Input Parameters:
2051: + dm - The DMPlex object
2052: . numPoints - The number of input points for the join
2053: - points - The input points

2055:   Output Parameters:
2056: + numCoveredPoints - The number of points in the join
2057: - coveredPoints - The points in the join

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

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

2065:   Level: intermediate

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

2079:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2080:   if (numCoveredPoints) *numCoveredPoints = 0;
2081:   return(0);
2082: }

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

2089:   Not Collective

2091:   Input Parameters:
2092: + dm - The DMPlex object
2093: . numPoints - The number of input points for the join
2094: - points - The input points

2096:   Output Parameters:
2097: + numCoveredPoints - The number of points in the join
2098: - coveredPoints - The points in the join

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

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

2106:   Level: intermediate

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


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

2134:   for (p = 0; p < numPoints; ++p) {
2135:     PetscInt closureSize;

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

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

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

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

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

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

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

2198:   Not Collective

2200:   Input Parameters:
2201: + dm - The DMPlex object
2202: . numPoints - The number of input points for the meet
2203: - points - The input points

2205:   Output Parameters:
2206: + numCoveredPoints - The number of points in the meet
2207: - coveredPoints - The points in the meet

2209:   Level: intermediate

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

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

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

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

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

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

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

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

2273:   Not Collective

2275:   Input Parameters:
2276: + dm - The DMPlex object
2277: . numPoints - The number of input points for the meet
2278: - points - The input points

2280:   Output Parameters:
2281: + numCoveredPoints - The number of points in the meet
2282: - coveredPoints - The points in the meet

2284:   Level: intermediate

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

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

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

2304:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2305:   if (numCoveredPoints) *numCoveredPoints = 0;
2306:   return(0);
2307: }

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

2314:   Not Collective

2316:   Input Parameters:
2317: + dm - The DMPlex object
2318: . numPoints - The number of input points for the meet
2319: - points - The input points

2321:   Output Parameters:
2322: + numCoveredPoints - The number of points in the meet
2323: - coveredPoints - The points in the meet

2325:   Level: intermediate

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

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

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


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

2359:   for (p = 0; p < numPoints; ++p) {
2360:     PetscInt closureSize;

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

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

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

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

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

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

2420: /*@C
2421:   DMPlexEqual - Determine if two DMs have the same topology

2423:   Not Collective

2425:   Input Parameters:
2426: + dmA - A DMPlex object
2427: - dmB - A DMPlex object

2429:   Output Parameters:
2430: . equal - PETSC_TRUE if the topologies are identical

2432:   Level: intermediate

2434:   Notes:
2435:   We are not solving graph isomorphism, so we do not permutation.

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


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

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

2487: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2488: {
2489:   MPI_Comm       comm;

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

2558: /*@
2559:   DMPlexLocalizeCoordinate - If a mesh is periodic (a torus with lengths L_i, some of which can be infinite), project the coordinate onto [0, L_i) in each dimension.

2561:   Input Parameters:
2562: + dm     - The DM
2563: - in     - The input coordinate point (dim numbers)

2565:   Output Parameter:
2566: . out - The localized coordinate point

2568:   Level: developer

2570: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeAddCoordinate()
2571: @*/
2572: PetscErrorCode DMPlexLocalizeCoordinate(DM dm, const PetscScalar in[], PetscScalar out[])
2573: {
2574:   PetscInt       dim, d;

2578:   DMGetCoordinateDim(dm, &dim);
2579:   if (!dm->maxCell) {
2580:     for (d = 0; d < dim; ++d) out[d] = in[d];
2581:   } else {
2582:     for (d = 0; d < dim; ++d) {
2583:       out[d] = in[d] - dm->L[d]*floor(PetscRealPart(in[d])/dm->L[d]);
2584:     }
2585:   }
2586:   return(0);
2587: }

2591: /*
2592:   DMPlexLocalizeCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.

2594:   Input Parameters:
2595: + dm     - The DM
2596: . dim    - The spatial dimension
2597: . anchor - The anchor point, the input point can be no more than maxCell away from it
2598: - in     - The input coordinate point (dim numbers)

2600:   Output Parameter:
2601: . out - The localized coordinate point

2603:   Level: developer

2605:   Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell

2607: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeAddCoordinate()
2608: */
2609: PetscErrorCode DMPlexLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2610: {
2611:   PetscInt d;

2614:   if (!dm->maxCell) {
2615:     for (d = 0; d < dim; ++d) out[d] = in[d];
2616:   } else {
2617:     for (d = 0; d < dim; ++d) {
2618:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2619:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2620:       } else {
2621:         out[d] = in[d];
2622:       }
2623:     }
2624:   }
2625:   return(0);
2626: }
2629: PetscErrorCode DMPlexLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
2630: {
2631:   PetscInt d;

2634:   if (!dm->maxCell) {
2635:     for (d = 0; d < dim; ++d) out[d] = in[d];
2636:   } else {
2637:     for (d = 0; d < dim; ++d) {
2638:       if (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d]) {
2639:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
2640:       } else {
2641:         out[d] = in[d];
2642:       }
2643:     }
2644:   }
2645:   return(0);
2646: }

2650: /*
2651:   DMPlexLocalizeAddCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.

2653:   Input Parameters:
2654: + dm     - The DM
2655: . dim    - The spatial dimension
2656: . anchor - The anchor point, the input point can be no more than maxCell away from it
2657: . in     - The input coordinate delta (dim numbers)
2658: - out    - The input coordinate point (dim numbers)

2660:   Output Parameter:
2661: . out    - The localized coordinate in + out

2663:   Level: developer

2665:   Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell

2667: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeCoordinate()
2668: */
2669: PetscErrorCode DMPlexLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2670: {
2671:   PetscInt d;

2674:   if (!dm->maxCell) {
2675:     for (d = 0; d < dim; ++d) out[d] += in[d];
2676:   } else {
2677:     for (d = 0; d < dim; ++d) {
2678:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2679:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2680:       } else {
2681:         out[d] += in[d];
2682:       }
2683:     }
2684:   }
2685:   return(0);
2686: }

2690: /*@
2691:   DMPlexLocalizeCoordinates - If a mesh is periodic, create local coordinates for each cell

2693:   Input Parameter:
2694: . dm - The DM

2696:   Level: developer

2698: .seealso: DMPlexLocalizeCoordinate(), DMPlexLocalizeAddCoordinate()
2699: @*/
2700: PetscErrorCode DMPlexLocalizeCoordinates(DM dm)
2701: {
2702:   PetscSection   coordSection, cSection;
2703:   Vec            coordinates,  cVec;
2704:   PetscScalar   *coords, *coords2, *anchor;
2705:   PetscInt       Nc, cStart, cEnd, c, vStart, vEnd, v, dof, d, off, off2, bs, coordSize;

2710:   if (!dm->maxCell) return(0);
2711:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2712:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2713:   DMGetCoordinatesLocal(dm, &coordinates);
2714:   DMGetCoordinateSection(dm, &coordSection);
2715:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
2716:   PetscSectionSetNumFields(cSection, 1);
2717:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
2718:   PetscSectionSetFieldComponents(cSection, 0, Nc);
2719:   PetscSectionSetChart(cSection, cStart, vEnd);
2720:   for (v = vStart; v < vEnd; ++v) {
2721:     PetscSectionGetDof(coordSection, v, &dof);
2722:     PetscSectionSetDof(cSection,     v,  dof);
2723:     PetscSectionSetFieldDof(cSection, v, 0, dof);
2724:   }
2725:   for (c = cStart; c < cEnd; ++c) {
2726:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, NULL);
2727:     PetscSectionSetDof(cSection, c, dof);
2728:     PetscSectionSetFieldDof(cSection, c, 0, dof);
2729:   }
2730:   PetscSectionSetUp(cSection);
2731:   PetscSectionGetStorageSize(cSection, &coordSize);
2732:   VecCreate(PetscObjectComm((PetscObject) dm), &cVec);
2733:   VecGetBlockSize(coordinates, &bs);
2734:   VecSetBlockSize(cVec,         bs);
2735:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
2736:   VecSetType(cVec,VECSTANDARD);
2737:   VecGetArray(coordinates, &coords);
2738:   VecGetArray(cVec,        &coords2);
2739:   for (v = vStart; v < vEnd; ++v) {
2740:     PetscSectionGetDof(coordSection, v, &dof);
2741:     PetscSectionGetOffset(coordSection, v, &off);
2742:     PetscSectionGetOffset(cSection,     v, &off2);
2743:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
2744:   }
2745:   DMGetWorkArray(dm, 3, PETSC_SCALAR, &anchor);
2746:   for (c = cStart; c < cEnd; ++c) {
2747:     PetscScalar *cellCoords = NULL;
2748:     PetscInt     b;

2750:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2751:     PetscSectionGetOffset(cSection, c, &off2);
2752:     for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
2753:     for (d = 0; d < dof/bs; ++d) {DMPlexLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
2754:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2755:   }
2756:   DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);
2757:   VecRestoreArray(coordinates, &coords);
2758:   VecRestoreArray(cVec,        &coords2);
2759:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
2760:   DMSetCoordinatesLocal(dm, cVec);
2761:   VecDestroy(&cVec);
2762:   PetscSectionDestroy(&cSection);
2763:   return(0);
2764: }

2768: /*@
2769:   DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point

2771:   Not Collective

2773:   Input Parameter:
2774: . dm    - The DMPlex object

2776:   Output Parameter:
2777: . depthLabel - The DMLabel recording point depth

2779:   Level: developer

2781: .keywords: mesh, points
2782: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2783: @*/
2784: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
2785: {
2786:   DM_Plex       *mesh = (DM_Plex*) dm->data;

2792:   if (!mesh->depthLabel) {DMPlexGetLabel(dm, "depth", &mesh->depthLabel);}
2793:   *depthLabel = mesh->depthLabel;
2794:   return(0);
2795: }

2799: /*@
2800:   DMPlexGetDepth - Get the depth of the DAG representing this mesh

2802:   Not Collective

2804:   Input Parameter:
2805: . dm    - The DMPlex object

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

2810:   Level: developer

2812: .keywords: mesh, points
2813: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2814: @*/
2815: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
2816: {
2817:   DMLabel        label;
2818:   PetscInt       d = 0;

2824:   DMPlexGetDepthLabel(dm, &label);
2825:   if (label) {DMLabelGetNumValues(label, &d);}
2826:   *depth = d-1;
2827:   return(0);
2828: }

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

2835:   Not Collective

2837:   Input Parameters:
2838: + dm           - The DMPlex object
2839: - stratumValue - The requested depth

2841:   Output Parameters:
2842: + start - The first point at this depth
2843: - end   - One beyond the last point at this depth

2845:   Level: developer

2847: .keywords: mesh, points
2848: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
2849: @*/
2850: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2851: {
2852:   DMLabel        label;
2853:   PetscInt       pStart, pEnd;

2860:   DMPlexGetChart(dm, &pStart, &pEnd);
2861:   if (pStart == pEnd) return(0);
2862:   if (stratumValue < 0) {
2863:     if (start) *start = pStart;
2864:     if (end)   *end   = pEnd;
2865:     return(0);
2866:   }
2867:   DMPlexGetDepthLabel(dm, &label);
2868:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2869:   DMLabelGetStratumBounds(label, stratumValue, start, end);
2870:   return(0);
2871: }

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

2878:   Not Collective

2880:   Input Parameters:
2881: + dm           - The DMPlex object
2882: - stratumValue - The requested height

2884:   Output Parameters:
2885: + start - The first point at this height
2886: - end   - One beyond the last point at this height

2888:   Level: developer

2890: .keywords: mesh, points
2891: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
2892: @*/
2893: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2894: {
2895:   DMLabel        label;
2896:   PetscInt       depth, pStart, pEnd;

2903:   DMPlexGetChart(dm, &pStart, &pEnd);
2904:   if (pStart == pEnd) return(0);
2905:   if (stratumValue < 0) {
2906:     if (start) *start = pStart;
2907:     if (end)   *end   = pEnd;
2908:     return(0);
2909:   }
2910:   DMPlexGetDepthLabel(dm, &label);
2911:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2912:   DMLabelGetNumValues(label, &depth);
2913:   DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
2914:   return(0);
2915: }

2919: /* Set the number of dof on each point and separate by fields */
2920: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
2921: {
2922:   PetscInt      *pMax;
2923:   PetscInt       depth, pStart = 0, pEnd = 0;
2924:   PetscInt       Nf, p, d, dep, f;
2925:   PetscBool     *isFE;

2929:   PetscMalloc1(numFields, &isFE);
2930:   DMGetNumFields(dm, &Nf);
2931:   for (f = 0; f < numFields; ++f) {
2932:     PetscObject  obj;
2933:     PetscClassId id;

2935:     isFE[f] = PETSC_FALSE;
2936:     if (f >= Nf) continue;
2937:     DMGetField(dm, f, &obj);
2938:     PetscObjectGetClassId(obj, &id);
2939:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
2940:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
2941:   }
2942:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
2943:   if (numFields > 0) {
2944:     PetscSectionSetNumFields(*section, numFields);
2945:     if (numComp) {
2946:       for (f = 0; f < numFields; ++f) {
2947:         PetscSectionSetFieldComponents(*section, f, numComp[f]);
2948:       }
2949:     }
2950:   }
2951:   DMPlexGetChart(dm, &pStart, &pEnd);
2952:   PetscSectionSetChart(*section, pStart, pEnd);
2953:   DMPlexGetDepth(dm, &depth);
2954:   PetscMalloc1(depth+1,&pMax);
2955:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2956:   for (dep = 0; dep <= depth; ++dep) {
2957:     d    = dim == depth ? dep : (!dep ? 0 : dim);
2958:     DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
2959:     pMax[dep] = pMax[dep] < 0 ? pEnd : pMax[dep];
2960:     for (p = pStart; p < pEnd; ++p) {
2961:       PetscInt tot = 0;

2963:       for (f = 0; f < numFields; ++f) {
2964:         if (isFE[f] && p >= pMax[dep]) continue;
2965:         PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
2966:         tot += numDof[f*(dim+1)+d];
2967:       }
2968:       PetscSectionSetDof(*section, p, tot);
2969:     }
2970:   }
2971:   PetscFree(pMax);
2972:   PetscFree(isFE);
2973:   return(0);
2974: }

2978: /* Set the number of dof on each point and separate by fields
2979:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
2980: */
2981: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2982: {
2983:   PetscInt       numFields;
2984:   PetscInt       bc;
2985:   PetscSection   aSec;

2989:   PetscSectionGetNumFields(section, &numFields);
2990:   for (bc = 0; bc < numBC; ++bc) {
2991:     PetscInt        field = 0;
2992:     const PetscInt *comp;
2993:     const PetscInt *idx;
2994:     PetscInt        Nc = -1, n, i;

2996:     if (numFields) field = bcField[bc];
2997:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2998:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2999:     ISGetLocalSize(bcPoints[bc], &n);
3000:     ISGetIndices(bcPoints[bc], &idx);
3001:     for (i = 0; i < n; ++i) {
3002:       const PetscInt p = idx[i];
3003:       PetscInt       numConst;

3005:       if (numFields) {
3006:         PetscSectionGetFieldDof(section, p, field, &numConst);
3007:       } else {
3008:         PetscSectionGetDof(section, p, &numConst);
3009:       }
3010:       /* If Nc < 0, constrain every dof on the point */
3011:       if (Nc > 0) numConst = PetscMin(numConst, Nc);
3012:       if (numFields) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
3013:       PetscSectionAddConstraintDof(section, p, numConst);
3014:     }
3015:     ISRestoreIndices(bcPoints[bc], &idx);
3016:     if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
3017:   }
3018:   DMPlexGetAnchors(dm, &aSec, NULL);
3019:   if (aSec) {
3020:     PetscInt aStart, aEnd, a;

3022:     PetscSectionGetChart(aSec, &aStart, &aEnd);
3023:     for (a = aStart; a < aEnd; a++) {
3024:       PetscInt dof, f;

3026:       PetscSectionGetDof(aSec, a, &dof);
3027:       if (dof) {
3028:         /* if there are point-to-point constraints, then all dofs are constrained */
3029:         PetscSectionGetDof(section, a, &dof);
3030:         PetscSectionSetConstraintDof(section, a, dof);
3031:         for (f = 0; f < numFields; f++) {
3032:           PetscSectionGetFieldDof(section, a, f, &dof);
3033:           PetscSectionSetFieldConstraintDof(section, a, f, dof);
3034:         }
3035:       }
3036:     }
3037:   }
3038:   return(0);
3039: }

3043: /* Set the constrained field indices on each point
3044:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
3045: */
3046: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
3047: {
3048:   PetscSection   aSec;
3049:   PetscInt      *indices;
3050:   PetscInt       numFields, maxDof, pStart, pEnd, p, bc, f, d;

3054:   PetscSectionGetNumFields(section, &numFields);
3055:   if (!numFields) return(0);
3056:   /* Initialize all field indices to -1 */
3057:   PetscSectionGetChart(section, &pStart, &pEnd);
3058:   PetscSectionGetMaxDof(section, &maxDof);
3059:   PetscMalloc1(maxDof, &indices);
3060:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
3061:   for (p = pStart; p < pEnd; ++p) for (f = 0; f < numFields; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
3062:   /* Handle BC constraints */
3063:   for (bc = 0; bc < numBC; ++bc) {
3064:     const PetscInt  field = bcField[bc];
3065:     const PetscInt *comp, *idx;
3066:     PetscInt        Nc = -1, n, i;

3068:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
3069:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
3070:     ISGetLocalSize(bcPoints[bc], &n);
3071:     ISGetIndices(bcPoints[bc], &idx);
3072:     for (i = 0; i < n; ++i) {
3073:       const PetscInt  p = idx[i];
3074:       const PetscInt *find;
3075:       PetscInt        fcdof, c;

3077:       PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
3078:       if (Nc < 0) {
3079:         for (d = 0; d < fcdof; ++d) indices[d] = d;
3080:       } else {
3081:         PetscSectionGetFieldConstraintIndices(section, p, field, &find);
3082:         for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
3083:         for (c = 0; c < Nc; ++c) indices[d+c] = comp[c];
3084:         PetscSortInt(d+Nc, indices);
3085:         for (c = d+Nc; c < fcdof; ++c) indices[c] = -1;
3086:       }
3087:       PetscSectionSetFieldConstraintIndices(section, p, field, indices);
3088:     }
3089:     if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
3090:     ISRestoreIndices(bcPoints[bc], &idx);
3091:   }
3092:   /* Handle anchors */
3093:   DMPlexGetAnchors(dm, &aSec, NULL);
3094:   if (aSec) {
3095:     PetscInt aStart, aEnd, a;

3097:     for (d = 0; d < maxDof; ++d) indices[d] = d;
3098:     PetscSectionGetChart(aSec, &aStart, &aEnd);
3099:     for (a = aStart; a < aEnd; a++) {
3100:       PetscInt dof, fdof, f;

3102:       PetscSectionGetDof(aSec, a, &dof);
3103:       if (dof) {
3104:         /* if there are point-to-point constraints, then all dofs are constrained */
3105:         for (f = 0; f < numFields; f++) {
3106:           PetscSectionGetFieldDof(section, a, f, &fdof);
3107:           PetscSectionSetFieldConstraintIndices(section, a, f, indices);
3108:         }
3109:       }
3110:     }
3111:   }
3112:   PetscFree(indices);
3113:   return(0);
3114: }

3118: /* Set the constrained indices on each point */
3119: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
3120: {
3121:   PetscInt      *indices;
3122:   PetscInt       numFields, maxDof, pStart, pEnd, p, f, d;

3126:   PetscSectionGetNumFields(section, &numFields);
3127:   PetscSectionGetMaxDof(section, &maxDof);
3128:   PetscSectionGetChart(section, &pStart, &pEnd);
3129:   PetscMalloc1(maxDof, &indices);
3130:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
3131:   for (p = pStart; p < pEnd; ++p) {
3132:     PetscInt cdof, d;

3134:     PetscSectionGetConstraintDof(section, p, &cdof);
3135:     if (cdof) {
3136:       if (numFields) {
3137:         PetscInt numConst = 0, foff = 0;

3139:         for (f = 0; f < numFields; ++f) {
3140:           const PetscInt *find;
3141:           PetscInt        fcdof, fdof;

3143:           PetscSectionGetFieldDof(section, p, f, &fdof);
3144:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
3145:           /* Change constraint numbering from field component to local dof number */
3146:           PetscSectionGetFieldConstraintIndices(section, p, f, &find);
3147:           for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
3148:           numConst += fcdof;
3149:           foff     += fdof;
3150:         }
3151:         if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
3152:       } else {
3153:         for (d = 0; d < cdof; ++d) indices[d] = d;
3154:       }
3155:       PetscSectionSetConstraintIndices(section, p, indices);
3156:     }
3157:   }
3158:   PetscFree(indices);
3159:   return(0);
3160: }

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

3167:   Not Collective

3169:   Input Parameters:
3170: + dm        - The DMPlex object
3171: . dim       - The spatial dimension of the problem
3172: . numFields - The number of fields in the problem
3173: . numComp   - An array of size numFields that holds the number of components for each field
3174: . numDof    - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
3175: . numBC     - The number of boundary conditions
3176: . bcField   - An array of size numBC giving the field number for each boundry condition
3177: . bcComps   - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
3178: . bcPoints  - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
3179: - perm      - Optional permutation of the chart, or NULL

3181:   Output Parameter:
3182: . section - The PetscSection object

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

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

3189:   Level: developer

3191:   Fortran Notes:
3192:   A Fortran 90 version is available as DMPlexCreateSectionF90()

3194: .keywords: mesh, elements
3195: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
3196: @*/
3197: 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)
3198: {
3199:   PetscSection   aSec;

3203:   DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
3204:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
3205:   if (perm) {PetscSectionSetPermutation(*section, perm);}
3206:   PetscSectionSetUp(*section);
3207:   DMPlexGetAnchors(dm,&aSec,NULL);
3208:   if (numBC || aSec) {
3209:     DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
3210:     DMPlexCreateSectionBCIndices(dm, *section);
3211:   }
3212:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
3213:   return(0);
3214: }

3218: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
3219: {
3220:   PetscSection   section;

3224:   DMClone(dm, cdm);
3225:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
3226:   DMSetDefaultSection(*cdm, section);
3227:   PetscSectionDestroy(&section);
3228:   return(0);
3229: }

3233: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
3234: {
3235:   DM_Plex *mesh = (DM_Plex*) dm->data;

3239:   if (section) *section = mesh->coneSection;
3240:   return(0);
3241: }

3245: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
3246: {
3247:   DM_Plex *mesh = (DM_Plex*) dm->data;

3251:   if (section) *section = mesh->supportSection;
3252:   return(0);
3253: }

3257: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
3258: {
3259:   DM_Plex *mesh = (DM_Plex*) dm->data;

3263:   if (cones) *cones = mesh->cones;
3264:   return(0);
3265: }

3269: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
3270: {
3271:   DM_Plex *mesh = (DM_Plex*) dm->data;

3275:   if (coneOrientations) *coneOrientations = mesh->coneOrientations;
3276:   return(0);
3277: }

3279: /******************************** FEM Support **********************************/

3283: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3284: {
3285:   PetscScalar    *array, *vArray;
3286:   const PetscInt *cone, *coneO;
3287:   PetscInt        pStart, pEnd, p, numPoints, size = 0, offset = 0;
3288:   PetscErrorCode  ierr;

3291:   PetscSectionGetChart(section, &pStart, &pEnd);
3292:   DMPlexGetConeSize(dm, point, &numPoints);
3293:   DMPlexGetCone(dm, point, &cone);
3294:   DMPlexGetConeOrientation(dm, point, &coneO);
3295:   if (!values || !*values) {
3296:     if ((point >= pStart) && (point < pEnd)) {
3297:       PetscInt dof;

3299:       PetscSectionGetDof(section, point, &dof);
3300:       size += dof;
3301:     }
3302:     for (p = 0; p < numPoints; ++p) {
3303:       const PetscInt cp = cone[p];
3304:       PetscInt       dof;

3306:       if ((cp < pStart) || (cp >= pEnd)) continue;
3307:       PetscSectionGetDof(section, cp, &dof);
3308:       size += dof;
3309:     }
3310:     if (!values) {
3311:       if (csize) *csize = size;
3312:       return(0);
3313:     }
3314:     DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
3315:   } else {
3316:     array = *values;
3317:   }
3318:   size = 0;
3319:   VecGetArray(v, &vArray);
3320:   if ((point >= pStart) && (point < pEnd)) {
3321:     PetscInt     dof, off, d;
3322:     PetscScalar *varr;

3324:     PetscSectionGetDof(section, point, &dof);
3325:     PetscSectionGetOffset(section, point, &off);
3326:     varr = &vArray[off];
3327:     for (d = 0; d < dof; ++d, ++offset) {
3328:       array[offset] = varr[d];
3329:     }
3330:     size += dof;
3331:   }
3332:   for (p = 0; p < numPoints; ++p) {
3333:     const PetscInt cp = cone[p];
3334:     PetscInt       o  = coneO[p];
3335:     PetscInt       dof, off, d;
3336:     PetscScalar   *varr;

3338:     if ((cp < pStart) || (cp >= pEnd)) continue;
3339:     PetscSectionGetDof(section, cp, &dof);
3340:     PetscSectionGetOffset(section, cp, &off);
3341:     varr = &vArray[off];
3342:     if (o >= 0) {
3343:       for (d = 0; d < dof; ++d, ++offset) {
3344:         array[offset] = varr[d];
3345:       }
3346:     } else {
3347:       for (d = dof-1; d >= 0; --d, ++offset) {
3348:         array[offset] = varr[d];
3349:       }
3350:     }
3351:     size += dof;
3352:   }
3353:   VecRestoreArray(v, &vArray);
3354:   if (!*values) {
3355:     if (csize) *csize = size;
3356:     *values = array;
3357:   } else {
3358:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3359:     *csize = size;
3360:   }
3361:   return(0);
3362: }

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

3372:   *size = 0;
3373:   for (p = 0; p < numPoints*2; p += 2) {
3374:     const PetscInt point = points[p];
3375:     const PetscInt o     = points[p+1];
3376:     PetscInt       dof, off, d;
3377:     const PetscScalar *varr;

3379:     PetscSectionGetDof(section, point, &dof);
3380:     PetscSectionGetOffset(section, point, &off);
3381:     varr = &vArray[off];
3382:     if (o >= 0) {
3383:       for (d = 0; d < dof; ++d, ++offset)    array[offset] = varr[d];
3384:     } else {
3385:       for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
3386:     }
3387:   }
3388:   *size = offset;
3389:   return(0);
3390: }

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

3400:   *size = 0;
3401:   for (f = 0; f < numFields; ++f) {
3402:     PetscInt fcomp, p;

3404:     PetscSectionGetFieldComponents(section, f, &fcomp);
3405:     for (p = 0; p < numPoints*2; p += 2) {
3406:       const PetscInt point = points[p];
3407:       const PetscInt o     = points[p+1];
3408:       PetscInt       fdof, foff, d, c;
3409:       const PetscScalar *varr;

3411:       PetscSectionGetFieldDof(section, point, f, &fdof);
3412:       PetscSectionGetFieldOffset(section, point, f, &foff);
3413:       varr = &vArray[foff];
3414:       if (o >= 0) {
3415:         for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
3416:       } else {
3417:         for (d = fdof/fcomp-1; d >= 0; --d) {
3418:           for (c = 0; c < fcomp; ++c, ++offset) {
3419:             array[offset] = varr[d*fcomp+c];
3420:           }
3421:         }
3422:       }
3423:     }
3424:   }
3425:   *size = offset;
3426:   return(0);
3427: }

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

3434:   Not collective

3436:   Input Parameters:
3437: + dm - The DM
3438: . section - The section describing the layout in v, or NULL to use the default section
3439: . v - The local vector
3440: - point - The sieve point in the DM

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

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

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

3452:   Level: intermediate

3454: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3455: @*/
3456: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3457: {
3458:   PetscSection    clSection;
3459:   IS              clPoints;
3460:   PetscScalar    *array, *vArray;
3461:   PetscInt       *points = NULL;
3462:   const PetscInt *clp;
3463:   PetscInt        depth, numFields, numPoints, size;
3464:   PetscErrorCode  ierr;

3468:   if (!section) {DMGetDefaultSection(dm, &section);}
3471:   DMPlexGetDepth(dm, &depth);
3472:   PetscSectionGetNumFields(section, &numFields);
3473:   if (depth == 1 && numFields < 2) {
3474:     DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
3475:     return(0);
3476:   }
3477:   /* Get points */
3478:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3479:   if (!clPoints) {
3480:     PetscInt pStart, pEnd, p, q;

3482:     PetscSectionGetChart(section, &pStart, &pEnd);
3483:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3484:     /* Compress out points not in the section */
3485:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3486:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3487:         points[q*2]   = points[p];
3488:         points[q*2+1] = points[p+1];
3489:         ++q;
3490:       }
3491:     }
3492:     numPoints = q;
3493:   } else {
3494:     PetscInt dof, off;

3496:     PetscSectionGetDof(clSection, point, &dof);
3497:     PetscSectionGetOffset(clSection, point, &off);
3498:     ISGetIndices(clPoints, &clp);
3499:     numPoints = dof/2;
3500:     points    = (PetscInt *) &clp[off];
3501:   }
3502:   /* Get array */
3503:   if (!values || !*values) {
3504:     PetscInt asize = 0, dof, p;

3506:     for (p = 0; p < numPoints*2; p += 2) {
3507:       PetscSectionGetDof(section, points[p], &dof);
3508:       asize += dof;
3509:     }
3510:     if (!values) {
3511:       if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3512:       else           {ISRestoreIndices(clPoints, &clp);}
3513:       if (csize) *csize = asize;
3514:       return(0);
3515:     }
3516:     DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
3517:   } else {
3518:     array = *values;
3519:   }
3520:   VecGetArray(v, &vArray);
3521:   /* Get values */
3522:   if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);}
3523:   else               {DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);}
3524:   /* Cleanup points */
3525:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3526:   else           {ISRestoreIndices(clPoints, &clp);}
3527:   /* Cleanup array */
3528:   VecRestoreArray(v, &vArray);
3529:   if (!*values) {
3530:     if (csize) *csize = size;
3531:     *values = array;
3532:   } else {
3533:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3534:     *csize = size;
3535:   }
3536:   return(0);
3537: }

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

3544:   Not collective

3546:   Input Parameters:
3547: + dm - The DM
3548: . section - The section describing the layout in v, or NULL to use the default section
3549: . v - The local vector
3550: . point - The sieve point in the DM
3551: . csize - The number of values in the closure, or NULL
3552: - values - The array of values, which is a borrowed array and should not be freed

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

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

3560:   Level: intermediate

3562: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3563: @*/
3564: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3565: {
3566:   PetscInt       size = 0;

3570:   /* Should work without recalculating size */
3571:   DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
3572:   return(0);
3573: }

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

3580: PETSC_STATIC_INLINE PetscErrorCode updatePoint_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3581: {
3582:   PetscInt        cdof;   /* The number of constraints on this point */
3583:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3584:   PetscScalar    *a;
3585:   PetscInt        off, cind = 0, k;
3586:   PetscErrorCode  ierr;

3589:   PetscSectionGetConstraintDof(section, point, &cdof);
3590:   PetscSectionGetOffset(section, point, &off);
3591:   a    = &array[off];
3592:   if (!cdof || setBC) {
3593:     if (orientation >= 0) {
3594:       for (k = 0; k < dof; ++k) {
3595:         fuse(&a[k], values[k]);
3596:       }
3597:     } else {
3598:       for (k = 0; k < dof; ++k) {
3599:         fuse(&a[k], values[dof-k-1]);
3600:       }
3601:     }
3602:   } else {
3603:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3604:     if (orientation >= 0) {
3605:       for (k = 0; k < dof; ++k) {
3606:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3607:         fuse(&a[k], values[k]);
3608:       }
3609:     } else {
3610:       for (k = 0; k < dof; ++k) {
3611:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3612:         fuse(&a[k], values[dof-k-1]);
3613:       }
3614:     }
3615:   }
3616:   return(0);
3617: }

3621: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3622: {
3623:   PetscInt        cdof;   /* The number of constraints on this point */
3624:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3625:   PetscScalar    *a;
3626:   PetscInt        off, cind = 0, k;
3627:   PetscErrorCode  ierr;

3630:   PetscSectionGetConstraintDof(section, point, &cdof);
3631:   PetscSectionGetOffset(section, point, &off);
3632:   a    = &array[off];
3633:   if (cdof) {
3634:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3635:     if (orientation >= 0) {
3636:       for (k = 0; k < dof; ++k) {
3637:         if ((cind < cdof) && (k == cdofs[cind])) {
3638:           fuse(&a[k], values[k]);
3639:           ++cind;
3640:         }
3641:       }
3642:     } else {
3643:       for (k = 0; k < dof; ++k) {
3644:         if ((cind < cdof) && (k == cdofs[cind])) {
3645:           fuse(&a[k], values[dof-k-1]);
3646:           ++cind;
3647:         }
3648:       }
3649:     }
3650:   }
3651:   return(0);
3652: }

3656: PETSC_STATIC_INLINE PetscErrorCode updatePointFields_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, const PetscScalar values[], PetscInt *offset, PetscScalar array[])
3657: {
3658:   PetscScalar    *a;
3659:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3660:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3661:   PetscInt        cind = 0, k, c;
3662:   PetscErrorCode  ierr;

3665:   PetscSectionGetFieldDof(section, point, f, &fdof);
3666:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3667:   PetscSectionGetFieldOffset(section, point, f, &foff);
3668:   a    = &array[foff];
3669:   if (!fcdof || setBC) {
3670:     if (o >= 0) {
3671:       for (k = 0; k < fdof; ++k) fuse(&a[k], values[foffset+k]);
3672:     } else {
3673:       for (k = fdof/fcomp-1; k >= 0; --k) {
3674:         for (c = 0; c < fcomp; ++c) {
3675:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3676:         }
3677:       }
3678:     }
3679:   } else {
3680:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3681:     if (o >= 0) {
3682:       for (k = 0; k < fdof; ++k) {
3683:         if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
3684:         fuse(&a[k], values[foffset+k]);
3685:       }
3686:     } else {
3687:       for (k = fdof/fcomp-1; k >= 0; --k) {
3688:         for (c = 0; c < fcomp; ++c) {
3689:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
3690:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3691:         }
3692:       }
3693:     }
3694:   }
3695:   *offset += fdof;
3696:   return(0);
3697: }

3701: PETSC_STATIC_INLINE PetscErrorCode updatePointFieldsBC_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), const PetscScalar values[], PetscInt *offset, PetscScalar array[])
3702: {
3703:   PetscScalar    *a;
3704:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3705:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3706:   PetscInt        cind = 0, k, c;
3707:   PetscErrorCode  ierr;

3710:   PetscSectionGetFieldDof(section, point, f, &fdof);
3711:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3712:   PetscSectionGetFieldOffset(section, point, f, &foff);
3713:   a    = &array[foff];
3714:   if (fcdof) {
3715:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3716:     if (o >= 0) {
3717:       for (k = 0; k < fdof; ++k) {
3718:         if ((cind < fcdof) && (k == fcdofs[cind])) {
3719:           fuse(&a[k], values[foffset+k]);
3720:           ++cind;
3721:         }
3722:       }
3723:     } else {
3724:       for (k = fdof/fcomp-1; k >= 0; --k) {
3725:         for (c = 0; c < fcomp; ++c) {
3726:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {
3727:             fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3728:             ++cind;
3729:           }
3730:         }
3731:       }
3732:     }
3733:   }
3734:   *offset += fdof;
3735:   return(0);
3736: }

3740: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3741: {
3742:   PetscScalar    *array;
3743:   const PetscInt *cone, *coneO;
3744:   PetscInt        pStart, pEnd, p, numPoints, off, dof;
3745:   PetscErrorCode  ierr;

3748:   PetscSectionGetChart(section, &pStart, &pEnd);
3749:   DMPlexGetConeSize(dm, point, &numPoints);
3750:   DMPlexGetCone(dm, point, &cone);
3751:   DMPlexGetConeOrientation(dm, point, &coneO);
3752:   VecGetArray(v, &array);
3753:   for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
3754:     const PetscInt cp = !p ? point : cone[p-1];
3755:     const PetscInt o  = !p ? 0     : coneO[p-1];

3757:     if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
3758:     PetscSectionGetDof(section, cp, &dof);
3759:     /* ADD_VALUES */
3760:     {
3761:       const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3762:       PetscScalar    *a;
3763:       PetscInt        cdof, coff, cind = 0, k;

3765:       PetscSectionGetConstraintDof(section, cp, &cdof);
3766:       PetscSectionGetOffset(section, cp, &coff);
3767:       a    = &array[coff];
3768:       if (!cdof) {
3769:         if (o >= 0) {
3770:           for (k = 0; k < dof; ++k) {
3771:             a[k] += values[off+k];
3772:           }
3773:         } else {
3774:           for (k = 0; k < dof; ++k) {
3775:             a[k] += values[off+dof-k-1];
3776:           }
3777:         }
3778:       } else {
3779:         PetscSectionGetConstraintIndices(section, cp, &cdofs);
3780:         if (o >= 0) {
3781:           for (k = 0; k < dof; ++k) {
3782:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3783:             a[k] += values[off+k];
3784:           }
3785:         } else {
3786:           for (k = 0; k < dof; ++k) {
3787:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3788:             a[k] += values[off+dof-k-1];
3789:           }
3790:         }
3791:       }
3792:     }
3793:   }
3794:   VecRestoreArray(v, &array);
3795:   return(0);
3796: }

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

3803:   Not collective

3805:   Input Parameters:
3806: + dm - The DM
3807: . section - The section describing the layout in v, or NULL to use the default section
3808: . v - The local vector
3809: . point - The sieve point in the DM
3810: . values - The array of values
3811: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

3816:   Level: intermediate

3818: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
3819: @*/
3820: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3821: {
3822:   PetscSection    clSection;
3823:   IS              clPoints;
3824:   PetscScalar    *array;
3825:   PetscInt       *points = NULL;
3826:   const PetscInt *clp;
3827:   PetscInt        depth, numFields, numPoints, p;
3828:   PetscErrorCode  ierr;

3832:   if (!section) {DMGetDefaultSection(dm, &section);}
3835:   DMPlexGetDepth(dm, &depth);
3836:   PetscSectionGetNumFields(section, &numFields);
3837:   if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
3838:     DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
3839:     return(0);
3840:   }
3841:   /* Get points */
3842:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3843:   if (!clPoints) {
3844:     PetscInt pStart, pEnd, q;

3846:     PetscSectionGetChart(section, &pStart, &pEnd);
3847:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3848:     /* Compress out points not in the section */
3849:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3850:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3851:         points[q*2]   = points[p];
3852:         points[q*2+1] = points[p+1];
3853:         ++q;
3854:       }
3855:     }
3856:     numPoints = q;
3857:   } else {
3858:     PetscInt dof, off;

3860:     PetscSectionGetDof(clSection, point, &dof);
3861:     PetscSectionGetOffset(clSection, point, &off);
3862:     ISGetIndices(clPoints, &clp);
3863:     numPoints = dof/2;
3864:     points    = (PetscInt *) &clp[off];
3865:   }
3866:   /* Get array */
3867:   VecGetArray(v, &array);
3868:   /* Get values */
3869:   if (numFields > 0) {
3870:     PetscInt offset = 0, fcomp, f;
3871:     for (f = 0; f < numFields; ++f) {
3872:       PetscSectionGetFieldComponents(section, f, &fcomp);
3873:       switch (mode) {
3874:       case INSERT_VALUES:
3875:         for (p = 0; p < numPoints*2; p += 2) {
3876:           const PetscInt point = points[p];
3877:           const PetscInt o     = points[p+1];
3878:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3879:         } break;
3880:       case INSERT_ALL_VALUES:
3881:         for (p = 0; p < numPoints*2; p += 2) {
3882:           const PetscInt point = points[p];
3883:           const PetscInt o     = points[p+1];
3884:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3885:         } break;
3886:       case INSERT_BC_VALUES:
3887:         for (p = 0; p < numPoints*2; p += 2) {
3888:           const PetscInt point = points[p];
3889:           const PetscInt o     = points[p+1];
3890:           updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3891:         } break;
3892:       case ADD_VALUES:
3893:         for (p = 0; p < numPoints*2; p += 2) {
3894:           const PetscInt point = points[p];
3895:           const PetscInt o     = points[p+1];
3896:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3897:         } break;
3898:       case ADD_ALL_VALUES:
3899:         for (p = 0; p < numPoints*2; p += 2) {
3900:           const PetscInt point = points[p];
3901:           const PetscInt o     = points[p+1];
3902:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3903:         } break;
3904:       default:
3905:         SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3906:       }
3907:     }
3908:   } else {
3909:     PetscInt dof, off;

3911:     switch (mode) {
3912:     case INSERT_VALUES:
3913:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3914:         PetscInt o = points[p+1];
3915:         PetscSectionGetDof(section, points[p], &dof);
3916:         updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
3917:       } break;
3918:     case INSERT_ALL_VALUES:
3919:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3920:         PetscInt o = points[p+1];
3921:         PetscSectionGetDof(section, points[p], &dof);
3922:         updatePoint_private(section, points[p], dof, insert, PETSC_TRUE,  o, &values[off], array);
3923:       } break;
3924:     case INSERT_BC_VALUES:
3925:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3926:         PetscInt o = points[p+1];
3927:         PetscSectionGetDof(section, points[p], &dof);
3928:         updatePointBC_private(section, points[p], dof, insert,  o, &values[off], array);
3929:       } break;
3930:     case ADD_VALUES:
3931:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3932:         PetscInt o = points[p+1];
3933:         PetscSectionGetDof(section, points[p], &dof);
3934:         updatePoint_private(section, points[p], dof, add,    PETSC_FALSE, o, &values[off], array);
3935:       } break;
3936:     case ADD_ALL_VALUES:
3937:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3938:         PetscInt o = points[p+1];
3939:         PetscSectionGetDof(section, points[p], &dof);
3940:         updatePoint_private(section, points[p], dof, add,    PETSC_TRUE,  o, &values[off], array);
3941:       } break;
3942:     default:
3943:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3944:     }
3945:   }
3946:   /* Cleanup points */
3947:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3948:   else           {ISRestoreIndices(clPoints, &clp);}
3949:   /* Cleanup array */
3950:   VecRestoreArray(v, &array);
3951:   return(0);
3952: }

3956: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
3957: {
3958:   PetscSection    clSection;
3959:   IS              clPoints;
3960:   PetscScalar    *array;
3961:   PetscInt       *points = NULL;
3962:   const PetscInt *clp;
3963:   PetscInt        numFields, numPoints, p;
3964:   PetscInt        offset = 0, fcomp, f;
3965:   PetscErrorCode  ierr;

3969:   if (!section) {DMGetDefaultSection(dm, &section);}
3972:   PetscSectionGetNumFields(section, &numFields);
3973:   /* Get points */
3974:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3975:   if (!clPoints) {
3976:     PetscInt pStart, pEnd, q;

3978:     PetscSectionGetChart(section, &pStart, &pEnd);
3979:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3980:     /* Compress out points not in the section */
3981:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3982:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3983:         points[q*2]   = points[p];
3984:         points[q*2+1] = points[p+1];
3985:         ++q;
3986:       }
3987:     }
3988:     numPoints = q;
3989:   } else {
3990:     PetscInt dof, off;

3992:     PetscSectionGetDof(clSection, point, &dof);
3993:     PetscSectionGetOffset(clSection, point, &off);
3994:     ISGetIndices(clPoints, &clp);
3995:     numPoints = dof/2;
3996:     points    = (PetscInt *) &clp[off];
3997:   }
3998:   /* Get array */
3999:   VecGetArray(v, &array);
4000:   /* Get values */
4001:   for (f = 0; f < numFields; ++f) {
4002:     PetscSectionGetFieldComponents(section, f, &fcomp);
4003:     if (!fieldActive[f]) {
4004:       for (p = 0; p < numPoints*2; p += 2) {
4005:         PetscInt fdof;
4006:         PetscSectionGetFieldDof(section, points[p], f, &fdof);
4007:         offset += fdof;
4008:       }
4009:       continue;
4010:     }
4011:     switch (mode) {
4012:     case INSERT_VALUES:
4013:       for (p = 0; p < numPoints*2; p += 2) {
4014:         const PetscInt point = points[p];
4015:         const PetscInt o     = points[p+1];
4016:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
4017:       } break;
4018:     case INSERT_ALL_VALUES:
4019:       for (p = 0; p < numPoints*2; p += 2) {
4020:         const PetscInt point = points[p];
4021:         const PetscInt o     = points[p+1];
4022:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
4023:         } break;
4024:     case INSERT_BC_VALUES:
4025:       for (p = 0; p < numPoints*2; p += 2) {
4026:         const PetscInt point = points[p];
4027:         const PetscInt o     = points[p+1];
4028:         updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
4029:       } break;
4030:     case ADD_VALUES:
4031:       for (p = 0; p < numPoints*2; p += 2) {
4032:         const PetscInt point = points[p];
4033:         const PetscInt o     = points[p+1];
4034:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
4035:       } break;
4036:     case ADD_ALL_VALUES:
4037:       for (p = 0; p < numPoints*2; p += 2) {
4038:         const PetscInt point = points[p];
4039:         const PetscInt o     = points[p+1];
4040:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
4041:       } break;
4042:     default:
4043:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
4044:     }
4045:   }
4046:   /* Cleanup points */
4047:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4048:   else           {ISRestoreIndices(clPoints, &clp);}
4049:   /* Cleanup array */
4050:   VecRestoreArray(v, &array);
4051:   return(0);
4052: }

4056: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
4057: {
4058:   PetscMPIInt    rank;
4059:   PetscInt       i, j;

4063:   MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
4064:   PetscViewerASCIIPrintf(viewer, "[%D]mat for sieve point %D\n", rank, point);
4065:   for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
4066:   for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
4067:   numCIndices = numCIndices ? numCIndices : numRIndices;
4068:   for (i = 0; i < numRIndices; i++) {
4069:     PetscViewerASCIIPrintf(viewer, "[%d]", rank);
4070:     for (j = 0; j < numCIndices; j++) {
4071: #if defined(PETSC_USE_COMPLEX)
4072:       PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
4073: #else
4074:       PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
4075: #endif
4076:     }
4077:     PetscViewerASCIIPrintf(viewer, "\n");
4078:   }
4079:   return(0);
4080: }

4084: /* . off - The global offset of this point */
4085: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
4086: {
4087:   PetscInt        dof;    /* The number of unknowns on this point */
4088:   PetscInt        cdof;   /* The number of constraints on this point */
4089:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
4090:   PetscInt        cind = 0, k;
4091:   PetscErrorCode  ierr;

4094:   PetscSectionGetDof(section, point, &dof);
4095:   PetscSectionGetConstraintDof(section, point, &cdof);
4096:   if (!cdof || setBC) {
4097:     if (orientation >= 0) {
4098:       for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
4099:     } else {
4100:       for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
4101:     }
4102:   } else {
4103:     PetscSectionGetConstraintIndices(section, point, &cdofs);
4104:     if (orientation >= 0) {
4105:       for (k = 0; k < dof; ++k) {
4106:         if ((cind < cdof) && (k == cdofs[cind])) {
4107:           /* Insert check for returning constrained indices */
4108:           indices[*loff+k] = -(off+k+1);
4109:           ++cind;
4110:         } else {
4111:           indices[*loff+k] = off+k-cind;
4112:         }
4113:       }
4114:     } else {
4115:       for (k = 0; k < dof; ++k) {
4116:         if ((cind < cdof) && (k == cdofs[cind])) {
4117:           /* Insert check for returning constrained indices */
4118:           indices[*loff+dof-k-1] = -(off+k+1);
4119:           ++cind;
4120:         } else {
4121:           indices[*loff+dof-k-1] = off+k-cind;
4122:         }
4123:       }
4124:     }
4125:   }
4126:   *loff += dof;
4127:   return(0);
4128: }

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

4139:   PetscSectionGetNumFields(section, &numFields);
4140:   for (f = 0, foff = 0; f < numFields; ++f) {
4141:     PetscInt        fdof, fcomp, cfdof;
4142:     const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4143:     PetscInt        cind = 0, k, c;

4145:     PetscSectionGetFieldComponents(section, f, &fcomp);
4146:     PetscSectionGetFieldDof(section, point, f, &fdof);
4147:     PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
4148:     if (!cfdof || setBC) {
4149:       if (orientation >= 0) {
4150:         for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
4151:       } else {
4152:         for (k = fdof/fcomp-1; k >= 0; --k) {
4153:           for (c = 0; c < fcomp; ++c) {
4154:             indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
4155:           }
4156:         }
4157:       }
4158:     } else {
4159:       PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
4160:       if (orientation >= 0) {
4161:         for (k = 0; k < fdof; ++k) {
4162:           if ((cind < cfdof) && (k == fcdofs[cind])) {
4163:             indices[foffs[f]+k] = -(off+foff+k+1);
4164:             ++cind;
4165:           } else {
4166:             indices[foffs[f]+k] = off+foff+k-cind;
4167:           }
4168:         }
4169:       } else {
4170:         for (k = fdof/fcomp-1; k >= 0; --k) {
4171:           for (c = 0; c < fcomp; ++c) {
4172:             if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
4173:               indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
4174:               ++cind;
4175:             } else {
4176:               indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
4177:             }
4178:           }
4179:         }
4180:       }
4181:     }
4182:     foff     += fdof - cfdof;
4183:     foffs[f] += fdof;
4184:   }
4185:   return(0);
4186: }

4190: static 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[])
4191: {
4192:   Mat             cMat;
4193:   PetscSection    aSec, cSec;
4194:   IS              aIS;
4195:   PetscInt        aStart = -1, aEnd = -1;
4196:   const PetscInt  *anchors;
4197:   PetscInt        numFields, f, p, q, newP = 0;
4198:   PetscInt        newNumPoints = 0, newNumIndices = 0;
4199:   PetscInt        *newPoints, *indices, *newIndices;
4200:   PetscInt        maxAnchor, maxDof;
4201:   PetscInt        newOffsets[32];
4202:   PetscInt        *pointMatOffsets[32];
4203:   PetscInt        *newPointOffsets[32];
4204:   PetscScalar     *pointMat[32];
4205:   PetscScalar     *newValues,*tmpValues;
4206:   PetscBool       anyConstrained = PETSC_FALSE;
4207:   PetscErrorCode  ierr;

4212:   PetscSectionGetNumFields(section, &numFields);

4214:   DMPlexGetAnchors(dm,&aSec,&aIS);
4215:   /* if there are point-to-point constraints */
4216:   if (aSec) {
4217:     PetscMemzero(newOffsets, 32 * sizeof(PetscInt));
4218:     ISGetIndices(aIS,&anchors);
4219:     PetscSectionGetChart(aSec,&aStart,&aEnd);
4220:     /* figure out how many points are going to be in the new element matrix
4221:      * (we allow double counting, because it's all just going to be summed
4222:      * into the global matrix anyway) */
4223:     for (p = 0; p < 2*numPoints; p+=2) {
4224:       PetscInt b    = points[p];
4225:       PetscInt bDof = 0;

4227:       if (b >= aStart && b < aEnd) {
4228:         PetscSectionGetDof(aSec,b,&bDof);
4229:       }
4230:       if (bDof) {
4231:         /* this point is constrained */
4232:         /* it is going to be replaced by its anchors */
4233:         PetscInt bOff, q;

4235:         anyConstrained = PETSC_TRUE;
4236:         newNumPoints  += bDof;
4237:         PetscSectionGetOffset(aSec,b,&bOff);
4238:         for (q = 0; q < bDof; q++) {
4239:           PetscInt a = anchors[bOff + q];
4240:           PetscInt aDof;

4242:           PetscSectionGetDof(section,a,&aDof);
4243:           newNumIndices += aDof;
4244:           for (f = 0; f < numFields; ++f) {
4245:             PetscInt fDof;

4247:             PetscSectionGetFieldDof(section, a, f, &fDof);
4248:             newOffsets[f+1] += fDof;
4249:           }
4250:         }
4251:       }
4252:       else {
4253:         /* this point is not constrained */
4254:         newNumPoints++;
4255:         PetscSectionGetDof(section,b,&bDof);
4256:         newNumIndices += bDof;
4257:         for (f = 0; f < numFields; ++f) {
4258:           PetscInt fDof;

4260:           PetscSectionGetFieldDof(section, b, f, &fDof);
4261:           newOffsets[f+1] += fDof;
4262:         }
4263:       }
4264:     }
4265:   }
4266:   if (!anyConstrained) {
4267:     *outNumPoints  = 0;
4268:     *outNumIndices = 0;
4269:     *outPoints     = NULL;
4270:     *outValues     = NULL;
4271:     if (aSec) {
4272:       ISRestoreIndices(aIS,&anchors);
4273:     }
4274:     return(0);
4275:   }

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

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

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

4283:   /* output arrays */
4284:   DMGetWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4285:   DMGetWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);

4287:   /* workspaces */
4288:   DMGetWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4289:   if (numFields) {
4290:     for (f = 0; f < numFields; f++) {
4291:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4292:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4293:     }
4294:   }
4295:   else {
4296:     DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4297:     DMGetWorkArray(dm,numPoints,PETSC_INT,&newPointOffsets[0]);
4298:   }

4300:   /* get workspaces for the point-to-point matrices */
4301:   if (numFields) {
4302:     for (p = 0; p < numPoints; p++) {
4303:       PetscInt b    = points[2*p];
4304:       PetscInt bDof = 0;

4306:       if (b >= aStart && b < aEnd) {
4307:         PetscSectionGetDof(aSec, b, &bDof);
4308:       }
4309:       if (bDof) {
4310:         for (f = 0; f < numFields; f++) {
4311:           PetscInt fDof, q, bOff, allFDof = 0;

4313:           PetscSectionGetFieldDof(section, b, f, &fDof);
4314:           PetscSectionGetOffset(aSec, b, &bOff);
4315:           for (q = 0; q < bDof; q++) {
4316:             PetscInt a = anchors[bOff + q];
4317:             PetscInt aFDof;

4319:             PetscSectionGetFieldDof(section, a, f, &aFDof);
4320:             allFDof += aFDof;
4321:           }
4322:           newPointOffsets[f][p+1] = allFDof;
4323:           pointMatOffsets[f][p+1] = fDof * allFDof;
4324:         }
4325:       }
4326:       else {
4327:         for (f = 0; f < numFields; f++) {
4328:           PetscInt fDof;

4330:           PetscSectionGetFieldDof(section, b, f, &fDof);
4331:           newPointOffsets[f][p+1] = fDof;
4332:           pointMatOffsets[f][p+1] = 0;
4333:         }
4334:       }
4335:     }
4336:     for (f = 0; f < numFields; f++) {
4337:       newPointOffsets[f][0] = 0;
4338:       pointMatOffsets[f][0] = 0;
4339:       for (p = 0; p < numPoints; p++) {
4340:         newPointOffsets[f][p+1] += newPointOffsets[f][p];
4341:         pointMatOffsets[f][p+1] += pointMatOffsets[f][p];
4342:       }
4343:       DMGetWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4344:     }
4345:   }
4346:   else {
4347:     for (p = 0; p < numPoints; p++) {
4348:       PetscInt b    = points[2*p];
4349:       PetscInt bDof = 0;

4351:       if (b >= aStart && b < aEnd) {
4352:         PetscSectionGetDof(aSec, b, &bDof);
4353:       }
4354:       if (bDof) {
4355:         PetscInt dof, bOff, q, allDof = 0;

4357:         PetscSectionGetDof(section, b, &dof);
4358:         PetscSectionGetOffset(aSec, b, &bOff);
4359:         for (q = 0; q < bDof; q++) {
4360:           PetscInt a = anchors[bOff + q], aDof;

4362:           PetscSectionGetDof(section, a, &aDof);
4363:           allDof += aDof;
4364:         }
4365:         newPointOffsets[0][p+1] = allDof;
4366:         pointMatOffsets[0][p+1] = dof * allDof;
4367:       }
4368:       else {
4369:         PetscInt dof;

4371:         PetscSectionGetDof(section, b, &dof);
4372:         newPointOffsets[0][p+1] = dof;
4373:         pointMatOffsets[0][p+1] = 0;
4374:       }
4375:     }
4376:     newPointOffsets[0][0] = 0;
4377:     pointMatOffsets[0][0] = 0;
4378:     for (p = 0; p < numPoints; p++) {
4379:       newPointOffsets[0][p+1] += newPointOffsets[0][p];
4380:       pointMatOffsets[0][p+1] += pointMatOffsets[0][p];
4381:     }
4382:     DMGetWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4383:   }

4385:   /* get the point-to-point matrices; construct newPoints */
4386:   PetscSectionGetMaxDof(aSec, &maxAnchor);
4387:   PetscSectionGetMaxDof(section, &maxDof);
4388:   DMGetWorkArray(dm,maxDof,PETSC_INT,&indices);
4389:   DMGetWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4390:   if (numFields) {
4391:     for (p = 0, newP = 0; p < numPoints; p++) {
4392:       PetscInt b    = points[2*p];
4393:       PetscInt o    = points[2*p+1];
4394:       PetscInt bDof = 0;

4396:       if (b >= aStart && b < aEnd) {
4397:         PetscSectionGetDof(aSec, b, &bDof);
4398:       }
4399:       if (bDof) {
4400:         PetscInt fStart[32], fEnd[32], fAnchorStart[32], fAnchorEnd[32], bOff, q;

4402:         fStart[0] = 0;
4403:         fEnd[0]   = 0;
4404:         for (f = 0; f < numFields; f++) {
4405:           PetscInt fDof;

4407:           PetscSectionGetFieldDof(cSec, b, f, &fDof);
4408:           fStart[f+1] = fStart[f] + fDof;
4409:           fEnd[f+1]   = fStart[f+1];
4410:         }
4411:         PetscSectionGetOffset(cSec, b, &bOff);
4412:         indicesPointFields_private(cSec, b, bOff, fEnd, PETSC_TRUE, o, indices);

4414:         fAnchorStart[0] = 0;
4415:         fAnchorEnd[0]   = 0;
4416:         for (f = 0; f < numFields; f++) {
4417:           PetscInt fDof = newPointOffsets[f][p + 1] - newPointOffsets[f][p];

4419:           fAnchorStart[f+1] = fAnchorStart[f] + fDof;
4420:           fAnchorEnd[f+1]   = fAnchorStart[f + 1];
4421:         }
4422:         PetscSectionGetOffset (aSec, b, &bOff);
4423:         for (q = 0; q < bDof; q++) {
4424:           PetscInt a = anchors[bOff + q], aOff;

4426:           /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4427:           newPoints[2*(newP + q)]     = a;
4428:           newPoints[2*(newP + q) + 1] = 0;
4429:           PetscSectionGetOffset(section, a, &aOff);
4430:           indicesPointFields_private(section, a, aOff, fAnchorEnd, PETSC_TRUE, 0, newIndices);
4431:         }
4432:         newP += bDof;

4434:         /* get the point-to-point submatrix */
4435:         for (f = 0; f < numFields; f++) {
4436:           MatGetValues(cMat,fEnd[f]-fStart[f],indices + fStart[f],fAnchorEnd[f] - fAnchorStart[f],newIndices + fAnchorStart[f],pointMat[f] + pointMatOffsets[f][p]);
4437:         }
4438:       }
4439:       else {
4440:         newPoints[2 * newP]     = b;
4441:         newPoints[2 * newP + 1] = o;
4442:         newP++;
4443:       }
4444:     }
4445:   } else {
4446:     for (p = 0; p < numPoints; p++) {
4447:       PetscInt b    = points[2*p];
4448:       PetscInt o    = points[2*p+1];
4449:       PetscInt bDof = 0;

4451:       if (b >= aStart && b < aEnd) {
4452:         PetscSectionGetDof(aSec, b, &bDof);
4453:       }
4454:       if (bDof) {
4455:         PetscInt bEnd = 0, bAnchorEnd = 0, bOff;

4457:         PetscSectionGetOffset(cSec, b, &bOff);
4458:         indicesPoint_private(cSec, b, bOff, &bEnd, PETSC_TRUE, o, indices);

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

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

4466:           newPoints[2*(newP + q)]     = a;
4467:           newPoints[2*(newP + q) + 1] = 0;
4468:           PetscSectionGetOffset(section, a, &aOff);
4469:           indicesPoint_private(section, a, aOff, &bAnchorEnd, PETSC_TRUE, 0, newIndices);
4470:         }
4471:         newP += bDof;

4473:         /* get the point-to-point submatrix */
4474:         MatGetValues(cMat,bEnd,indices,bAnchorEnd,newIndices,pointMat[0] + pointMatOffsets[0][p]);
4475:       }
4476:       else {
4477:         newPoints[2 * newP]     = b;
4478:         newPoints[2 * newP + 1] = o;
4479:         newP++;
4480:       }
4481:     }
4482:   }

4484:   PetscMemzero(tmpValues,newNumIndices*numIndices*sizeof(*tmpValues));
4485:   /* multiply constraints on the right */
4486:   if (numFields) {
4487:     for (f = 0; f < numFields; f++) {
4488:       PetscInt oldOff = offsets[f];

4490:       for (p = 0; p < numPoints; p++) {
4491:         PetscInt cStart = newPointOffsets[f][p];
4492:         PetscInt b      = points[2 * p];
4493:         PetscInt c, r, k;
4494:         PetscInt dof;

4496:         PetscSectionGetFieldDof(section,b,f,&dof);
4497:         if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4498:           PetscInt nCols         = newPointOffsets[f][p+1]-cStart;
4499:           const PetscScalar *mat = pointMat[f] + pointMatOffsets[f][p];

4501:           for (r = 0; r < numIndices; r++) {
4502:             for (c = 0; c < nCols; c++) {
4503:               for (k = 0; k < dof; k++) {
4504:                 tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4505:               }
4506:             }
4507:           }
4508:         }
4509:         else {
4510:           /* copy this column as is */
4511:           for (r = 0; r < numIndices; r++) {
4512:             for (c = 0; c < dof; c++) {
4513:               tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4514:             }
4515:           }
4516:         }
4517:         oldOff += dof;
4518:       }
4519:     }
4520:   }
4521:   else {
4522:     PetscInt oldOff = 0;
4523:     for (p = 0; p < numPoints; p++) {
4524:       PetscInt cStart = newPointOffsets[0][p];
4525:       PetscInt b      = points[2 * p];
4526:       PetscInt c, r, k;
4527:       PetscInt dof;

4529:       PetscSectionGetDof(section,b,&dof);
4530:       if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4531:         PetscInt nCols         = newPointOffsets[0][p+1]-cStart;
4532:         const PetscScalar *mat = pointMat[0] + pointMatOffsets[0][p];

4534:         for (r = 0; r < numIndices; r++) {
4535:           for (c = 0; c < nCols; c++) {
4536:             for (k = 0; k < dof; k++) {
4537:               tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4538:             }
4539:           }
4540:         }
4541:       }
4542:       else {
4543:         /* copy this column as is */
4544:         for (r = 0; r < numIndices; r++) {
4545:           for (c = 0; c < dof; c++) {
4546:             tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4547:           }
4548:         }
4549:       }
4550:       oldOff += dof;
4551:     }
4552:   }

4554:   PetscMemzero(newValues,newNumIndices*newNumIndices*sizeof(*newValues));
4555:   /* multiply constraints transpose on the left */
4556:   if (numFields) {
4557:     for (f = 0; f < numFields; f++) {
4558:       PetscInt oldOff = offsets[f];

4560:       for (p = 0; p < numPoints; p++) {
4561:         PetscInt rStart = newPointOffsets[f][p];
4562:         PetscInt b      = points[2 * p];
4563:         PetscInt c, r, k;
4564:         PetscInt dof;

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

4571:           for (r = 0; r < nRows; r++) {
4572:             for (c = 0; c < newNumIndices; c++) {
4573:               for (k = 0; k < dof; k++) {
4574:                 newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4575:               }
4576:             }
4577:           }
4578:         }
4579:         else {
4580:           /* copy this row as is */
4581:           for (r = 0; r < dof; r++) {
4582:             for (c = 0; c < newNumIndices; c++) {
4583:               newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4584:             }
4585:           }
4586:         }
4587:         oldOff += dof;
4588:       }
4589:     }
4590:   }
4591:   else {
4592:     PetscInt oldOff = 0;

4594:     for (p = 0; p < numPoints; p++) {
4595:       PetscInt rStart = newPointOffsets[0][p];
4596:       PetscInt b      = points[2 * p];
4597:       PetscInt c, r, k;
4598:       PetscInt dof;

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

4605:         for (r = 0; r < nRows; r++) {
4606:           for (c = 0; c < newNumIndices; c++) {
4607:             for (k = 0; k < dof; k++) {
4608:               newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4609:             }
4610:           }
4611:         }
4612:       }
4613:       else {
4614:         /* copy this row as is */
4615:         for (r = 0; r < dof; c++) {
4616:           for (c = 0; c < newNumIndices; c++) {
4617:             newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4618:           }
4619:         }
4620:       }
4621:       oldOff += dof;
4622:     }
4623:   }

4625:   /* clean up */
4626:   DMRestoreWorkArray(dm,maxDof,PETSC_INT,&indices);
4627:   DMRestoreWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4628:   if (numFields) {
4629:     for (f = 0; f < numFields; f++) {
4630:       DMRestoreWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4631:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4632:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4633:     }
4634:   }
4635:   else {
4636:     DMRestoreWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4637:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4638:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[0]);
4639:   }
4640:   DMRestoreWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4641:   ISRestoreIndices(aIS,&anchors);

4643:   /* output */
4644:   *outNumPoints  = newNumPoints;
4645:   *outNumIndices = newNumIndices;
4646:   *outPoints     = newPoints;
4647:   *outValues     = newValues;
4648:   for (f = 0; f < numFields; f++) {
4649:     offsets[f] = newOffsets[f];
4650:   }
4651:   return(0);
4652: }

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

4659:   Not collective

4661:   Input Parameters:
4662: + dm - The DM
4663: . section - The section describing the layout in v, or NULL to use the default section
4664: . globalSection - The section describing the layout in v, or NULL to use the default global section
4665: . A - The matrix
4666: . point - The sieve point in the DM
4667: . values - The array of values
4668: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

4673:   Level: intermediate

4675: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
4676: @*/
4677: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4678: {
4679:   DM_Plex        *mesh   = (DM_Plex*) dm->data;
4680:   PetscSection    clSection;
4681:   IS              clPoints;
4682:   PetscInt       *points = NULL, *newPoints;
4683:   const PetscInt *clp;
4684:   PetscInt       *indices;
4685:   PetscInt        offsets[32];
4686:   PetscInt        numFields, numPoints, newNumPoints, numIndices, newNumIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
4687:   PetscScalar    *newValues;
4688:   PetscErrorCode  ierr;

4692:   if (!section) {DMGetDefaultSection(dm, &section);}
4694:   if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
4697:   PetscSectionGetNumFields(section, &numFields);
4698:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4699:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
4700:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4701:   if (!clPoints) {
4702:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4703:     /* Compress out points not in the section */
4704:     PetscSectionGetChart(section, &pStart, &pEnd);
4705:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4706:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4707:         points[q*2]   = points[p];
4708:         points[q*2+1] = points[p+1];
4709:         ++q;
4710:       }
4711:     }
4712:     numPoints = q;
4713:   } else {
4714:     PetscInt dof, off;

4716:     PetscSectionGetDof(clSection, point, &dof);
4717:     numPoints = dof/2;
4718:     PetscSectionGetOffset(clSection, point, &off);
4719:     ISGetIndices(clPoints, &clp);
4720:     points = (PetscInt *) &clp[off];
4721:   }
4722:   for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
4723:     PetscInt fdof;

4725:     PetscSectionGetDof(section, points[p], &dof);
4726:     for (f = 0; f < numFields; ++f) {
4727:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
4728:       offsets[f+1] += fdof;
4729:     }
4730:     numIndices += dof;
4731:   }
4732:   for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];

4734:   if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
4735:   DMPlexAnchorsModifyMat(dm,section,numPoints,numIndices,points,values,&newNumPoints,&newNumIndices,&newPoints,&newValues,offsets);
4736:   if (newNumPoints) {
4737:     if (!clPoints) {
4738:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4739:     } else {
4740:       ISRestoreIndices(clPoints, &clp);
4741:     }
4742:     numPoints  = newNumPoints;
4743:     numIndices = newNumIndices;
4744:     points     = newPoints;
4745:     values     = newValues;
4746:   }
4747:   DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
4748:   if (numFields) {
4749:     for (p = 0; p < numPoints*2; p += 2) {
4750:       PetscInt o = points[p+1];
4751:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4752:       indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
4753:     }
4754:   } else {
4755:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
4756:       PetscInt o = points[p+1];
4757:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4758:       indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
4759:     }
4760:   }
4761:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
4762:   MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
4763:   if (mesh->printFEM > 1) {
4764:     PetscInt i;
4765:     PetscPrintf(PETSC_COMM_SELF, "  Indices:");
4766:     for (i = 0; i < numIndices; ++i) {PetscPrintf(PETSC_COMM_SELF, " %d", indices[i]);}
4767:     PetscPrintf(PETSC_COMM_SELF, "\n");
4768:   }
4769:   if (ierr) {
4770:     PetscMPIInt    rank;
4771:     PetscErrorCode ierr2;

4773:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4774:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4775:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
4776:     ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
4777: 
4778:   }
4779:   if (newNumPoints) {
4780:     DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4781:     DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4782:   }
4783:   else {
4784:     if (!clPoints) {
4785:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4786:     } else {
4787:       ISRestoreIndices(clPoints, &clp);
4788:     }
4789:   }
4790:   DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
4791:   return(0);
4792: }

4796: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4797: {
4798:   DM_Plex        *mesh   = (DM_Plex*) dmf->data;
4799:   PetscInt       *fpoints = NULL, *ftotpoints = NULL;
4800:   PetscInt       *cpoints = NULL;
4801:   PetscInt       *findices, *cindices;
4802:   PetscInt        foffsets[32], coffsets[32];
4803:   CellRefiner     cellRefiner;
4804:   PetscInt        numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4805:   PetscErrorCode  ierr;

4810:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4812:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
4814:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4816:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4819:   PetscSectionGetNumFields(fsection, &numFields);
4820:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4821:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4822:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4823:   /* Column indices */
4824:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4825:   maxFPoints = numCPoints;
4826:   /* Compress out points not in the section */
4827:   /*   TODO: Squeeze out points with 0 dof as well */
4828:   PetscSectionGetChart(csection, &pStart, &pEnd);
4829:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4830:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4831:       cpoints[q*2]   = cpoints[p];
4832:       cpoints[q*2+1] = cpoints[p+1];
4833:       ++q;
4834:     }
4835:   }
4836:   numCPoints = q;
4837:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4838:     PetscInt fdof;

4840:     PetscSectionGetDof(csection, cpoints[p], &dof);
4841:     if (!dof) continue;
4842:     for (f = 0; f < numFields; ++f) {
4843:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4844:       coffsets[f+1] += fdof;
4845:     }
4846:     numCIndices += dof;
4847:   }
4848:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4849:   /* Row indices */
4850:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4851:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4852:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4853:   for (r = 0, q = 0; r < numSubcells; ++r) {
4854:     /* TODO Map from coarse to fine cells */
4855:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4856:     /* Compress out points not in the section */
4857:     PetscSectionGetChart(fsection, &pStart, &pEnd);
4858:     for (p = 0; p < numFPoints*2; p += 2) {
4859:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4860:         PetscSectionGetDof(fsection, fpoints[p], &dof);
4861:         if (!dof) continue;
4862:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4863:         if (s < q) continue;
4864:         ftotpoints[q*2]   = fpoints[p];
4865:         ftotpoints[q*2+1] = fpoints[p+1];
4866:         ++q;
4867:       }
4868:     }
4869:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4870:   }
4871:   numFPoints = q;
4872:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4873:     PetscInt fdof;

4875:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4876:     if (!dof) continue;
4877:     for (f = 0; f < numFields; ++f) {
4878:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4879:       foffsets[f+1] += fdof;
4880:     }
4881:     numFIndices += dof;
4882:   }
4883:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

4885:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4886:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4887:   DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4888:   DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4889:   if (numFields) {
4890:     for (p = 0; p < numFPoints*2; p += 2) {
4891:       PetscInt o = ftotpoints[p+1];
4892:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4893:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4894:     }
4895:     for (p = 0; p < numCPoints*2; p += 2) {
4896:       PetscInt o = cpoints[p+1];
4897:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4898:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4899:     }
4900:   } else {
4901:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4902:       PetscInt o = ftotpoints[p+1];
4903:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4904:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4905:     }
4906:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4907:       PetscInt o = cpoints[p+1];
4908:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4909:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4910:     }
4911:   }
4912:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
4913:   MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
4914:   if (ierr) {
4915:     PetscMPIInt    rank;
4916:     PetscErrorCode ierr2;

4918:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4919:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4920:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
4921:     ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
4922:     ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
4923: 
4924:   }
4925:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4926:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4927:   DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4928:   DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4929:   return(0);
4930: }

4934: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
4935: {
4936:   PetscInt      *fpoints = NULL, *ftotpoints = NULL;
4937:   PetscInt      *cpoints = NULL;
4938:   PetscInt       foffsets[32], coffsets[32];
4939:   CellRefiner    cellRefiner;
4940:   PetscInt       numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;

4946:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4948:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
4950:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4952:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4954:   PetscSectionGetNumFields(fsection, &numFields);
4955:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4956:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4957:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4958:   /* Column indices */
4959:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4960:   maxFPoints = numCPoints;
4961:   /* Compress out points not in the section */
4962:   /*   TODO: Squeeze out points with 0 dof as well */
4963:   PetscSectionGetChart(csection, &pStart, &pEnd);
4964:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4965:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4966:       cpoints[q*2]   = cpoints[p];
4967:       cpoints[q*2+1] = cpoints[p+1];
4968:       ++q;
4969:     }
4970:   }
4971:   numCPoints = q;
4972:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4973:     PetscInt fdof;

4975:     PetscSectionGetDof(csection, cpoints[p], &dof);
4976:     if (!dof) continue;
4977:     for (f = 0; f < numFields; ++f) {
4978:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4979:       coffsets[f+1] += fdof;
4980:     }
4981:     numCIndices += dof;
4982:   }
4983:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4984:   /* Row indices */
4985:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4986:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4987:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4988:   for (r = 0, q = 0; r < numSubcells; ++r) {
4989:     /* TODO Map from coarse to fine cells */
4990:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4991:     /* Compress out points not in the section */
4992:     PetscSectionGetChart(fsection, &pStart, &pEnd);
4993:     for (p = 0; p < numFPoints*2; p += 2) {
4994:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4995:         PetscSectionGetDof(fsection, fpoints[p], &dof);
4996:         if (!dof) continue;
4997:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4998:         if (s < q) continue;
4999:         ftotpoints[q*2]   = fpoints[p];
5000:         ftotpoints[q*2+1] = fpoints[p+1];
5001:         ++q;
5002:       }
5003:     }
5004:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
5005:   }
5006:   numFPoints = q;
5007:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
5008:     PetscInt fdof;

5010:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
5011:     if (!dof) continue;
5012:     for (f = 0; f < numFields; ++f) {
5013:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
5014:       foffsets[f+1] += fdof;
5015:     }
5016:     numFIndices += dof;
5017:   }
5018:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

5020:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
5021:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
5022:   if (numFields) {
5023:     for (p = 0; p < numFPoints*2; p += 2) {
5024:       PetscInt o = ftotpoints[p+1];
5025:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5026:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
5027:     }
5028:     for (p = 0; p < numCPoints*2; p += 2) {
5029:       PetscInt o = cpoints[p+1];
5030:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5031:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5032:     }
5033:   } else {
5034:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5035:       PetscInt o = ftotpoints[p+1];
5036:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5037:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5038:     }
5039:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5040:       PetscInt o = cpoints[p+1];
5041:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5042:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5043:     }
5044:   }
5045:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5046:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5047:   return(0);
5048: }

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

5055:   Input Parameter:
5056: . dm - The DMPlex object

5058:   Output Parameters:
5059: + cMax - The first hybrid cell
5060: . fMax - The first hybrid face
5061: . eMax - The first hybrid edge
5062: - vMax - The first hybrid vertex

5064:   Level: developer

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

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

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

5089:   Input Parameters:
5090: . dm   - The DMPlex object
5091: . cMax - The first hybrid cell
5092: . fMax - The first hybrid face
5093: . eMax - The first hybrid edge
5094: - vMax - The first hybrid vertex

5096:   Level: developer

5098: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5099: @*/
5100: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5101: {
5102:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5103:   PetscInt       dim;

5108:   DMGetDimension(dm, &dim);
5109:   if (cMax >= 0) mesh->hybridPointMax[dim]   = cMax;
5110:   if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5111:   if (eMax >= 0) mesh->hybridPointMax[1]     = eMax;
5112:   if (vMax >= 0) mesh->hybridPointMax[0]     = vMax;
5113:   return(0);
5114: }

5118: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5119: {
5120:   DM_Plex *mesh = (DM_Plex*) dm->data;

5125:   *cellHeight = mesh->vtkCellHeight;
5126:   return(0);
5127: }

5131: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5132: {
5133:   DM_Plex *mesh = (DM_Plex*) dm->data;

5137:   mesh->vtkCellHeight = cellHeight;
5138:   return(0);
5139: }

5143: /* We can easily have a form that takes an IS instead */
5144: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5145: {
5146:   PetscSection   section, globalSection;
5147:   PetscInt      *numbers, p;

5151:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
5152:   PetscSectionSetChart(section, pStart, pEnd);
5153:   for (p = pStart; p < pEnd; ++p) {
5154:     PetscSectionSetDof(section, p, 1);
5155:   }
5156:   PetscSectionSetUp(section);
5157:   PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
5158:   PetscMalloc1(pEnd - pStart, &numbers);
5159:   for (p = pStart; p < pEnd; ++p) {
5160:     PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5161:     if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5162:     else                       numbers[p-pStart] += shift;
5163:   }
5164:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5165:   if (globalSize) {
5166:     PetscLayout layout;
5167:     PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5168:     PetscLayoutGetSize(layout, globalSize);
5169:     PetscLayoutDestroy(&layout);
5170:   }
5171:   PetscSectionDestroy(&section);
5172:   PetscSectionDestroy(&globalSection);
5173:   return(0);
5174: }

5178: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5179: {
5180:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5181:   PetscInt       cellHeight, cStart, cEnd, cMax;

5186:   if (!mesh->globalCellNumbers) {
5187:     DMPlexGetVTKCellHeight(dm, &cellHeight);
5188:     DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5189:     DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5190:     if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
5191:     DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, &mesh->globalCellNumbers);
5192:   }
5193:   *globalCellNumbers = mesh->globalCellNumbers;
5194:   return(0);
5195: }

5199: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5200: {
5201:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5202:   PetscInt       vStart, vEnd, vMax;

5207:   if (!mesh->globalVertexNumbers) {
5208:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5209:     DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5210:     if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
5211:     DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, &mesh->globalVertexNumbers);
5212:   }
5213:   *globalVertexNumbers = mesh->globalVertexNumbers;
5214:   return(0);
5215: }

5219: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5220: {
5221:   IS             nums[4];
5222:   PetscInt       depths[4];
5223:   PetscInt       depth, d, shift = 0;

5228:   DMPlexGetDepth(dm, &depth);
5229:   /* For unstratified meshes use dim instead of depth */
5230:   if (depth < 0) {DMGetDimension(dm, &depth);}
5231:   depths[0] = depth; depths[1] = 0;
5232:   for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5233:   for (d = 0; d <= depth; ++d) {
5234:     PetscInt pStart, pEnd, gsize;

5236:     DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5237:     DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5238:     shift += gsize;
5239:   }
5240:   ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5241:   for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5242:   return(0);
5243: }


5248: /*@C
5249:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
5250:   the local section and an SF describing the section point overlap.

5252:   Input Parameters:
5253:   + s - The PetscSection for the local field layout
5254:   . sf - The SF describing parallel layout of the section points
5255:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
5256:   . label - The label specifying the points
5257:   - labelValue - The label stratum specifying the points

5259:   Output Parameter:
5260:   . gsection - The PetscSection for the global field layout

5262:   Note: This gives negative sizes and offsets to points not owned by this process

5264:   Level: developer

5266: .seealso: PetscSectionCreate()
5267: @*/
5268: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
5269: {
5270:   PetscInt      *neg = NULL, *tmpOff = NULL;
5271:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

5275:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
5276:   PetscSectionGetChart(s, &pStart, &pEnd);
5277:   PetscSectionSetChart(*gsection, pStart, pEnd);
5278:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
5279:   if (nroots >= 0) {
5280:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
5281:     PetscCalloc1(nroots, &neg);
5282:     if (nroots > pEnd-pStart) {
5283:       PetscCalloc1(nroots, &tmpOff);
5284:     } else {
5285:       tmpOff = &(*gsection)->atlasDof[-pStart];
5286:     }
5287:   }
5288:   /* Mark ghost points with negative dof */
5289:   for (p = pStart; p < pEnd; ++p) {
5290:     PetscInt value;

5292:     DMLabelGetValue(label, p, &value);
5293:     if (value != labelValue) continue;
5294:     PetscSectionGetDof(s, p, &dof);
5295:     PetscSectionSetDof(*gsection, p, dof);
5296:     PetscSectionGetConstraintDof(s, p, &cdof);
5297:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
5298:     if (neg) neg[p] = -(dof+1);
5299:   }
5300:   PetscSectionSetUpBC(*gsection);
5301:   if (nroots >= 0) {
5302:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5303:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5304:     if (nroots > pEnd-pStart) {
5305:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
5306:     }
5307:   }
5308:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
5309:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5310:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
5311:     (*gsection)->atlasOff[p] = off;
5312:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
5313:   }
5314:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
5315:   globalOff -= off;
5316:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5317:     (*gsection)->atlasOff[p] += globalOff;
5318:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
5319:   }
5320:   /* Put in negative offsets for ghost points */
5321:   if (nroots >= 0) {
5322:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5323:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5324:     if (nroots > pEnd-pStart) {
5325:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
5326:     }
5327:   }
5328:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
5329:   PetscFree(neg);
5330:   return(0);
5331: }

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

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

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

5343:   Level: developer

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

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

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

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

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

5431:   Level: developer

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

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

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

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

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

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

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

5491:   Level: developer

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

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

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

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

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

5562:   /*
5563:   Loop over coarse cells
5564:     Loop over coarse basis functions
5565:       Loop over fine cells in coarse cell
5566:         Loop over fine dual basis functions
5567:           Evaluate coarse basis on fine dual basis quad points
5568:           Sum
5569:           Update local element matrix
5570:     Accumulate to interpolation matrix

5572:    Can extend PetscFEIntegrateJacobian_Basic() to do a specialized cell loop
5573:   */
5574:   DMGetDefaultGlobalSection(dmFine, &gsf);
5575:   PetscSectionGetConstrainedStorageSize(gsf, &m);
5576:   DMGetDefaultGlobalSection(dmCoarse, &gsc);
5577:   PetscSectionGetConstrainedStorageSize(gsc, &n);
5578:   /* We need to preallocate properly */
5579:   MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
5580:   MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
5581:   MatSetType(*interpolation, dmCoarse->mattype);
5582:   DMGetApplicationContext(dmFine, &ctx);
5583:   DMPlexComputeInterpolatorFEM(dmCoarse, dmFine, *interpolation, ctx);
5584:   /* Use naive scaling */
5585:   DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
5586:   return(0);
5587: }

5591: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat)
5592: {
5594:   VecScatter     ctx;

5597:   DMPlexComputeInjectorFEM(dmCoarse, dmFine, &ctx, NULL);
5598:   MatCreateScatter(PetscObjectComm((PetscObject)ctx), ctx, mat);
5599:   VecScatterDestroy(&ctx);
5600:   return(0);
5601: }

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

5616:   DMGetNumFields(dm, &numFields);
5617:   /* FE and FV boundary conditions are handled slightly differently */
5618:   PetscMalloc1(numFields, &isFE);
5619:   for (f = 0; f < numFields; ++f) {
5620:     PetscObject  obj;
5621:     PetscClassId id;

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

5639:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL, NULL);
5640:     if (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:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
5665:     if (!isFE[field]) continue;
5666:     DMPlexGetLabel(dm, bdLabel, &label);
5667:     /* Only want to modify label once */
5668:     for (bd2 = 0; bd2 < bd; ++bd2) {
5669:       const char *bdname;
5670:       DMPlexGetBoundary(dm, 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:       DMPlexLabelComplete(dm, label);
5676:       DMPlexLabelAddCells(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:         DMPlexGetStratumIS(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:         DMPlexGetStratumIS(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:   DMPlexGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

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

5775:   Output Parameter:
5776: . cdm - The coarse DM

5778:   Level: intermediate

5780: .seealso: DMPlexSetCoarseDM()
5781: @*/
5782: PetscErrorCode DMPlexGetCoarseDM(DM dm, DM *cdm)
5783: {
5787:   *cdm = ((DM_Plex *) dm->data)->coarseMesh;
5788:   return(0);
5789: }

5793: /*@
5794:   DMPlexSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

5796:   Input Parameters:
5797: + dm - The DMPlex object
5798: - cdm - The coarse DM

5800:   Level: intermediate

5802: .seealso: DMPlexGetCoarseDM()
5803: @*/
5804: PetscErrorCode DMPlexSetCoarseDM(DM dm, DM cdm)
5805: {
5806:   DM_Plex       *mesh;

5812:   mesh = (DM_Plex *) dm->data;
5813:   DMDestroy(&mesh->coarseMesh);
5814:   mesh->coarseMesh = cdm;
5815:   PetscObjectReference((PetscObject) mesh->coarseMesh);
5816:   return(0);
5817: }

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

5826:   not collective

5828:   Input Parameters:
5829: . dm - The DMPlex object

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


5836:   Level: intermediate

5838: .seealso: DMPlexSetAnchors(), DMGetConstraints(), DMSetConstraints()
5839: @*/
5840: PetscErrorCode DMPlexGetAnchors(DM dm, PetscSection *anchorSection, IS *anchorIS)
5841: {
5842:   DM_Plex *plex = (DM_Plex *)dm->data;

5847:   if (!plex->anchorSection && !plex->anchorIS && plex->createanchors) {(*plex->createanchors)(dm);}
5848:   if (anchorSection) *anchorSection = plex->anchorSection;
5849:   if (anchorIS) *anchorIS = plex->anchorIS;
5850:   return(0);
5851: }

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

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

5863:   collective on dm

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

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

5872:   Level: intermediate

5874: .seealso: DMPlexGetAnchors(), DMGetConstraints(), DMSetConstraints()
5875: @*/
5876: PetscErrorCode DMPlexSetAnchors(DM dm, PetscSection anchorSection, IS anchorIS)
5877: {
5878:   DM_Plex        *plex = (DM_Plex *)dm->data;
5879:   PetscMPIInt    result;

5884:   if (anchorSection) {
5886:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorSection),&result);
5887:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor section must have local communicator");
5888:   }
5889:   if (anchorIS) {
5891:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorIS),&result);
5892:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor IS must have local communicator");
5893:   }

5895:   PetscObjectReference((PetscObject)anchorSection);
5896:   PetscSectionDestroy(&plex->anchorSection);
5897:   plex->anchorSection = anchorSection;

5899:   PetscObjectReference((PetscObject)anchorIS);
5900:   ISDestroy(&plex->anchorIS);
5901:   plex->anchorIS = anchorIS;

5903: #if defined(PETSC_USE_DEBUG)
5904:   if (anchorIS && anchorSection) {
5905:     PetscInt size, a, pStart, pEnd;
5906:     const PetscInt *anchors;

5908:     PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5909:     ISGetLocalSize(anchorIS,&size);
5910:     ISGetIndices(anchorIS,&anchors);
5911:     for (a = 0; a < size; a++) {
5912:       PetscInt p;

5914:       p = anchors[a];
5915:       if (p >= pStart && p < pEnd) {
5916:         PetscInt dof;

5918:         PetscSectionGetDof(anchorSection,p,&dof);
5919:         if (dof) {
5920:           PetscErrorCode ierr2;

5922:           ierr2 = ISRestoreIndices(anchorIS,&anchors);CHKERRQ(ierr2);
5923:           SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point %d cannot be constrained and an anchor",p);
5924:         }
5925:       }
5926:     }
5927:     ISRestoreIndices(anchorIS,&anchors);
5928:   }
5929: #endif
5930:   /* reset the generic constraints */
5931:   DMSetDefaultConstraints(dm,NULL,NULL);
5932:   return(0);
5933: }

5937: static PetscErrorCode DMPlexCreateConstraintSection_Anchors(DM dm, PetscSection section, PetscSection *cSec)
5938: {
5939:   PetscSection anchorSection;
5940:   PetscInt pStart, pEnd, p, dof, numFields, f;

5945:   DMPlexGetAnchors(dm,&anchorSection,NULL);
5946:   PetscSectionCreate(PETSC_COMM_SELF,cSec);
5947:   PetscSectionGetNumFields(section,&numFields);
5948:   PetscSectionSetNumFields(*cSec,numFields);
5949:   PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5950:   PetscSectionSetChart(*cSec,pStart,pEnd);
5951:   for (p = pStart; p < pEnd; p++) {
5952:     PetscSectionGetDof(anchorSection,p,&dof);
5953:     if (dof) {
5954:       PetscSectionGetDof(section,p,&dof);
5955:       PetscSectionSetDof(*cSec,p,dof);
5956:       for (f = 0; f < numFields; f++) {
5957:         PetscSectionGetFieldDof(section,p,f,&dof);
5958:         PetscSectionSetFieldDof(*cSec,p,f,dof);
5959:       }
5960:     }
5961:   }
5962:   PetscSectionSetUp(*cSec);
5963:   return(0);
5964: }

5968: static PetscErrorCode DMPlexCreateConstraintMatrix_Anchors(DM dm, PetscSection section, PetscSection cSec, Mat *cMat)
5969: {
5970:   PetscSection aSec;
5971:   PetscInt pStart, pEnd, p, dof, aDof, aOff, off, nnz, annz, m, n, q, a, offset, *i, *j;
5972:   const PetscInt *anchors;
5973:   PetscInt numFields, f;
5974:   IS aIS;

5979:   PetscSectionGetStorageSize(cSec, &m);
5980:   PetscSectionGetStorageSize(section, &n);
5981:   MatCreate(PETSC_COMM_SELF,cMat);
5982:   MatSetSizes(*cMat,m,n,m,n);
5983:   MatSetType(*cMat,MATSEQAIJ);
5984:   DMPlexGetAnchors(dm,&aSec,&aIS);
5985:   ISGetIndices(aIS,&anchors);
5986:   PetscSectionGetChart(aSec,&pStart,&pEnd);
5987:   PetscMalloc1(m+1,&i);
5988:   i[0] = 0;
5989:   PetscSectionGetNumFields(section,&numFields);
5990:   for (p = pStart; p < pEnd; p++) {
5991:     PetscSectionGetDof(aSec,p,&dof);
5992:     if (!dof) continue;
5993:     PetscSectionGetOffset(aSec,p,&off);
5994:     if (numFields) {
5995:       for (f = 0; f < numFields; f++) {
5996:         annz = 0;
5997:         for (q = 0; q < dof; q++) {
5998:           a = anchors[off + q];
5999:           PetscSectionGetFieldDof(section,a,f,&aDof);
6000:           annz += aDof;
6001:         }
6002:         PetscSectionGetFieldDof(cSec,p,f,&dof);
6003:         PetscSectionGetFieldOffset(cSec,p,f,&off);
6004:         for (q = 0; q < dof; q++) {
6005:           i[off + q + 1] = i[off + q] + annz;
6006:         }
6007:       }
6008:     }
6009:     else {
6010:       annz = 0;
6011:       for (q = 0; q < dof; q++) {
6012:         a = anchors[off + q];
6013:         PetscSectionGetDof(section,a,&aDof);
6014:         annz += aDof;
6015:       }
6016:       PetscSectionGetDof(cSec,p,&dof);
6017:       PetscSectionGetOffset(cSec,p,&off);
6018:       for (q = 0; q < dof; q++) {
6019:         i[off + q + 1] = i[off + q] + annz;
6020:       }
6021:     }
6022:   }
6023:   nnz = i[m];
6024:   PetscMalloc1(nnz,&j);
6025:   offset = 0;
6026:   for (p = pStart; p < pEnd; p++) {
6027:     if (numFields) {
6028:       for (f = 0; f < numFields; f++) {
6029:         PetscSectionGetFieldDof(cSec,p,f,&dof);
6030:         for (q = 0; q < dof; q++) {
6031:           PetscInt rDof, rOff, r;
6032:           PetscSectionGetDof(aSec,p,&rDof);
6033:           PetscSectionGetOffset(aSec,p,&rOff);
6034:           for (r = 0; r < rDof; r++) {
6035:             PetscInt s;

6037:             a = anchors[rOff + r];

6039:             PetscSectionGetFieldDof(section,a,f,&aDof);
6040:             PetscSectionGetFieldOffset(section,a,f,&aOff);
6041:             for (s = 0; s < aDof; s++) {
6042:               j[offset++] = aOff + s;
6043:             }
6044:           }
6045:         }
6046:       }
6047:     }
6048:     else {
6049:       PetscSectionGetDof(cSec,p,&dof);
6050:       for (q = 0; q < dof; q++) {
6051:         PetscInt rDof, rOff, r;
6052:         PetscSectionGetDof(aSec,p,&rDof);
6053:         PetscSectionGetOffset(aSec,p,&rOff);
6054:         for (r = 0; r < rDof; r++) {
6055:           PetscInt s;

6057:           a = anchors[rOff + r];

6059:           PetscSectionGetDof(section,a,&aDof);
6060:           PetscSectionGetOffset(section,a,&aOff);
6061:           for (s = 0; s < aDof; s++) {
6062:             j[offset++] = aOff + s;
6063:           }
6064:         }
6065:       }
6066:     }
6067:   }
6068:   MatSeqAIJSetPreallocationCSR(*cMat,i,j,NULL);
6069:   PetscFree(i);
6070:   PetscFree(j);
6071:   ISRestoreIndices(aIS,&anchors);
6072:   return(0);
6073: }

6077: PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm)
6078: {
6079:   DM_Plex        *plex = (DM_Plex *)dm->data;
6080:   PetscSection   anchorSection, section, cSec;
6081:   Mat            cMat;

6086:   DMPlexGetAnchors(dm,&anchorSection,NULL);
6087:   if (anchorSection) {
6088:     PetscDS  ds;
6089:     PetscInt nf;

6091:     DMGetDefaultSection(dm,&section);
6092:     DMPlexCreateConstraintSection_Anchors(dm,section,&cSec);
6093:     DMPlexCreateConstraintMatrix_Anchors(dm,section,cSec,&cMat);
6094:     DMGetDS(dm,&ds);
6095:     PetscDSGetNumFields(ds,&nf);
6096:     if (nf && plex->computeanchormatrix) {(*plex->computeanchormatrix)(dm,section,cSec,cMat);}
6097:     DMSetDefaultConstraints(dm,cSec,cMat);
6098:     PetscSectionDestroy(&cSec);
6099:     MatDestroy(&cMat);
6100:   }
6101:   return(0);
6102: }