Actual source code: plex.c

petsc-master 2015-08-27
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:     PetscViewerASCIIPushSynchronized(viewer);
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:     PetscViewerASCIIPushSynchronized(viewer);
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", (double)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:     DMPlexGetCoarseDM(dm, &cdm);
653:     if (cdm) {
654:       PetscViewerASCIIPushTab(viewer);
655:       DMPlexView_Ascii(cdm, viewer);
656:       PetscViewerASCIIPopTab(viewer);
657:     }
658:   }
659:   return(0);
660: }

664: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
665: {
666:   PetscBool      iascii, ishdf5, isvtk;

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

694: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
695: {
696:   PetscBool      isbinary, ishdf5;

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

717: static PetscErrorCode BoundaryDestroy(DMBoundary *boundary)
718: {
719:   DMBoundary     b, next;

723:   if (!boundary) return(0);
724:   b = *boundary;
725:   *boundary = NULL;
726:   for (; b; b = next) {
727:     next = b->next;
728:     PetscFree(b->comps);
729:     PetscFree(b->ids);
730:     PetscFree(b->name);
731:     PetscFree(b->labelname);
732:     PetscFree(b);
733:   }
734:   return(0);
735: }

739: PetscErrorCode DMDestroy_Plex(DM dm)
740: {
741:   DM_Plex       *mesh = (DM_Plex*) dm->data;
742:   PlexLabel      next = mesh->labels;

746:   if (--mesh->refct > 0) return(0);
747:   PetscSectionDestroy(&mesh->coneSection);
748:   PetscFree(mesh->cones);
749:   PetscFree(mesh->coneOrientations);
750:   PetscSectionDestroy(&mesh->supportSection);
751:   PetscFree(mesh->supports);
752:   PetscFree(mesh->facesTmp);
753:   PetscFree(mesh->tetgenOpts);
754:   PetscFree(mesh->triangleOpts);
755:   PetscPartitionerDestroy(&mesh->partitioner);
756:   while (next) {
757:     PlexLabel tmp = next->next;

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

784: PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
785: {
786:   PetscSection   sectionGlobal;
787:   PetscInt       bs = -1;
788:   PetscInt       localSize;
789:   PetscBool      isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock;
791:   MatType        mtype;
792:   ISLocalToGlobalMapping ltog;

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

815:     if (bs < 0) {
816:       if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
817:         PetscInt pStart, pEnd, p, dof, cdof;

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

851:     /* Set localtoglobalmapping on the matrix for MatSetValuesLocal() to work */
852:     DMGetLocalToGlobalMapping(dm,&ltog);
853:     MatSetLocalToGlobalMapping(*J,ltog,ltog);
854:   }
855:   return(0);
856: }

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

863:   Not collective

865:   Input Parameter:
866: . mesh - The DMPlex

868:   Output Parameters:
869: + pStart - The first mesh point
870: - pEnd   - The upper bound for mesh points

872:   Level: beginner

874: .seealso: DMPlexCreate(), DMPlexSetChart()
875: @*/
876: PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
877: {
878:   DM_Plex       *mesh = (DM_Plex*) dm->data;

883:   PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
884:   return(0);
885: }

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

892:   Not collective

894:   Input Parameters:
895: + mesh - The DMPlex
896: . pStart - The first mesh point
897: - pEnd   - The upper bound for mesh points

899:   Output Parameters:

901:   Level: beginner

903: .seealso: DMPlexCreate(), DMPlexGetChart()
904: @*/
905: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
906: {
907:   DM_Plex       *mesh = (DM_Plex*) dm->data;

912:   PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
913:   PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
914:   return(0);
915: }

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

922:   Not collective

924:   Input Parameters:
925: + mesh - The DMPlex
926: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

928:   Output Parameter:
929: . size - The cone size for point p

931:   Level: beginner

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

943:   PetscSectionGetDof(mesh->coneSection, p, size);
944:   return(0);
945: }

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

952:   Not collective

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

959:   Output Parameter:

961:   Note:
962:   This should be called after DMPlexSetChart().

964:   Level: beginner

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

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

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

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

986:   Not collective

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

993:   Output Parameter:

995:   Note:
996:   This should be called after DMPlexSetChart().

998:   Level: beginner

1000: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart()
1001: @*/
1002: PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size)
1003: {
1004:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1005:   PetscInt       csize;

1010:   PetscSectionAddDof(mesh->coneSection, p, size);
1011:   PetscSectionGetDof(mesh->coneSection, p, &csize);

1013:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
1014:   return(0);
1015: }

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

1022:   Not collective

1024:   Input Parameters:
1025: + mesh - The DMPlex
1026: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1031:   Level: beginner

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

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

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

1050:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1051:   *cone = &mesh->cones[off];
1052:   return(0);
1053: }

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

1060:   Not collective

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

1067:   Output Parameter:

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

1072:   Level: beginner

1074: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1075: @*/
1076: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
1077: {
1078:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1079:   PetscInt       pStart, pEnd;
1080:   PetscInt       dof, off, c;

1085:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1086:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1088:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1089:   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);
1090:   for (c = 0; c < dof; ++c) {
1091:     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);
1092:     mesh->cones[off+c] = cone[c];
1093:   }
1094:   return(0);
1095: }

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

1102:   Not collective

1104:   Input Parameters:
1105: + mesh - The DMPlex
1106: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1114:   Level: beginner

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

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

1122: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
1123: @*/
1124: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
1125: {
1126:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1127:   PetscInt       off;

1132: #if defined(PETSC_USE_DEBUG)
1133:   {
1134:     PetscInt dof;
1135:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1137:   }
1138: #endif
1139:   PetscSectionGetOffset(mesh->coneSection, p, &off);

1141:   *coneOrientation = &mesh->coneOrientations[off];
1142:   return(0);
1143: }

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

1150:   Not collective

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

1160:   Output Parameter:

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

1165:   Level: beginner

1167: .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1168: @*/
1169: PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
1170: {
1171:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1172:   PetscInt       pStart, pEnd;
1173:   PetscInt       dof, off, c;

1178:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1179:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1181:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1182:   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);
1183:   for (c = 0; c < dof; ++c) {
1184:     PetscInt cdof, o = coneOrientation[c];

1186:     PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1187:     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);
1188:     mesh->coneOrientations[off+c] = o;
1189:   }
1190:   return(0);
1191: }

1195: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
1196: {
1197:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1198:   PetscInt       pStart, pEnd;
1199:   PetscInt       dof, off;

1204:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1205:   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);
1206:   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);
1207:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1208:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1209:   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);
1210:   mesh->cones[off+conePos] = conePoint;
1211:   return(0);
1212: }

1216: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1217: {
1218:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1219:   PetscInt       pStart, pEnd;
1220:   PetscInt       dof, off;

1225:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1226:   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);
1227:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1228:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1229:   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);
1230:   mesh->coneOrientations[off+conePos] = coneOrientation;
1231:   return(0);
1232: }

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

1239:   Not collective

1241:   Input Parameters:
1242: + mesh - The DMPlex
1243: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1245:   Output Parameter:
1246: . size - The support size for point p

1248:   Level: beginner

1250: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1251: @*/
1252: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1253: {
1254:   DM_Plex       *mesh = (DM_Plex*) dm->data;

1260:   PetscSectionGetDof(mesh->supportSection, p, size);
1261:   return(0);
1262: }

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

1269:   Not collective

1271:   Input Parameters:
1272: + mesh - The DMPlex
1273: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1274: - size - The support size for point p

1276:   Output Parameter:

1278:   Note:
1279:   This should be called after DMPlexSetChart().

1281:   Level: beginner

1283: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1284: @*/
1285: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1286: {
1287:   DM_Plex       *mesh = (DM_Plex*) dm->data;

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

1294:   mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1295:   return(0);
1296: }

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

1303:   Not collective

1305:   Input Parameters:
1306: + mesh - The DMPlex
1307: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1312:   Level: beginner

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

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

1320: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1321: @*/
1322: PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1323: {
1324:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1325:   PetscInt       off;

1331:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1332:   *support = &mesh->supports[off];
1333:   return(0);
1334: }

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

1341:   Not collective

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

1348:   Output Parameter:

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

1353:   Level: beginner

1355: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1356: @*/
1357: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1358: {
1359:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1360:   PetscInt       pStart, pEnd;
1361:   PetscInt       dof, off, c;

1366:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1367:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1369:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1370:   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);
1371:   for (c = 0; c < dof; ++c) {
1372:     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);
1373:     mesh->supports[off+c] = support[c];
1374:   }
1375:   return(0);
1376: }

1380: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1381: {
1382:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1383:   PetscInt       pStart, pEnd;
1384:   PetscInt       dof, off;

1389:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1390:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1391:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1392:   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);
1393:   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);
1394:   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);
1395:   mesh->supports[off+supportPos] = supportPoint;
1396:   return(0);
1397: }

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

1404:   Not collective

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

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

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

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

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

1425:   Level: beginner

1427: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1428: @*/
1429: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1430: {
1431:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1432:   PetscInt       *closure, *fifo;
1433:   const PetscInt *tmp = NULL, *tmpO = NULL;
1434:   PetscInt        tmpSize, t;
1435:   PetscInt        depth       = 0, maxSize;
1436:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1437:   PetscErrorCode  ierr;

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

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

1487:     closure[closureSize]   = cp;
1488:     closure[closureSize+1] = co;
1489:     fifo[fifoSize]         = cp;
1490:     fifo[fifoSize+1]       = co;
1491:   }
1492:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1493:   while (fifoSize - fifoStart) {
1494:     const PetscInt q   = fifo[fifoStart];
1495:     const PetscInt o   = fifo[fifoStart+1];
1496:     const PetscInt rev = o >= 0 ? 0 : 1;
1497:     const PetscInt off = rev ? -(o+1) : o;

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

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

1546: /*@C
1547:   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

1549:   Not collective

1551:   Input Parameters:
1552: + mesh - The DMPlex
1553: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1554: . orientation - The orientation of the point
1555: . useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
1556: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used

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

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

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

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

1571:   Level: beginner

1573: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1574: @*/
1575: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1576: {
1577:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1578:   PetscInt       *closure, *fifo;
1579:   const PetscInt *tmp = NULL, *tmpO = NULL;
1580:   PetscInt        tmpSize, t;
1581:   PetscInt        depth       = 0, maxSize;
1582:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1583:   PetscErrorCode  ierr;

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

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

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

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

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

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

1703:   Not collective

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

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

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

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

1721:   Level: beginner

1723: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1724: @*/
1725: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1726: {

1733:   DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1734:   if (numPoints) *numPoints = 0;
1735:   return(0);
1736: }

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

1743:   Not collective

1745:   Input Parameter:
1746: . mesh - The DMPlex

1748:   Output Parameters:
1749: + maxConeSize - The maximum number of in-edges
1750: - maxSupportSize - The maximum number of out-edges

1752:   Level: beginner

1754: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1755: @*/
1756: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1757: {
1758:   DM_Plex *mesh = (DM_Plex*) dm->data;

1762:   if (maxConeSize)    *maxConeSize    = mesh->maxConeSize;
1763:   if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1764:   return(0);
1765: }

1769: PetscErrorCode DMSetUp_Plex(DM dm)
1770: {
1771:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1772:   PetscInt       size;

1777:   PetscSectionSetUp(mesh->coneSection);
1778:   PetscSectionGetStorageSize(mesh->coneSection, &size);
1779:   PetscMalloc1(size, &mesh->cones);
1780:   PetscCalloc1(size, &mesh->coneOrientations);
1781:   if (mesh->maxSupportSize) {
1782:     PetscSectionSetUp(mesh->supportSection);
1783:     PetscSectionGetStorageSize(mesh->supportSection, &size);
1784:     PetscMalloc1(size, &mesh->supports);
1785:   }
1786:   return(0);
1787: }

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

1796:   if (subdm) {DMClone(dm, subdm);}
1797:   DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1798:   return(0);
1799: }

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

1806:   Not collective

1808:   Input Parameter:
1809: . mesh - The DMPlex

1811:   Output Parameter:

1813:   Note:
1814:   This should be called after all calls to DMPlexSetCone()

1816:   Level: beginner

1818: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1819: @*/
1820: PetscErrorCode DMPlexSymmetrize(DM dm)
1821: {
1822:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1823:   PetscInt      *offsets;
1824:   PetscInt       supportSize;
1825:   PetscInt       pStart, pEnd, p;

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

1836:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1837:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1838:     for (c = off; c < off+dof; ++c) {
1839:       PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1840:     }
1841:   }
1842:   for (p = pStart; p < pEnd; ++p) {
1843:     PetscInt dof;

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

1847:     mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1848:   }
1849:   PetscSectionSetUp(mesh->supportSection);
1850:   /* Calculate supports */
1851:   PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1852:   PetscMalloc1(supportSize, &mesh->supports);
1853:   PetscCalloc1(pEnd - pStart, &offsets);
1854:   for (p = pStart; p < pEnd; ++p) {
1855:     PetscInt dof, off, c;

1857:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1858:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1859:     for (c = off; c < off+dof; ++c) {
1860:       const PetscInt q = mesh->cones[c];
1861:       PetscInt       offS;

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

1865:       mesh->supports[offS+offsets[q]] = p;
1866:       ++offsets[q];
1867:     }
1868:   }
1869:   PetscFree(offsets);
1870:   return(0);
1871: }

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

1881:   Not collective

1883:   Input Parameter:
1884: . mesh - The DMPlex

1886:   Output Parameter:

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

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

1896:   Level: beginner

1898: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1899: @*/
1900: PetscErrorCode DMPlexStratify(DM dm)
1901: {
1902:   DMLabel        label;
1903:   PetscInt       pStart, pEnd, p;
1904:   PetscInt       numRoots = 0, numLeaves = 0;

1909:   PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1910:   /* Calculate depth */
1911:   DMPlexGetChart(dm, &pStart, &pEnd);
1912:   DMPlexCreateLabel(dm, "depth");
1913:   DMPlexGetDepthLabel(dm, &label);
1914:   /* Initialize roots and count leaves */
1915:   for (p = pStart; p < pEnd; ++p) {
1916:     PetscInt coneSize, supportSize;

1918:     DMPlexGetConeSize(dm, p, &coneSize);
1919:     DMPlexGetSupportSize(dm, p, &supportSize);
1920:     if (!coneSize && supportSize) {
1921:       ++numRoots;
1922:       DMLabelSetValue(label, p, 0);
1923:     } else if (!supportSize && coneSize) {
1924:       ++numLeaves;
1925:     } else if (!supportSize && !coneSize) {
1926:       /* Isolated points */
1927:       DMLabelSetValue(label, p, 0);
1928:     }
1929:   }
1930:   if (numRoots + numLeaves == (pEnd - pStart)) {
1931:     for (p = pStart; p < pEnd; ++p) {
1932:       PetscInt coneSize, supportSize;

1934:       DMPlexGetConeSize(dm, p, &coneSize);
1935:       DMPlexGetSupportSize(dm, p, &supportSize);
1936:       if (!supportSize && coneSize) {
1937:         DMLabelSetValue(label, p, 1);
1938:       }
1939:     }
1940:   } else {
1941:     IS       pointIS;
1942:     PetscInt numPoints = 0, level = 0;

1944:     DMLabelGetStratumIS(label, level, &pointIS);
1945:     if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1946:     while (numPoints) {
1947:       const PetscInt *points;
1948:       const PetscInt  newLevel = level+1;

1950:       ISGetIndices(pointIS, &points);
1951:       for (p = 0; p < numPoints; ++p) {
1952:         const PetscInt  point = points[p];
1953:         const PetscInt *support;
1954:         PetscInt        supportSize, s;

1956:         DMPlexGetSupportSize(dm, point, &supportSize);
1957:         DMPlexGetSupport(dm, point, &support);
1958:         for (s = 0; s < supportSize; ++s) {
1959:           DMLabelSetValue(label, support[s], newLevel);
1960:         }
1961:       }
1962:       ISRestoreIndices(pointIS, &points);
1963:       ++level;
1964:       ISDestroy(&pointIS);
1965:       DMLabelGetStratumIS(label, level, &pointIS);
1966:       if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1967:       else         {numPoints = 0;}
1968:     }
1969:     ISDestroy(&pointIS);
1970:   }
1971:   PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1972:   return(0);
1973: }

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

1980:   Not Collective

1982:   Input Parameters:
1983: + dm - The DMPlex object
1984: . numPoints - The number of input points for the join
1985: - points - The input points

1987:   Output Parameters:
1988: + numCoveredPoints - The number of points in the join
1989: - coveredPoints - The points in the join

1991:   Level: intermediate

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

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

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

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

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

2029:     PetscSectionGetDof(mesh->supportSection, points[p], &dof);
2030:     PetscSectionGetOffset(mesh->supportSection, points[p], &off);
2031:     for (c = 0; c < dof; ++c) {
2032:       const PetscInt point = mesh->supports[off+c];

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

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

2055:   Not Collective

2057:   Input Parameters:
2058: + dm - The DMPlex object
2059: . numPoints - The number of input points for the join
2060: - points - The input points

2062:   Output Parameters:
2063: + numCoveredPoints - The number of points in the join
2064: - coveredPoints - The points in the join

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

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

2072:   Level: intermediate

2074: .keywords: mesh
2075: .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet()
2076: @*/
2077: PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2078: {

2086:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2087:   if (numCoveredPoints) *numCoveredPoints = 0;
2088:   return(0);
2089: }

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

2096:   Not Collective

2098:   Input Parameters:
2099: + dm - The DMPlex object
2100: . numPoints - The number of input points for the join
2101: - points - The input points

2103:   Output Parameters:
2104: + numCoveredPoints - The number of points in the join
2105: - coveredPoints - The points in the join

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

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

2113:   Level: intermediate

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


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

2141:   for (p = 0; p < numPoints; ++p) {
2142:     PetscInt closureSize;

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

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

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

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

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

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

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

2205:   Not Collective

2207:   Input Parameters:
2208: + dm - The DMPlex object
2209: . numPoints - The number of input points for the meet
2210: - points - The input points

2212:   Output Parameters:
2213: + numCoveredPoints - The number of points in the meet
2214: - coveredPoints - The points in the meet

2216:   Level: intermediate

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

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

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

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

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

2254:     PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2255:     PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2256:     for (c = 0; c < dof; ++c) {
2257:       const PetscInt point = mesh->cones[off+c];

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

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

2280:   Not Collective

2282:   Input Parameters:
2283: + dm - The DMPlex object
2284: . numPoints - The number of input points for the meet
2285: - points - The input points

2287:   Output Parameters:
2288: + numCoveredPoints - The number of points in the meet
2289: - coveredPoints - The points in the meet

2291:   Level: intermediate

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

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

2299: .keywords: mesh
2300: .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin()
2301: @*/
2302: PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2303: {

2311:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2312:   if (numCoveredPoints) *numCoveredPoints = 0;
2313:   return(0);
2314: }

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

2321:   Not Collective

2323:   Input Parameters:
2324: + dm - The DMPlex object
2325: . numPoints - The number of input points for the meet
2326: - points - The input points

2328:   Output Parameters:
2329: + numCoveredPoints - The number of points in the meet
2330: - coveredPoints - The points in the meet

2332:   Level: intermediate

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

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

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


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

2366:   for (p = 0; p < numPoints; ++p) {
2367:     PetscInt closureSize;

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

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

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

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

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

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

2427: /*@C
2428:   DMPlexEqual - Determine if two DMs have the same topology

2430:   Not Collective

2432:   Input Parameters:
2433: + dmA - A DMPlex object
2434: - dmB - A DMPlex object

2436:   Output Parameters:
2437: . equal - PETSC_TRUE if the topologies are identical

2439:   Level: intermediate

2441:   Notes:
2442:   We are not solving graph isomorphism, so we do not permutation.

2444: .keywords: mesh
2445: .seealso: DMPlexGetCone()
2446: @*/
2447: PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal)
2448: {
2449:   PetscInt       depth, depthB, pStart, pEnd, pStartB, pEndB, p;


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

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

2494: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2495: {
2496:   MPI_Comm       comm;

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

2565: /*@
2566:   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.

2568:   Input Parameters:
2569: + dm     - The DM
2570: - in     - The input coordinate point (dim numbers)

2572:   Output Parameter:
2573: . out - The localized coordinate point

2575:   Level: developer

2577: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeAddCoordinate()
2578: @*/
2579: PetscErrorCode DMPlexLocalizeCoordinate(DM dm, const PetscScalar in[], PetscScalar out[])
2580: {
2581:   PetscInt       dim, d;

2585:   DMGetCoordinateDim(dm, &dim);
2586:   if (!dm->maxCell) {
2587:     for (d = 0; d < dim; ++d) out[d] = in[d];
2588:   } else {
2589:     for (d = 0; d < dim; ++d) {
2590:       out[d] = in[d] - dm->L[d]*floor(PetscRealPart(in[d])/dm->L[d]);
2591:     }
2592:   }
2593:   return(0);
2594: }

2598: /*
2599:   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.

2601:   Input Parameters:
2602: + dm     - The DM
2603: . dim    - The spatial dimension
2604: . anchor - The anchor point, the input point can be no more than maxCell away from it
2605: - in     - The input coordinate point (dim numbers)

2607:   Output Parameter:
2608: . out - The localized coordinate point

2610:   Level: developer

2612:   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

2614: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeAddCoordinate()
2615: */
2616: PetscErrorCode DMPlexLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2617: {
2618:   PetscInt d;

2621:   if (!dm->maxCell) {
2622:     for (d = 0; d < dim; ++d) out[d] = in[d];
2623:   } else {
2624:     for (d = 0; d < dim; ++d) {
2625:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2626:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2627:       } else {
2628:         out[d] = in[d];
2629:       }
2630:     }
2631:   }
2632:   return(0);
2633: }
2636: PetscErrorCode DMPlexLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
2637: {
2638:   PetscInt d;

2641:   if (!dm->maxCell) {
2642:     for (d = 0; d < dim; ++d) out[d] = in[d];
2643:   } else {
2644:     for (d = 0; d < dim; ++d) {
2645:       if (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d]) {
2646:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
2647:       } else {
2648:         out[d] = in[d];
2649:       }
2650:     }
2651:   }
2652:   return(0);
2653: }

2657: /*
2658:   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.

2660:   Input Parameters:
2661: + dm     - The DM
2662: . dim    - The spatial dimension
2663: . anchor - The anchor point, the input point can be no more than maxCell away from it
2664: . in     - The input coordinate delta (dim numbers)
2665: - out    - The input coordinate point (dim numbers)

2667:   Output Parameter:
2668: . out    - The localized coordinate in + out

2670:   Level: developer

2672:   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

2674: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeCoordinate()
2675: */
2676: PetscErrorCode DMPlexLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2677: {
2678:   PetscInt d;

2681:   if (!dm->maxCell) {
2682:     for (d = 0; d < dim; ++d) out[d] += in[d];
2683:   } else {
2684:     for (d = 0; d < dim; ++d) {
2685:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2686:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2687:       } else {
2688:         out[d] += in[d];
2689:       }
2690:     }
2691:   }
2692:   return(0);
2693: }

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

2700:   Input Parameter:
2701: . dm - The DM

2703:   Level: developer

2705: .seealso: DMPlexLocalizeCoordinate(), DMPlexLocalizeAddCoordinate()
2706: @*/
2707: PetscErrorCode DMPlexLocalizeCoordinates(DM dm)
2708: {
2709:   PetscSection   coordSection, cSection;
2710:   Vec            coordinates,  cVec;
2711:   PetscScalar   *coords, *coords2, *anchor;
2712:   PetscInt       Nc, cStart, cEnd, c, vStart, vEnd, v, dof, d, off, off2, bs, coordSize;

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

2757:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2758:     PetscSectionGetOffset(cSection, c, &off2);
2759:     for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
2760:     for (d = 0; d < dof/bs; ++d) {DMPlexLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
2761:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2762:   }
2763:   DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);
2764:   VecRestoreArray(coordinates, &coords);
2765:   VecRestoreArray(cVec,        &coords2);
2766:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
2767:   DMSetCoordinatesLocal(dm, cVec);
2768:   VecDestroy(&cVec);
2769:   PetscSectionDestroy(&cSection);
2770:   return(0);
2771: }

2775: /*@
2776:   DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point

2778:   Not Collective

2780:   Input Parameter:
2781: . dm    - The DMPlex object

2783:   Output Parameter:
2784: . depthLabel - The DMLabel recording point depth

2786:   Level: developer

2788: .keywords: mesh, points
2789: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2790: @*/
2791: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
2792: {
2793:   DM_Plex       *mesh = (DM_Plex*) dm->data;

2799:   if (!mesh->depthLabel) {DMPlexGetLabel(dm, "depth", &mesh->depthLabel);}
2800:   *depthLabel = mesh->depthLabel;
2801:   return(0);
2802: }

2806: /*@
2807:   DMPlexGetDepth - Get the depth of the DAG representing this mesh

2809:   Not Collective

2811:   Input Parameter:
2812: . dm    - The DMPlex object

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

2817:   Level: developer

2819: .keywords: mesh, points
2820: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2821: @*/
2822: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
2823: {
2824:   DMLabel        label;
2825:   PetscInt       d = 0;

2831:   DMPlexGetDepthLabel(dm, &label);
2832:   if (label) {DMLabelGetNumValues(label, &d);}
2833:   *depth = d-1;
2834:   return(0);
2835: }

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

2842:   Not Collective

2844:   Input Parameters:
2845: + dm           - The DMPlex object
2846: - stratumValue - The requested depth

2848:   Output Parameters:
2849: + start - The first point at this depth
2850: - end   - One beyond the last point at this depth

2852:   Level: developer

2854: .keywords: mesh, points
2855: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
2856: @*/
2857: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2858: {
2859:   DMLabel        label;
2860:   PetscInt       pStart, pEnd;

2867:   DMPlexGetChart(dm, &pStart, &pEnd);
2868:   if (pStart == pEnd) return(0);
2869:   if (stratumValue < 0) {
2870:     if (start) *start = pStart;
2871:     if (end)   *end   = pEnd;
2872:     return(0);
2873:   }
2874:   DMPlexGetDepthLabel(dm, &label);
2875:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2876:   DMLabelGetStratumBounds(label, stratumValue, start, end);
2877:   return(0);
2878: }

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

2885:   Not Collective

2887:   Input Parameters:
2888: + dm           - The DMPlex object
2889: - stratumValue - The requested height

2891:   Output Parameters:
2892: + start - The first point at this height
2893: - end   - One beyond the last point at this height

2895:   Level: developer

2897: .keywords: mesh, points
2898: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
2899: @*/
2900: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2901: {
2902:   DMLabel        label;
2903:   PetscInt       depth, pStart, pEnd;

2910:   DMPlexGetChart(dm, &pStart, &pEnd);
2911:   if (pStart == pEnd) return(0);
2912:   if (stratumValue < 0) {
2913:     if (start) *start = pStart;
2914:     if (end)   *end   = pEnd;
2915:     return(0);
2916:   }
2917:   DMPlexGetDepthLabel(dm, &label);
2918:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2919:   DMLabelGetNumValues(label, &depth);
2920:   DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
2921:   return(0);
2922: }

2926: /* Set the number of dof on each point and separate by fields */
2927: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
2928: {
2929:   PetscInt      *pMax;
2930:   PetscInt       depth, pStart = 0, pEnd = 0;
2931:   PetscInt       Nf, p, d, dep, f;
2932:   PetscBool     *isFE;

2936:   PetscMalloc1(numFields, &isFE);
2937:   DMGetNumFields(dm, &Nf);
2938:   for (f = 0; f < numFields; ++f) {
2939:     PetscObject  obj;
2940:     PetscClassId id;

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

2970:       for (f = 0; f < numFields; ++f) {
2971:         if (isFE[f] && p >= pMax[dep]) continue;
2972:         PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
2973:         tot += numDof[f*(dim+1)+d];
2974:       }
2975:       PetscSectionSetDof(*section, p, tot);
2976:     }
2977:   }
2978:   PetscFree(pMax);
2979:   PetscFree(isFE);
2980:   return(0);
2981: }

2985: /* Set the number of dof on each point and separate by fields
2986:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
2987: */
2988: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2989: {
2990:   PetscInt       numFields;
2991:   PetscInt       bc;
2992:   PetscSection   aSec;

2996:   PetscSectionGetNumFields(section, &numFields);
2997:   for (bc = 0; bc < numBC; ++bc) {
2998:     PetscInt        field = 0;
2999:     const PetscInt *comp;
3000:     const PetscInt *idx;
3001:     PetscInt        Nc = -1, n, i;

3003:     if (numFields) field = bcField[bc];
3004:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
3005:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
3006:     ISGetLocalSize(bcPoints[bc], &n);
3007:     ISGetIndices(bcPoints[bc], &idx);
3008:     for (i = 0; i < n; ++i) {
3009:       const PetscInt p = idx[i];
3010:       PetscInt       numConst;

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

3029:     PetscSectionGetChart(aSec, &aStart, &aEnd);
3030:     for (a = aStart; a < aEnd; a++) {
3031:       PetscInt dof, f;

3033:       PetscSectionGetDof(aSec, a, &dof);
3034:       if (dof) {
3035:         /* if there are point-to-point constraints, then all dofs are constrained */
3036:         PetscSectionGetDof(section, a, &dof);
3037:         PetscSectionSetConstraintDof(section, a, dof);
3038:         for (f = 0; f < numFields; f++) {
3039:           PetscSectionGetFieldDof(section, a, f, &dof);
3040:           PetscSectionSetFieldConstraintDof(section, a, f, dof);
3041:         }
3042:       }
3043:     }
3044:   }
3045:   return(0);
3046: }

3050: /* Set the constrained field indices on each point
3051:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
3052: */
3053: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
3054: {
3055:   PetscSection   aSec;
3056:   PetscInt      *indices;
3057:   PetscInt       numFields, maxDof, pStart, pEnd, p, bc, f, d;

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

3075:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
3076:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
3077:     ISGetLocalSize(bcPoints[bc], &n);
3078:     ISGetIndices(bcPoints[bc], &idx);
3079:     for (i = 0; i < n; ++i) {
3080:       const PetscInt  p = idx[i];
3081:       const PetscInt *find;
3082:       PetscInt        fcdof, c;

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

3104:     for (d = 0; d < maxDof; ++d) indices[d] = d;
3105:     PetscSectionGetChart(aSec, &aStart, &aEnd);
3106:     for (a = aStart; a < aEnd; a++) {
3107:       PetscInt dof, fdof, f;

3109:       PetscSectionGetDof(aSec, a, &dof);
3110:       if (dof) {
3111:         /* if there are point-to-point constraints, then all dofs are constrained */
3112:         for (f = 0; f < numFields; f++) {
3113:           PetscSectionGetFieldDof(section, a, f, &fdof);
3114:           PetscSectionSetFieldConstraintIndices(section, a, f, indices);
3115:         }
3116:       }
3117:     }
3118:   }
3119:   PetscFree(indices);
3120:   return(0);
3121: }

3125: /* Set the constrained indices on each point */
3126: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
3127: {
3128:   PetscInt      *indices;
3129:   PetscInt       numFields, maxDof, pStart, pEnd, p, f, d;

3133:   PetscSectionGetNumFields(section, &numFields);
3134:   PetscSectionGetMaxDof(section, &maxDof);
3135:   PetscSectionGetChart(section, &pStart, &pEnd);
3136:   PetscMalloc1(maxDof, &indices);
3137:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
3138:   for (p = pStart; p < pEnd; ++p) {
3139:     PetscInt cdof, d;

3141:     PetscSectionGetConstraintDof(section, p, &cdof);
3142:     if (cdof) {
3143:       if (numFields) {
3144:         PetscInt numConst = 0, foff = 0;

3146:         for (f = 0; f < numFields; ++f) {
3147:           const PetscInt *find;
3148:           PetscInt        fcdof, fdof;

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

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

3174:   Not Collective

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

3188:   Output Parameter:
3189: . section - The PetscSection object

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

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

3196:   Level: developer

3198:   Fortran Notes:
3199:   A Fortran 90 version is available as DMPlexCreateSectionF90()

3201: .keywords: mesh, elements
3202: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
3203: @*/
3204: 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)
3205: {
3206:   PetscSection   aSec;

3210:   DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
3211:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
3212:   if (perm) {PetscSectionSetPermutation(*section, perm);}
3213:   PetscSectionSetUp(*section);
3214:   DMPlexGetAnchors(dm,&aSec,NULL);
3215:   if (numBC || aSec) {
3216:     DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
3217:     DMPlexCreateSectionBCIndices(dm, *section);
3218:   }
3219:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
3220:   return(0);
3221: }

3225: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
3226: {
3227:   PetscSection   section;

3231:   DMClone(dm, cdm);
3232:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
3233:   DMSetDefaultSection(*cdm, section);
3234:   PetscSectionDestroy(&section);
3235:   return(0);
3236: }

3240: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
3241: {
3242:   DM_Plex *mesh = (DM_Plex*) dm->data;

3246:   if (section) *section = mesh->coneSection;
3247:   return(0);
3248: }

3252: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
3253: {
3254:   DM_Plex *mesh = (DM_Plex*) dm->data;

3258:   if (section) *section = mesh->supportSection;
3259:   return(0);
3260: }

3264: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
3265: {
3266:   DM_Plex *mesh = (DM_Plex*) dm->data;

3270:   if (cones) *cones = mesh->cones;
3271:   return(0);
3272: }

3276: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
3277: {
3278:   DM_Plex *mesh = (DM_Plex*) dm->data;

3282:   if (coneOrientations) *coneOrientations = mesh->coneOrientations;
3283:   return(0);
3284: }

3286: /******************************** FEM Support **********************************/

3290: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3291: {
3292:   PetscScalar    *array, *vArray;
3293:   const PetscInt *cone, *coneO;
3294:   PetscInt        pStart, pEnd, p, numPoints, size = 0, offset = 0;
3295:   PetscErrorCode  ierr;

3298:   PetscSectionGetChart(section, &pStart, &pEnd);
3299:   DMPlexGetConeSize(dm, point, &numPoints);
3300:   DMPlexGetCone(dm, point, &cone);
3301:   DMPlexGetConeOrientation(dm, point, &coneO);
3302:   if (!values || !*values) {
3303:     if ((point >= pStart) && (point < pEnd)) {
3304:       PetscInt dof;

3306:       PetscSectionGetDof(section, point, &dof);
3307:       size += dof;
3308:     }
3309:     for (p = 0; p < numPoints; ++p) {
3310:       const PetscInt cp = cone[p];
3311:       PetscInt       dof;

3313:       if ((cp < pStart) || (cp >= pEnd)) continue;
3314:       PetscSectionGetDof(section, cp, &dof);
3315:       size += dof;
3316:     }
3317:     if (!values) {
3318:       if (csize) *csize = size;
3319:       return(0);
3320:     }
3321:     DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
3322:   } else {
3323:     array = *values;
3324:   }
3325:   size = 0;
3326:   VecGetArray(v, &vArray);
3327:   if ((point >= pStart) && (point < pEnd)) {
3328:     PetscInt     dof, off, d;
3329:     PetscScalar *varr;

3331:     PetscSectionGetDof(section, point, &dof);
3332:     PetscSectionGetOffset(section, point, &off);
3333:     varr = &vArray[off];
3334:     for (d = 0; d < dof; ++d, ++offset) {
3335:       array[offset] = varr[d];
3336:     }
3337:     size += dof;
3338:   }
3339:   for (p = 0; p < numPoints; ++p) {
3340:     const PetscInt cp = cone[p];
3341:     PetscInt       o  = coneO[p];
3342:     PetscInt       dof, off, d;
3343:     PetscScalar   *varr;

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

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

3379:   *size = 0;
3380:   for (p = 0; p < numPoints*2; p += 2) {
3381:     const PetscInt point = points[p];
3382:     const PetscInt o     = points[p+1];
3383:     PetscInt       dof, off, d;
3384:     const PetscScalar *varr;

3386:     PetscSectionGetDof(section, point, &dof);
3387:     PetscSectionGetOffset(section, point, &off);
3388:     varr = &vArray[off];
3389:     if (o >= 0) {
3390:       for (d = 0; d < dof; ++d, ++offset)    array[offset] = varr[d];
3391:     } else {
3392:       for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
3393:     }
3394:   }
3395:   *size = offset;
3396:   return(0);
3397: }

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

3407:   *size = 0;
3408:   for (f = 0; f < numFields; ++f) {
3409:     PetscInt fcomp, p;

3411:     PetscSectionGetFieldComponents(section, f, &fcomp);
3412:     for (p = 0; p < numPoints*2; p += 2) {
3413:       const PetscInt point = points[p];
3414:       const PetscInt o     = points[p+1];
3415:       PetscInt       fdof, foff, d, c;
3416:       const PetscScalar *varr;

3418:       PetscSectionGetFieldDof(section, point, f, &fdof);
3419:       PetscSectionGetFieldOffset(section, point, f, &foff);
3420:       varr = &vArray[foff];
3421:       if (o >= 0) {
3422:         for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
3423:       } else {
3424:         for (d = fdof/fcomp-1; d >= 0; --d) {
3425:           for (c = 0; c < fcomp; ++c, ++offset) {
3426:             array[offset] = varr[d*fcomp+c];
3427:           }
3428:         }
3429:       }
3430:     }
3431:   }
3432:   *size = offset;
3433:   return(0);
3434: }

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

3441:   Not collective

3443:   Input Parameters:
3444: + dm - The DM
3445: . section - The section describing the layout in v, or NULL to use the default section
3446: . v - The local vector
3447: - point - The sieve point in the DM

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

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

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

3459:   Level: intermediate

3461: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3462: @*/
3463: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3464: {
3465:   PetscSection    clSection;
3466:   IS              clPoints;
3467:   PetscScalar    *array, *vArray;
3468:   PetscInt       *points = NULL;
3469:   const PetscInt *clp;
3470:   PetscInt        depth, numFields, numPoints, size;
3471:   PetscErrorCode  ierr;

3475:   if (!section) {DMGetDefaultSection(dm, &section);}
3478:   DMPlexGetDepth(dm, &depth);
3479:   PetscSectionGetNumFields(section, &numFields);
3480:   if (depth == 1 && numFields < 2) {
3481:     DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
3482:     return(0);
3483:   }
3484:   /* Get points */
3485:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3486:   if (!clPoints) {
3487:     PetscInt pStart, pEnd, p, q;

3489:     PetscSectionGetChart(section, &pStart, &pEnd);
3490:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3491:     /* Compress out points not in the section */
3492:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3493:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3494:         points[q*2]   = points[p];
3495:         points[q*2+1] = points[p+1];
3496:         ++q;
3497:       }
3498:     }
3499:     numPoints = q;
3500:   } else {
3501:     PetscInt dof, off;

3503:     PetscSectionGetDof(clSection, point, &dof);
3504:     PetscSectionGetOffset(clSection, point, &off);
3505:     ISGetIndices(clPoints, &clp);
3506:     numPoints = dof/2;
3507:     points    = (PetscInt *) &clp[off];
3508:   }
3509:   /* Get array */
3510:   if (!values || !*values) {
3511:     PetscInt asize = 0, dof, p;

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

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

3551:   Not collective

3553:   Input Parameters:
3554: + dm - The DM
3555: . section - The section describing the layout in v, or NULL to use the default section
3556: . v - The local vector
3557: . point - The sieve point in the DM
3558: . csize - The number of values in the closure, or NULL
3559: - values - The array of values, which is a borrowed array and should not be freed

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

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

3567:   Level: intermediate

3569: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3570: @*/
3571: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3572: {
3573:   PetscInt       size = 0;

3577:   /* Should work without recalculating size */
3578:   DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
3579:   return(0);
3580: }

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

3587: 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[])
3588: {
3589:   PetscInt        cdof;   /* The number of constraints on this point */
3590:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3591:   PetscScalar    *a;
3592:   PetscInt        off, cind = 0, k;
3593:   PetscErrorCode  ierr;

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

3628: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3629: {
3630:   PetscInt        cdof;   /* The number of constraints on this point */
3631:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3632:   PetscScalar    *a;
3633:   PetscInt        off, cind = 0, k;
3634:   PetscErrorCode  ierr;

3637:   PetscSectionGetConstraintDof(section, point, &cdof);
3638:   PetscSectionGetOffset(section, point, &off);
3639:   a    = &array[off];
3640:   if (cdof) {
3641:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3642:     if (orientation >= 0) {
3643:       for (k = 0; k < dof; ++k) {
3644:         if ((cind < cdof) && (k == cdofs[cind])) {
3645:           fuse(&a[k], values[k]);
3646:           ++cind;
3647:         }
3648:       }
3649:     } else {
3650:       for (k = 0; k < dof; ++k) {
3651:         if ((cind < cdof) && (k == cdofs[cind])) {
3652:           fuse(&a[k], values[dof-k-1]);
3653:           ++cind;
3654:         }
3655:       }
3656:     }
3657:   }
3658:   return(0);
3659: }

3663: 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[])
3664: {
3665:   PetscScalar    *a;
3666:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3667:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3668:   PetscInt        cind = 0, k, c;
3669:   PetscErrorCode  ierr;

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

3708: 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[])
3709: {
3710:   PetscScalar    *a;
3711:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3712:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3713:   PetscInt        cind = 0, k, c;
3714:   PetscErrorCode  ierr;

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

3747: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3748: {
3749:   PetscScalar    *array;
3750:   const PetscInt *cone, *coneO;
3751:   PetscInt        pStart, pEnd, p, numPoints, off, dof;
3752:   PetscErrorCode  ierr;

3755:   PetscSectionGetChart(section, &pStart, &pEnd);
3756:   DMPlexGetConeSize(dm, point, &numPoints);
3757:   DMPlexGetCone(dm, point, &cone);
3758:   DMPlexGetConeOrientation(dm, point, &coneO);
3759:   VecGetArray(v, &array);
3760:   for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
3761:     const PetscInt cp = !p ? point : cone[p-1];
3762:     const PetscInt o  = !p ? 0     : coneO[p-1];

3764:     if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
3765:     PetscSectionGetDof(section, cp, &dof);
3766:     /* ADD_VALUES */
3767:     {
3768:       const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3769:       PetscScalar    *a;
3770:       PetscInt        cdof, coff, cind = 0, k;

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

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

3810:   Not collective

3812:   Input Parameters:
3813: + dm - The DM
3814: . section - The section describing the layout in v, or NULL to use the default section
3815: . v - The local vector
3816: . point - The sieve point in the DM
3817: . values - The array of values
3818: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

3823:   Level: intermediate

3825: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
3826: @*/
3827: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3828: {
3829:   PetscSection    clSection;
3830:   IS              clPoints;
3831:   PetscScalar    *array;
3832:   PetscInt       *points = NULL;
3833:   const PetscInt *clp;
3834:   PetscInt        depth, numFields, numPoints, p;
3835:   PetscErrorCode  ierr;

3839:   if (!section) {DMGetDefaultSection(dm, &section);}
3842:   DMPlexGetDepth(dm, &depth);
3843:   PetscSectionGetNumFields(section, &numFields);
3844:   if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
3845:     DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
3846:     return(0);
3847:   }
3848:   /* Get points */
3849:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3850:   if (!clPoints) {
3851:     PetscInt pStart, pEnd, q;

3853:     PetscSectionGetChart(section, &pStart, &pEnd);
3854:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3855:     /* Compress out points not in the section */
3856:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3857:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3858:         points[q*2]   = points[p];
3859:         points[q*2+1] = points[p+1];
3860:         ++q;
3861:       }
3862:     }
3863:     numPoints = q;
3864:   } else {
3865:     PetscInt dof, off;

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

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

3963: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
3964: {
3965:   PetscSection    clSection;
3966:   IS              clPoints;
3967:   PetscScalar    *array;
3968:   PetscInt       *points = NULL;
3969:   const PetscInt *clp;
3970:   PetscInt        numFields, numPoints, p;
3971:   PetscInt        offset = 0, fcomp, f;
3972:   PetscErrorCode  ierr;

3976:   if (!section) {DMGetDefaultSection(dm, &section);}
3979:   PetscSectionGetNumFields(section, &numFields);
3980:   /* Get points */
3981:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3982:   if (!clPoints) {
3983:     PetscInt pStart, pEnd, q;

3985:     PetscSectionGetChart(section, &pStart, &pEnd);
3986:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3987:     /* Compress out points not in the section */
3988:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3989:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3990:         points[q*2]   = points[p];
3991:         points[q*2+1] = points[p+1];
3992:         ++q;
3993:       }
3994:     }
3995:     numPoints = q;
3996:   } else {
3997:     PetscInt dof, off;

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

4063: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
4064: {
4065:   PetscMPIInt    rank;
4066:   PetscInt       i, j;

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

4091: /* . off - The global offset of this point */
4092: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
4093: {
4094:   PetscInt        dof;    /* The number of unknowns on this point */
4095:   PetscInt        cdof;   /* The number of constraints on this point */
4096:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
4097:   PetscInt        cind = 0, k;
4098:   PetscErrorCode  ierr;

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

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

4146:   PetscSectionGetNumFields(section, &numFields);
4147:   for (f = 0, foff = 0; f < numFields; ++f) {
4148:     PetscInt        fdof, fcomp, cfdof;
4149:     const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4150:     PetscInt        cind = 0, k, c;

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

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

4219:   PetscSectionGetNumFields(section, &numFields);

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

4234:       if (b >= aStart && b < aEnd) {
4235:         PetscSectionGetDof(aSec,b,&bDof);
4236:       }
4237:       if (bDof) {
4238:         /* this point is constrained */
4239:         /* it is going to be replaced by its anchors */
4240:         PetscInt bOff, q;

4242:         anyConstrained = PETSC_TRUE;
4243:         newNumPoints  += bDof;
4244:         PetscSectionGetOffset(aSec,b,&bOff);
4245:         for (q = 0; q < bDof; q++) {
4246:           PetscInt a = anchors[bOff + q];
4247:           PetscInt aDof;

4249:           PetscSectionGetDof(section,a,&aDof);
4250:           newNumIndices += aDof;
4251:           for (f = 0; f < numFields; ++f) {
4252:             PetscInt fDof;

4254:             PetscSectionGetFieldDof(section, a, f, &fDof);
4255:             newOffsets[f+1] += fDof;
4256:           }
4257:         }
4258:       }
4259:       else {
4260:         /* this point is not constrained */
4261:         newNumPoints++;
4262:         PetscSectionGetDof(section,b,&bDof);
4263:         newNumIndices += bDof;
4264:         for (f = 0; f < numFields; ++f) {
4265:           PetscInt fDof;

4267:           PetscSectionGetFieldDof(section, b, f, &fDof);
4268:           newOffsets[f+1] += fDof;
4269:         }
4270:       }
4271:     }
4272:   }
4273:   if (!anyConstrained) {
4274:     if (outNumPoints)  *outNumPoints  = 0;
4275:     if (outNumIndices) *outNumIndices = 0;
4276:     if (outPoints)     *outPoints     = NULL;
4277:     if (outValues)     *outValues     = NULL;
4278:     if (aSec) {ISRestoreIndices(aIS,&anchors);}
4279:     return(0);
4280:   }

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

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

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

4288:   /* output arrays */
4289:   DMGetWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4290:   DMGetWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);

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

4305:   /* get workspaces for the point-to-point matrices */
4306:   if (numFields) {
4307:     for (p = 0; p < numPoints; p++) {
4308:       PetscInt b    = points[2*p];
4309:       PetscInt bDof = 0;

4311:       if (b >= aStart && b < aEnd) {
4312:         PetscSectionGetDof(aSec, b, &bDof);
4313:       }
4314:       if (bDof) {
4315:         for (f = 0; f < numFields; f++) {
4316:           PetscInt fDof, q, bOff, allFDof = 0;

4318:           PetscSectionGetFieldDof(section, b, f, &fDof);
4319:           PetscSectionGetOffset(aSec, b, &bOff);
4320:           for (q = 0; q < bDof; q++) {
4321:             PetscInt a = anchors[bOff + q];
4322:             PetscInt aFDof;

4324:             PetscSectionGetFieldDof(section, a, f, &aFDof);
4325:             allFDof += aFDof;
4326:           }
4327:           newPointOffsets[f][p+1] = allFDof;
4328:           pointMatOffsets[f][p+1] = fDof * allFDof;
4329:         }
4330:       }
4331:       else {
4332:         for (f = 0; f < numFields; f++) {
4333:           PetscInt fDof;

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

4356:       if (b >= aStart && b < aEnd) {
4357:         PetscSectionGetDof(aSec, b, &bDof);
4358:       }
4359:       if (bDof) {
4360:         PetscInt dof, bOff, q, allDof = 0;

4362:         PetscSectionGetDof(section, b, &dof);
4363:         PetscSectionGetOffset(aSec, b, &bOff);
4364:         for (q = 0; q < bDof; q++) {
4365:           PetscInt a = anchors[bOff + q], aDof;

4367:           PetscSectionGetDof(section, a, &aDof);
4368:           allDof += aDof;
4369:         }
4370:         newPointOffsets[0][p+1] = allDof;
4371:         pointMatOffsets[0][p+1] = dof * allDof;
4372:       }
4373:       else {
4374:         PetscInt dof;

4376:         PetscSectionGetDof(section, b, &dof);
4377:         newPointOffsets[0][p+1] = dof;
4378:         pointMatOffsets[0][p+1] = 0;
4379:       }
4380:     }
4381:     newPointOffsets[0][0] = 0;
4382:     pointMatOffsets[0][0] = 0;
4383:     for (p = 0; p < numPoints; p++) {
4384:       newPointOffsets[0][p+1] += newPointOffsets[0][p];
4385:       pointMatOffsets[0][p+1] += pointMatOffsets[0][p];
4386:     }
4387:     DMGetWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4388:   }

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

4401:       if (b >= aStart && b < aEnd) {
4402:         PetscSectionGetDof(aSec, b, &bDof);
4403:       }
4404:       if (bDof) {
4405:         PetscInt fStart[32], fEnd[32], fAnchorStart[32], fAnchorEnd[32], bOff, q;

4407:         fStart[0] = 0;
4408:         fEnd[0]   = 0;
4409:         for (f = 0; f < numFields; f++) {
4410:           PetscInt fDof;

4412:           PetscSectionGetFieldDof(cSec, b, f, &fDof);
4413:           fStart[f+1] = fStart[f] + fDof;
4414:           fEnd[f+1]   = fStart[f+1];
4415:         }
4416:         PetscSectionGetOffset(cSec, b, &bOff);
4417:         indicesPointFields_private(cSec, b, bOff, fEnd, PETSC_TRUE, o, indices);

4419:         fAnchorStart[0] = 0;
4420:         fAnchorEnd[0]   = 0;
4421:         for (f = 0; f < numFields; f++) {
4422:           PetscInt fDof = newPointOffsets[f][p + 1] - newPointOffsets[f][p];

4424:           fAnchorStart[f+1] = fAnchorStart[f] + fDof;
4425:           fAnchorEnd[f+1]   = fAnchorStart[f + 1];
4426:         }
4427:         PetscSectionGetOffset (aSec, b, &bOff);
4428:         for (q = 0; q < bDof; q++) {
4429:           PetscInt a = anchors[bOff + q], aOff;

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

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

4456:       if (b >= aStart && b < aEnd) {
4457:         PetscSectionGetDof(aSec, b, &bDof);
4458:       }
4459:       if (bDof) {
4460:         PetscInt bEnd = 0, bAnchorEnd = 0, bOff;

4462:         PetscSectionGetOffset(cSec, b, &bOff);
4463:         indicesPoint_private(cSec, b, bOff, &bEnd, PETSC_TRUE, o, indices);

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

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

4471:           newPoints[2*(newP + q)]     = a;
4472:           newPoints[2*(newP + q) + 1] = 0;
4473:           PetscSectionGetOffset(section, a, &aOff);
4474:           indicesPoint_private(section, a, aOff, &bAnchorEnd, PETSC_TRUE, 0, newIndices);
4475:         }
4476:         newP += bDof;

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

4489:   PetscMemzero(tmpValues,newNumIndices*numIndices*sizeof(*tmpValues));
4490:   /* multiply constraints on the right */
4491:   if (numFields) {
4492:     for (f = 0; f < numFields; f++) {
4493:       PetscInt oldOff = offsets[f];

4495:       for (p = 0; p < numPoints; p++) {
4496:         PetscInt cStart = newPointOffsets[f][p];
4497:         PetscInt b      = points[2 * p];
4498:         PetscInt c, r, k;
4499:         PetscInt dof;

4501:         PetscSectionGetFieldDof(section,b,f,&dof);
4502:         if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4503:           PetscInt nCols         = newPointOffsets[f][p+1]-cStart;
4504:           const PetscScalar *mat = pointMat[f] + pointMatOffsets[f][p];

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

4534:       PetscSectionGetDof(section,b,&dof);
4535:       if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4536:         PetscInt nCols         = newPointOffsets[0][p+1]-cStart;
4537:         const PetscScalar *mat = pointMat[0] + pointMatOffsets[0][p];

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

4559:   PetscMemzero(newValues,newNumIndices*newNumIndices*sizeof(*newValues));
4560:   /* multiply constraints transpose on the left */
4561:   if (numFields) {
4562:     for (f = 0; f < numFields; f++) {
4563:       PetscInt oldOff = offsets[f];

4565:       for (p = 0; p < numPoints; p++) {
4566:         PetscInt rStart = newPointOffsets[f][p];
4567:         PetscInt b      = points[2 * p];
4568:         PetscInt c, r, k;
4569:         PetscInt dof;

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

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

4599:     for (p = 0; p < numPoints; p++) {
4600:       PetscInt rStart = newPointOffsets[0][p];
4601:       PetscInt b      = points[2 * p];
4602:       PetscInt c, r, k;
4603:       PetscInt dof;

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

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

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

4648:   /* output */
4649:   *outNumPoints  = newNumPoints;
4650:   *outNumIndices = newNumIndices;
4651:   *outPoints     = newPoints;
4652:   *outValues     = newValues;
4653:   for (f = 0; f < numFields; f++) {
4654:     offsets[f] = newOffsets[f];
4655:   }
4656:   return(0);
4657: }

4661: PetscErrorCode DMPlexGetClosureIndices(DM dm, PetscSection section, PetscSection globalSection, PetscInt point, PetscInt *numIndices, PetscInt **indices)
4662: {
4663:   PetscSection    clSection;
4664:   IS              clPoints;
4665:   const PetscInt *clp;
4666:   PetscInt       *points = NULL, *pointsNew;
4667:   PetscInt        numPoints, numPointsNew;
4668:   PetscInt        offsets[32];
4669:   PetscInt        Nf, Nind, NindNew, off, globalOff, f, p;
4670:   PetscErrorCode  ierr;

4678:   PetscSectionGetNumFields(section, &Nf);
4679:   if (Nf > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", Nf);
4680:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
4681:   /* Get points in closure */
4682:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4683:   if (!clPoints) {
4684:     PetscInt pStart, pEnd, q;

4686:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4687:     /* Compress out points not in the section */
4688:     PetscSectionGetChart(section, &pStart, &pEnd);
4689:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4690:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4691:         points[q*2]   = points[p];
4692:         points[q*2+1] = points[p+1];
4693:         ++q;
4694:       }
4695:     }
4696:     numPoints = q;
4697:   } else {
4698:     PetscInt dof, off;

4700:     PetscSectionGetDof(clSection, point, &dof);
4701:     numPoints = dof/2;
4702:     PetscSectionGetOffset(clSection, point, &off);
4703:     ISGetIndices(clPoints, &clp);
4704:     points = (PetscInt *) &clp[off];
4705:   }
4706:   /* Get number of indices and indices per field */
4707:   for (p = 0, Nind = 0; p < numPoints*2; p += 2) {
4708:     PetscInt dof, fdof;

4710:     PetscSectionGetDof(section, points[p], &dof);
4711:     for (f = 0; f < Nf; ++f) {
4712:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
4713:       offsets[f+1] += fdof;
4714:     }
4715:     Nind += dof;
4716:   }
4717:   for (f = 1; f < Nf; ++f) offsets[f+1] += offsets[f];
4718:   if (Nf && offsets[Nf] != Nind) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[Nf], Nind);
4719:   /* Correct for hanging node constraints */
4720:   {
4721:     DMPlexAnchorsModifyMat(dm, section, numPoints, Nind, points, NULL, &numPointsNew, &NindNew, &pointsNew, NULL, offsets);
4722:     if (numPointsNew) {
4723:       if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4724:       else           {ISRestoreIndices(clPoints, &clp);}
4725:       numPoints = numPointsNew;
4726:       Nind      = NindNew;
4727:       points    = pointsNew;
4728:     }
4729:   }
4730:   /* Calculate indices */
4731:   DMGetWorkArray(dm, Nind, PETSC_INT, indices);
4732:   if (Nf) {
4733:     for (p = 0; p < numPoints*2; p += 2) {
4734:       PetscInt o = points[p+1];
4735:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4736:       indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, *indices);
4737:     }
4738:   } else {
4739:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
4740:       PetscInt o = points[p+1];
4741:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4742:       indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, *indices);
4743:     }
4744:   }
4745:   /* Cleanup points */
4746:   if (numPointsNew) {
4747:     DMRestoreWorkArray(dm, 2*numPointsNew, PETSC_INT, &pointsNew);
4748:   } else {
4749:     if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4750:     else           {ISRestoreIndices(clPoints, &clp);}
4751:   }
4752:   if (numIndices) *numIndices = Nind;
4753:   return(0);
4754: }

4758: PetscErrorCode DMPlexRestoreClosureIndices(DM dm, PetscSection section, PetscSection globalSection, PetscInt point, PetscInt *numIndices, PetscInt **indices)
4759: {

4765:   DMRestoreWorkArray(dm, 0, PETSC_INT, indices);
4766:   return(0);
4767: }

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

4774:   Not collective

4776:   Input Parameters:
4777: + dm - The DM
4778: . section - The section describing the layout in v, or NULL to use the default section
4779: . globalSection - The section describing the layout in v, or NULL to use the default global section
4780: . A - The matrix
4781: . point - The sieve point in the DM
4782: . values - The array of values
4783: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

4788:   Level: intermediate

4790: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
4791: @*/
4792: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4793: {
4794:   DM_Plex        *mesh   = (DM_Plex*) dm->data;
4795:   PetscSection    clSection;
4796:   IS              clPoints;
4797:   PetscInt       *points = NULL, *newPoints;
4798:   const PetscInt *clp;
4799:   PetscInt       *indices;
4800:   PetscInt        offsets[32];
4801:   PetscInt        numFields, numPoints, newNumPoints, numIndices, newNumIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
4802:   PetscScalar    *newValues;
4803:   PetscErrorCode  ierr;

4807:   if (!section) {DMGetDefaultSection(dm, &section);}
4809:   if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
4812:   PetscSectionGetNumFields(section, &numFields);
4813:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4814:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
4815:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4816:   if (!clPoints) {
4817:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4818:     /* Compress out points not in the section */
4819:     PetscSectionGetChart(section, &pStart, &pEnd);
4820:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4821:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4822:         points[q*2]   = points[p];
4823:         points[q*2+1] = points[p+1];
4824:         ++q;
4825:       }
4826:     }
4827:     numPoints = q;
4828:   } else {
4829:     PetscInt dof, off;

4831:     PetscSectionGetDof(clSection, point, &dof);
4832:     numPoints = dof/2;
4833:     PetscSectionGetOffset(clSection, point, &off);
4834:     ISGetIndices(clPoints, &clp);
4835:     points = (PetscInt *) &clp[off];
4836:   }
4837:   for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
4838:     PetscInt fdof;

4840:     PetscSectionGetDof(section, points[p], &dof);
4841:     for (f = 0; f < numFields; ++f) {
4842:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
4843:       offsets[f+1] += fdof;
4844:     }
4845:     numIndices += dof;
4846:   }
4847:   for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];

4849:   if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
4850:   DMPlexAnchorsModifyMat(dm,section,numPoints,numIndices,points,values,&newNumPoints,&newNumIndices,&newPoints,&newValues,offsets);
4851:   if (newNumPoints) {
4852:     if (!clPoints) {
4853:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4854:     } else {
4855:       ISRestoreIndices(clPoints, &clp);
4856:     }
4857:     numPoints  = newNumPoints;
4858:     numIndices = newNumIndices;
4859:     points     = newPoints;
4860:     values     = newValues;
4861:   }
4862:   DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
4863:   if (numFields) {
4864:     for (p = 0; p < numPoints*2; p += 2) {
4865:       PetscInt o = points[p+1];
4866:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4867:       indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
4868:     }
4869:   } else {
4870:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
4871:       PetscInt o = points[p+1];
4872:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4873:       indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
4874:     }
4875:   }
4876:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
4877:   MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
4878:   if (mesh->printFEM > 1) {
4879:     PetscInt i;
4880:     PetscPrintf(PETSC_COMM_SELF, "  Indices:");
4881:     for (i = 0; i < numIndices; ++i) {PetscPrintf(PETSC_COMM_SELF, " %d", indices[i]);}
4882:     PetscPrintf(PETSC_COMM_SELF, "\n");
4883:   }
4884:   if (ierr) {
4885:     PetscMPIInt    rank;
4886:     PetscErrorCode ierr2;

4888:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4889:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4890:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
4891:     ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
4892: 
4893:   }
4894:   if (newNumPoints) {
4895:     DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4896:     DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4897:   }
4898:   else {
4899:     if (!clPoints) {
4900:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4901:     } else {
4902:       ISRestoreIndices(clPoints, &clp);
4903:     }
4904:   }
4905:   DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
4906:   return(0);
4907: }

4911: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4912: {
4913:   DM_Plex        *mesh   = (DM_Plex*) dmf->data;
4914:   PetscInt       *fpoints = NULL, *ftotpoints = NULL;
4915:   PetscInt       *cpoints = NULL;
4916:   PetscInt       *findices, *cindices;
4917:   PetscInt        foffsets[32], coffsets[32];
4918:   CellRefiner     cellRefiner;
4919:   PetscInt        numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4920:   PetscErrorCode  ierr;

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

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

4990:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4991:     if (!dof) continue;
4992:     for (f = 0; f < numFields; ++f) {
4993:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4994:       foffsets[f+1] += fdof;
4995:     }
4996:     numFIndices += dof;
4997:   }
4998:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

5000:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
5001:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
5002:   DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
5003:   DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
5004:   if (numFields) {
5005:     for (p = 0; p < numFPoints*2; p += 2) {
5006:       PetscInt o = ftotpoints[p+1];
5007:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5008:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
5009:     }
5010:     for (p = 0; p < numCPoints*2; p += 2) {
5011:       PetscInt o = cpoints[p+1];
5012:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5013:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5014:     }
5015:   } else {
5016:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5017:       PetscInt o = ftotpoints[p+1];
5018:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5019:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5020:     }
5021:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5022:       PetscInt o = cpoints[p+1];
5023:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5024:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5025:     }
5026:   }
5027:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
5028:   MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
5029:   if (ierr) {
5030:     PetscMPIInt    rank;
5031:     PetscErrorCode ierr2;

5033:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
5034:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
5035:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
5036:     ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
5037:     ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
5038: 
5039:   }
5040:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5041:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5042:   DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
5043:   DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
5044:   return(0);
5045: }

5049: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
5050: {
5051:   PetscInt      *fpoints = NULL, *ftotpoints = NULL;
5052:   PetscInt      *cpoints = NULL;
5053:   PetscInt       foffsets[32], coffsets[32];
5054:   CellRefiner    cellRefiner;
5055:   PetscInt       numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;

5061:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
5063:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
5065:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
5067:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
5069:   PetscSectionGetNumFields(fsection, &numFields);
5070:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
5071:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
5072:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
5073:   /* Column indices */
5074:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5075:   maxFPoints = numCPoints;
5076:   /* Compress out points not in the section */
5077:   /*   TODO: Squeeze out points with 0 dof as well */
5078:   PetscSectionGetChart(csection, &pStart, &pEnd);
5079:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
5080:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
5081:       cpoints[q*2]   = cpoints[p];
5082:       cpoints[q*2+1] = cpoints[p+1];
5083:       ++q;
5084:     }
5085:   }
5086:   numCPoints = q;
5087:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
5088:     PetscInt fdof;

5090:     PetscSectionGetDof(csection, cpoints[p], &dof);
5091:     if (!dof) continue;
5092:     for (f = 0; f < numFields; ++f) {
5093:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
5094:       coffsets[f+1] += fdof;
5095:     }
5096:     numCIndices += dof;
5097:   }
5098:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
5099:   /* Row indices */
5100:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
5101:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
5102:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
5103:   for (r = 0, q = 0; r < numSubcells; ++r) {
5104:     /* TODO Map from coarse to fine cells */
5105:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
5106:     /* Compress out points not in the section */
5107:     PetscSectionGetChart(fsection, &pStart, &pEnd);
5108:     for (p = 0; p < numFPoints*2; p += 2) {
5109:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
5110:         PetscSectionGetDof(fsection, fpoints[p], &dof);
5111:         if (!dof) continue;
5112:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
5113:         if (s < q) continue;
5114:         ftotpoints[q*2]   = fpoints[p];
5115:         ftotpoints[q*2+1] = fpoints[p+1];
5116:         ++q;
5117:       }
5118:     }
5119:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
5120:   }
5121:   numFPoints = q;
5122:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
5123:     PetscInt fdof;

5125:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
5126:     if (!dof) continue;
5127:     for (f = 0; f < numFields; ++f) {
5128:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
5129:       foffsets[f+1] += fdof;
5130:     }
5131:     numFIndices += dof;
5132:   }
5133:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

5135:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
5136:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
5137:   if (numFields) {
5138:     for (p = 0; p < numFPoints*2; p += 2) {
5139:       PetscInt o = ftotpoints[p+1];
5140:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5141:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
5142:     }
5143:     for (p = 0; p < numCPoints*2; p += 2) {
5144:       PetscInt o = cpoints[p+1];
5145:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5146:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5147:     }
5148:   } else {
5149:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5150:       PetscInt o = ftotpoints[p+1];
5151:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5152:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5153:     }
5154:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5155:       PetscInt o = cpoints[p+1];
5156:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5157:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5158:     }
5159:   }
5160:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5161:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5162:   return(0);
5163: }

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

5170:   Input Parameter:
5171: . dm - The DMPlex object

5173:   Output Parameters:
5174: + cMax - The first hybrid cell
5175: . fMax - The first hybrid face
5176: . eMax - The first hybrid edge
5177: - vMax - The first hybrid vertex

5179:   Level: developer

5181: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
5182: @*/
5183: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
5184: {
5185:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5186:   PetscInt       dim;

5191:   DMGetDimension(dm, &dim);
5192:   if (cMax) *cMax = mesh->hybridPointMax[dim];
5193:   if (fMax) *fMax = mesh->hybridPointMax[dim-1];
5194:   if (eMax) *eMax = mesh->hybridPointMax[1];
5195:   if (vMax) *vMax = mesh->hybridPointMax[0];
5196:   return(0);
5197: }

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

5204:   Input Parameters:
5205: . dm   - The DMPlex object
5206: . cMax - The first hybrid cell
5207: . fMax - The first hybrid face
5208: . eMax - The first hybrid edge
5209: - vMax - The first hybrid vertex

5211:   Level: developer

5213: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5214: @*/
5215: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5216: {
5217:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5218:   PetscInt       dim;

5223:   DMGetDimension(dm, &dim);
5224:   if (cMax >= 0) mesh->hybridPointMax[dim]   = cMax;
5225:   if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5226:   if (eMax >= 0) mesh->hybridPointMax[1]     = eMax;
5227:   if (vMax >= 0) mesh->hybridPointMax[0]     = vMax;
5228:   return(0);
5229: }

5233: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5234: {
5235:   DM_Plex *mesh = (DM_Plex*) dm->data;

5240:   *cellHeight = mesh->vtkCellHeight;
5241:   return(0);
5242: }

5246: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5247: {
5248:   DM_Plex *mesh = (DM_Plex*) dm->data;

5252:   mesh->vtkCellHeight = cellHeight;
5253:   return(0);
5254: }

5258: /* We can easily have a form that takes an IS instead */
5259: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5260: {
5261:   PetscSection   section, globalSection;
5262:   PetscInt      *numbers, p;

5266:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
5267:   PetscSectionSetChart(section, pStart, pEnd);
5268:   for (p = pStart; p < pEnd; ++p) {
5269:     PetscSectionSetDof(section, p, 1);
5270:   }
5271:   PetscSectionSetUp(section);
5272:   PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
5273:   PetscMalloc1(pEnd - pStart, &numbers);
5274:   for (p = pStart; p < pEnd; ++p) {
5275:     PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5276:     if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5277:     else                       numbers[p-pStart] += shift;
5278:   }
5279:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5280:   if (globalSize) {
5281:     PetscLayout layout;
5282:     PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5283:     PetscLayoutGetSize(layout, globalSize);
5284:     PetscLayoutDestroy(&layout);
5285:   }
5286:   PetscSectionDestroy(&section);
5287:   PetscSectionDestroy(&globalSection);
5288:   return(0);
5289: }

5293: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5294: {
5295:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5296:   PetscInt       cellHeight, cStart, cEnd, cMax;

5301:   if (!mesh->globalCellNumbers) {
5302:     DMPlexGetVTKCellHeight(dm, &cellHeight);
5303:     DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5304:     DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5305:     if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
5306:     DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, &mesh->globalCellNumbers);
5307:   }
5308:   *globalCellNumbers = mesh->globalCellNumbers;
5309:   return(0);
5310: }

5314: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5315: {
5316:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5317:   PetscInt       vStart, vEnd, vMax;

5322:   if (!mesh->globalVertexNumbers) {
5323:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5324:     DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5325:     if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
5326:     DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, &mesh->globalVertexNumbers);
5327:   }
5328:   *globalVertexNumbers = mesh->globalVertexNumbers;
5329:   return(0);
5330: }

5334: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5335: {
5336:   IS             nums[4];
5337:   PetscInt       depths[4];
5338:   PetscInt       depth, d, shift = 0;

5343:   DMPlexGetDepth(dm, &depth);
5344:   /* For unstratified meshes use dim instead of depth */
5345:   if (depth < 0) {DMGetDimension(dm, &depth);}
5346:   depths[0] = depth; depths[1] = 0;
5347:   for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5348:   for (d = 0; d <= depth; ++d) {
5349:     PetscInt pStart, pEnd, gsize;

5351:     DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5352:     DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5353:     shift += gsize;
5354:   }
5355:   ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5356:   for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5357:   return(0);
5358: }


5363: /*@C
5364:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
5365:   the local section and an SF describing the section point overlap.

5367:   Input Parameters:
5368:   + s - The PetscSection for the local field layout
5369:   . sf - The SF describing parallel layout of the section points
5370:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
5371:   . label - The label specifying the points
5372:   - labelValue - The label stratum specifying the points

5374:   Output Parameter:
5375:   . gsection - The PetscSection for the global field layout

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

5379:   Level: developer

5381: .seealso: PetscSectionCreate()
5382: @*/
5383: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
5384: {
5385:   PetscInt      *neg = NULL, *tmpOff = NULL;
5386:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

5390:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
5391:   PetscSectionGetChart(s, &pStart, &pEnd);
5392:   PetscSectionSetChart(*gsection, pStart, pEnd);
5393:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
5394:   if (nroots >= 0) {
5395:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
5396:     PetscCalloc1(nroots, &neg);
5397:     if (nroots > pEnd-pStart) {
5398:       PetscCalloc1(nroots, &tmpOff);
5399:     } else {
5400:       tmpOff = &(*gsection)->atlasDof[-pStart];
5401:     }
5402:   }
5403:   /* Mark ghost points with negative dof */
5404:   for (p = pStart; p < pEnd; ++p) {
5405:     PetscInt value;

5407:     DMLabelGetValue(label, p, &value);
5408:     if (value != labelValue) continue;
5409:     PetscSectionGetDof(s, p, &dof);
5410:     PetscSectionSetDof(*gsection, p, dof);
5411:     PetscSectionGetConstraintDof(s, p, &cdof);
5412:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
5413:     if (neg) neg[p] = -(dof+1);
5414:   }
5415:   PetscSectionSetUpBC(*gsection);
5416:   if (nroots >= 0) {
5417:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5418:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5419:     if (nroots > pEnd-pStart) {
5420:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
5421:     }
5422:   }
5423:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
5424:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5425:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
5426:     (*gsection)->atlasOff[p] = off;
5427:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
5428:   }
5429:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
5430:   globalOff -= off;
5431:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5432:     (*gsection)->atlasOff[p] += globalOff;
5433:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
5434:   }
5435:   /* Put in negative offsets for ghost points */
5436:   if (nroots >= 0) {
5437:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5438:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5439:     if (nroots > pEnd-pStart) {
5440:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
5441:     }
5442:   }
5443:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
5444:   PetscFree(neg);
5445:   return(0);
5446: }

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

5453:   Input Parameters:
5454:   + dm - The DMPlex object

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

5458:   Level: developer

5460: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5461: @*/
5462: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5463: {
5464:   PetscSection    coneSection, supportSection;
5465:   const PetscInt *cone, *support;
5466:   PetscInt        coneSize, c, supportSize, s;
5467:   PetscInt        pStart, pEnd, p, csize, ssize;
5468:   PetscErrorCode  ierr;

5472:   DMPlexGetConeSection(dm, &coneSection);
5473:   DMPlexGetSupportSection(dm, &supportSection);
5474:   /* Check that point p is found in the support of its cone points, and vice versa */
5475:   DMPlexGetChart(dm, &pStart, &pEnd);
5476:   for (p = pStart; p < pEnd; ++p) {
5477:     DMPlexGetConeSize(dm, p, &coneSize);
5478:     DMPlexGetCone(dm, p, &cone);
5479:     for (c = 0; c < coneSize; ++c) {
5480:       PetscBool dup = PETSC_FALSE;
5481:       PetscInt  d;
5482:       for (d = c-1; d >= 0; --d) {
5483:         if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5484:       }
5485:       DMPlexGetSupportSize(dm, cone[c], &supportSize);
5486:       DMPlexGetSupport(dm, cone[c], &support);
5487:       for (s = 0; s < supportSize; ++s) {
5488:         if (support[s] == p) break;
5489:       }
5490:       if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5491:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5492:         for (s = 0; s < coneSize; ++s) {
5493:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
5494:         }
5495:         PetscPrintf(PETSC_COMM_SELF, "\n");
5496:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
5497:         for (s = 0; s < supportSize; ++s) {
5498:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
5499:         }
5500:         PetscPrintf(PETSC_COMM_SELF, "\n");
5501:         if (dup) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
5502:         else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
5503:       }
5504:     }
5505:     DMPlexGetSupportSize(dm, p, &supportSize);
5506:     DMPlexGetSupport(dm, p, &support);
5507:     for (s = 0; s < supportSize; ++s) {
5508:       DMPlexGetConeSize(dm, support[s], &coneSize);
5509:       DMPlexGetCone(dm, support[s], &cone);
5510:       for (c = 0; c < coneSize; ++c) {
5511:         if (cone[c] == p) break;
5512:       }
5513:       if (c >= coneSize) {
5514:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
5515:         for (c = 0; c < supportSize; ++c) {
5516:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
5517:         }
5518:         PetscPrintf(PETSC_COMM_SELF, "\n");
5519:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
5520:         for (c = 0; c < coneSize; ++c) {
5521:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
5522:         }
5523:         PetscPrintf(PETSC_COMM_SELF, "\n");
5524:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
5525:       }
5526:     }
5527:   }
5528:   PetscSectionGetStorageSize(coneSection, &csize);
5529:   PetscSectionGetStorageSize(supportSection, &ssize);
5530:   if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
5531:   return(0);
5532: }

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

5539:   Input Parameters:
5540: + dm - The DMPlex object
5541: . isSimplex - Are the cells simplices or tensor products
5542: - cellHeight - Normally 0

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

5546:   Level: developer

5548: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
5549: @*/
5550: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5551: {
5552:   PetscInt       dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;

5557:   DMGetDimension(dm, &dim);
5558:   switch (dim) {
5559:   case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
5560:   case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
5561:   case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
5562:   default:
5563:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
5564:   }
5565:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5566:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5567:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5568:   cMax = cMax >= 0 ? cMax : cEnd;
5569:   for (c = cStart; c < cMax; ++c) {
5570:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5572:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5573:     for (cl = 0; cl < closureSize*2; cl += 2) {
5574:       const PetscInt p = closure[cl];
5575:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5576:     }
5577:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5578:     if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has  %d vertices != %d", c, coneSize, numCorners);
5579:   }
5580:   for (c = cMax; c < cEnd; ++c) {
5581:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5583:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5584:     for (cl = 0; cl < closureSize*2; cl += 2) {
5585:       const PetscInt p = closure[cl];
5586:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5587:     }
5588:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5589:     if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has  %d vertices > %d", c, coneSize, numHybridCorners);
5590:   }
5591:   return(0);
5592: }

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

5599:   Input Parameters:
5600: + dm - The DMPlex object
5601: . isSimplex - Are the cells simplices or tensor products
5602: - cellHeight - Normally 0

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

5606:   Level: developer

5608: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
5609: @*/
5610: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5611: {
5612:   PetscInt       pMax[4];
5613:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, h;

5618:   DMGetDimension(dm, &dim);
5619:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5620:   DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
5621:   for (h = cellHeight; h < dim; ++h) {
5622:     DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
5623:     for (c = cStart; c < cEnd; ++c) {
5624:       const PetscInt *cone, *ornt, *faces;
5625:       PetscInt        numFaces, faceSize, coneSize,f;
5626:       PetscInt       *closure = NULL, closureSize, cl, numCorners = 0;

5628:       if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
5629:       DMPlexGetConeSize(dm, c, &coneSize);
5630:       DMPlexGetCone(dm, c, &cone);
5631:       DMPlexGetConeOrientation(dm, c, &ornt);
5632:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5633:       for (cl = 0; cl < closureSize*2; cl += 2) {
5634:         const PetscInt p = closure[cl];
5635:         if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
5636:       }
5637:       DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
5638:       if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
5639:       for (f = 0; f < numFaces; ++f) {
5640:         PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;

5642:         DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
5643:         for (cl = 0; cl < fclosureSize*2; cl += 2) {
5644:           const PetscInt p = fclosure[cl];
5645:           if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
5646:         }
5647:         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);
5648:         for (v = 0; v < fnumCorners; ++v) {
5649:           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]);
5650:         }
5651:         DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
5652:       }
5653:       DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
5654:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5655:     }
5656:   }
5657:   return(0);
5658: }

5662: /* Pointwise interpolation
5663:      Just code FEM for now
5664:      u^f = I u^c
5665:      sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
5666:      u^f_i = sum_j psi^f_i I phi^c_j u^c_j
5667:      I_{ij} = psi^f_i phi^c_j
5668: */
5669: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
5670: {
5671:   PetscSection   gsc, gsf;
5672:   PetscInt       m, n;
5673:   void          *ctx;
5674:   DM             cdm;
5675:   PetscBool      regular;

5679:   DMGetDefaultGlobalSection(dmFine, &gsf);
5680:   PetscSectionGetConstrainedStorageSize(gsf, &m);
5681:   DMGetDefaultGlobalSection(dmCoarse, &gsc);
5682:   PetscSectionGetConstrainedStorageSize(gsc, &n);

5684:   MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
5685:   MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
5686:   MatSetType(*interpolation, dmCoarse->mattype);
5687:   DMGetApplicationContext(dmFine, &ctx);

5689:   DMPlexGetCoarseDM(dmFine, &cdm);
5690:   DMPlexGetRegularRefinement(dmFine, &regular);
5691:   if (regular && cdm == dmCoarse) {DMPlexComputeInterpolatorNested(dmCoarse, dmFine, *interpolation, ctx);}
5692:   else                            {DMPlexComputeInterpolatorGeneral(dmCoarse, dmFine, *interpolation, ctx);}
5693:   MatViewFromOptions(*interpolation, NULL, "-interp_mat_view");
5694:   /* Use naive scaling */
5695:   DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
5696:   return(0);
5697: }

5701: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat)
5702: {
5704:   VecScatter     ctx;

5707:   DMPlexComputeInjectorFEM(dmCoarse, dmFine, &ctx, NULL);
5708:   MatCreateScatter(PetscObjectComm((PetscObject)ctx), ctx, mat);
5709:   VecScatterDestroy(&ctx);
5710:   return(0);
5711: }

5715: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
5716: {
5717:   PetscSection   section;
5718:   IS            *bcPoints, *bcComps;
5719:   PetscBool     *isFE;
5720:   PetscInt      *bcFields, *numComp, *numDof;
5721:   PetscInt       depth, dim, numBd, numBC = 0, numFields, bd, bc = 0, f;
5722:   PetscInt       cStart, cEnd, cEndInterior;

5726:   DMGetNumFields(dm, &numFields);
5727:   /* FE and FV boundary conditions are handled slightly differently */
5728:   PetscMalloc1(numFields, &isFE);
5729:   for (f = 0; f < numFields; ++f) {
5730:     PetscObject  obj;
5731:     PetscClassId id;

5733:     DMGetField(dm, f, &obj);
5734:     PetscObjectGetClassId(obj, &id);
5735:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
5736:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
5737:     else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
5738:   }
5739:   /* Allocate boundary point storage for FEM boundaries */
5740:   DMPlexGetDepth(dm, &depth);
5741:   DMGetDimension(dm, &dim);
5742:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
5743:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
5744:   DMPlexGetNumBoundary(dm, &numBd);
5745:   for (bd = 0; bd < numBd; ++bd) {
5746:     PetscInt  field;
5747:     PetscBool isEssential;

5749:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL, NULL);
5750:     if (isFE[field] && isEssential) ++numBC;
5751:   }
5752:   /* Add ghost cell boundaries for FVM */
5753:   for (f = 0; f < numFields; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
5754:   PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);
5755:   /* Constrain ghost cells for FV */
5756:   for (f = 0; f < numFields; ++f) {
5757:     PetscInt *newidx, c;

5759:     if (isFE[f] || cEndInterior < 0) continue;
5760:     PetscMalloc1(cEnd-cEndInterior,&newidx);
5761:     for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
5762:     bcFields[bc] = f;
5763:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5764:   }
5765:   /* Handle FEM Dirichlet boundaries */
5766:   for (bd = 0; bd < numBd; ++bd) {
5767:     const char     *bdLabel;
5768:     DMLabel         label;
5769:     const PetscInt *comps;
5770:     const PetscInt *values;
5771:     PetscInt        bd2, field, numComps, numValues;
5772:     PetscBool       isEssential, duplicate = PETSC_FALSE;

5774:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
5775:     if (!isFE[field]) continue;
5776:     DMPlexGetLabel(dm, bdLabel, &label);
5777:     /* Only want to modify label once */
5778:     for (bd2 = 0; bd2 < bd; ++bd2) {
5779:       const char *bdname;
5780:       DMPlexGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
5781:       PetscStrcmp(bdname, bdLabel, &duplicate);
5782:       if (duplicate) break;
5783:     }
5784:     if (!duplicate && (isFE[field])) {
5785:       DMPlexLabelComplete(dm, label);
5786:       DMPlexLabelAddCells(dm, label);
5787:     }
5788:     /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
5789:     if (isEssential) {
5790:       PetscInt       *newidx;
5791:       PetscInt        n, newn = 0, p, v;

5793:       bcFields[bc] = field;
5794:       if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
5795:       for (v = 0; v < numValues; ++v) {
5796:         IS              tmp;
5797:         const PetscInt *idx;

5799:         DMPlexGetStratumIS(dm, bdLabel, values[v], &tmp);
5800:         if (!tmp) continue;
5801:         ISGetLocalSize(tmp, &n);
5802:         ISGetIndices(tmp, &idx);
5803:         if (isFE[field]) {
5804:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
5805:         } else {
5806:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
5807:         }
5808:         ISRestoreIndices(tmp, &idx);
5809:         ISDestroy(&tmp);
5810:       }
5811:       PetscMalloc1(newn,&newidx);
5812:       newn = 0;
5813:       for (v = 0; v < numValues; ++v) {
5814:         IS              tmp;
5815:         const PetscInt *idx;

5817:         DMPlexGetStratumIS(dm, bdLabel, values[v], &tmp);
5818:         if (!tmp) continue;
5819:         ISGetLocalSize(tmp, &n);
5820:         ISGetIndices(tmp, &idx);
5821:         if (isFE[field]) {
5822:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
5823:         } else {
5824:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
5825:         }
5826:         ISRestoreIndices(tmp, &idx);
5827:         ISDestroy(&tmp);
5828:       }
5829:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5830:     }
5831:   }
5832:   /* Handle discretization */
5833:   PetscCalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
5834:   for (f = 0; f < numFields; ++f) {
5835:     PetscObject obj;

5837:     DMGetField(dm, f, &obj);
5838:     if (isFE[f]) {
5839:       PetscFE         fe = (PetscFE) obj;
5840:       const PetscInt *numFieldDof;
5841:       PetscInt        d;

5843:       PetscFEGetNumComponents(fe, &numComp[f]);
5844:       PetscFEGetNumDof(fe, &numFieldDof);
5845:       for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
5846:     } else {
5847:       PetscFV fv = (PetscFV) obj;

5849:       PetscFVGetNumComponents(fv, &numComp[f]);
5850:       numDof[f*(dim+1)+dim] = numComp[f];
5851:     }
5852:   }
5853:   for (f = 0; f < numFields; ++f) {
5854:     PetscInt d;
5855:     for (d = 1; d < dim; ++d) {
5856:       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.");
5857:     }
5858:   }
5859:   DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section);
5860:   for (f = 0; f < numFields; ++f) {
5861:     PetscFE     fe;
5862:     const char *name;

5864:     DMGetField(dm, f, (PetscObject *) &fe);
5865:     PetscObjectGetName((PetscObject) fe, &name);
5866:     PetscSectionSetFieldName(section, f, name);
5867:   }
5868:   DMSetDefaultSection(dm, section);
5869:   PetscSectionDestroy(&section);
5870:   for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
5871:   PetscFree3(bcFields,bcPoints,bcComps);
5872:   PetscFree2(numComp,numDof);
5873:   PetscFree(isFE);
5874:   return(0);
5875: }

5879: /*@
5880:   DMPlexGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

5882:   Input Parameter:
5883: . dm - The DMPlex object

5885:   Output Parameter:
5886: . cdm - The coarse DM

5888:   Level: intermediate

5890: .seealso: DMPlexSetCoarseDM()
5891: @*/
5892: PetscErrorCode DMPlexGetCoarseDM(DM dm, DM *cdm)
5893: {
5897:   *cdm = ((DM_Plex *) dm->data)->coarseMesh;
5898:   return(0);
5899: }

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

5906:   Input Parameters:
5907: + dm - The DMPlex object
5908: - cdm - The coarse DM

5910:   Level: intermediate

5912: .seealso: DMPlexGetCoarseDM()
5913: @*/
5914: PetscErrorCode DMPlexSetCoarseDM(DM dm, DM cdm)
5915: {
5916:   DM_Plex       *mesh;

5922:   mesh = (DM_Plex *) dm->data;
5923:   DMDestroy(&mesh->coarseMesh);
5924:   mesh->coarseMesh = cdm;
5925:   PetscObjectReference((PetscObject) mesh->coarseMesh);
5926:   return(0);
5927: }

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

5934:   Input Parameter:
5935: . dm - The DMPlex object

5937:   Output Parameter:
5938: . regular - The flag

5940:   Level: intermediate

5942: .seealso: DMPlexSetRegularRefinement()
5943: @*/
5944: PetscErrorCode DMPlexGetRegularRefinement(DM dm, PetscBool *regular)
5945: {
5949:   *regular = ((DM_Plex *) dm->data)->regularRefinement;
5950:   return(0);
5951: }

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

5958:   Input Parameters:
5959: + dm - The DMPlex object
5960: - regular - The flag

5962:   Level: intermediate

5964: .seealso: DMPlexGetRegularRefinement()
5965: @*/
5966: PetscErrorCode DMPlexSetRegularRefinement(DM dm, PetscBool regular)
5967: {
5970:   ((DM_Plex *) dm->data)->regularRefinement = regular;
5971:   return(0);
5972: }

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

5981:   not collective

5983:   Input Parameters:
5984: . dm - The DMPlex object

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


5991:   Level: intermediate

5993: .seealso: DMPlexSetAnchors(), DMGetConstraints(), DMSetConstraints()
5994: @*/
5995: PetscErrorCode DMPlexGetAnchors(DM dm, PetscSection *anchorSection, IS *anchorIS)
5996: {
5997:   DM_Plex *plex = (DM_Plex *)dm->data;

6002:   if (!plex->anchorSection && !plex->anchorIS && plex->createanchors) {(*plex->createanchors)(dm);}
6003:   if (anchorSection) *anchorSection = plex->anchorSection;
6004:   if (anchorIS) *anchorIS = plex->anchorIS;
6005:   return(0);
6006: }

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

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

6018:   collective on dm

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

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

6027:   Level: intermediate

6029: .seealso: DMPlexGetAnchors(), DMGetConstraints(), DMSetConstraints()
6030: @*/
6031: PetscErrorCode DMPlexSetAnchors(DM dm, PetscSection anchorSection, IS anchorIS)
6032: {
6033:   DM_Plex        *plex = (DM_Plex *)dm->data;
6034:   PetscMPIInt    result;

6039:   if (anchorSection) {
6041:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorSection),&result);
6042:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor section must have local communicator");
6043:   }
6044:   if (anchorIS) {
6046:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorIS),&result);
6047:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor IS must have local communicator");
6048:   }

6050:   PetscObjectReference((PetscObject)anchorSection);
6051:   PetscSectionDestroy(&plex->anchorSection);
6052:   plex->anchorSection = anchorSection;

6054:   PetscObjectReference((PetscObject)anchorIS);
6055:   ISDestroy(&plex->anchorIS);
6056:   plex->anchorIS = anchorIS;

6058: #if defined(PETSC_USE_DEBUG)
6059:   if (anchorIS && anchorSection) {
6060:     PetscInt size, a, pStart, pEnd;
6061:     const PetscInt *anchors;

6063:     PetscSectionGetChart(anchorSection,&pStart,&pEnd);
6064:     ISGetLocalSize(anchorIS,&size);
6065:     ISGetIndices(anchorIS,&anchors);
6066:     for (a = 0; a < size; a++) {
6067:       PetscInt p;

6069:       p = anchors[a];
6070:       if (p >= pStart && p < pEnd) {
6071:         PetscInt dof;

6073:         PetscSectionGetDof(anchorSection,p,&dof);
6074:         if (dof) {
6075:           PetscErrorCode ierr2;

6077:           ierr2 = ISRestoreIndices(anchorIS,&anchors);CHKERRQ(ierr2);
6078:           SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point %d cannot be constrained and an anchor",p);
6079:         }
6080:       }
6081:     }
6082:     ISRestoreIndices(anchorIS,&anchors);
6083:   }
6084: #endif
6085:   /* reset the generic constraints */
6086:   DMSetDefaultConstraints(dm,NULL,NULL);
6087:   return(0);
6088: }

6092: static PetscErrorCode DMPlexCreateConstraintSection_Anchors(DM dm, PetscSection section, PetscSection *cSec)
6093: {
6094:   PetscSection anchorSection;
6095:   PetscInt pStart, pEnd, p, dof, numFields, f;

6100:   DMPlexGetAnchors(dm,&anchorSection,NULL);
6101:   PetscSectionCreate(PETSC_COMM_SELF,cSec);
6102:   PetscSectionGetNumFields(section,&numFields);
6103:   PetscSectionSetNumFields(*cSec,numFields);
6104:   PetscSectionGetChart(anchorSection,&pStart,&pEnd);
6105:   PetscSectionSetChart(*cSec,pStart,pEnd);
6106:   for (p = pStart; p < pEnd; p++) {
6107:     PetscSectionGetDof(anchorSection,p,&dof);
6108:     if (dof) {
6109:       PetscSectionGetDof(section,p,&dof);
6110:       PetscSectionSetDof(*cSec,p,dof);
6111:       for (f = 0; f < numFields; f++) {
6112:         PetscSectionGetFieldDof(section,p,f,&dof);
6113:         PetscSectionSetFieldDof(*cSec,p,f,dof);
6114:       }
6115:     }
6116:   }
6117:   PetscSectionSetUp(*cSec);
6118:   return(0);
6119: }

6123: static PetscErrorCode DMPlexCreateConstraintMatrix_Anchors(DM dm, PetscSection section, PetscSection cSec, Mat *cMat)
6124: {
6125:   PetscSection aSec;
6126:   PetscInt pStart, pEnd, p, dof, aDof, aOff, off, nnz, annz, m, n, q, a, offset, *i, *j;
6127:   const PetscInt *anchors;
6128:   PetscInt numFields, f;
6129:   IS aIS;

6134:   PetscSectionGetStorageSize(cSec, &m);
6135:   PetscSectionGetStorageSize(section, &n);
6136:   MatCreate(PETSC_COMM_SELF,cMat);
6137:   MatSetSizes(*cMat,m,n,m,n);
6138:   MatSetType(*cMat,MATSEQAIJ);
6139:   DMPlexGetAnchors(dm,&aSec,&aIS);
6140:   ISGetIndices(aIS,&anchors);
6141:   PetscSectionGetChart(aSec,&pStart,&pEnd);
6142:   PetscMalloc1(m+1,&i);
6143:   i[0] = 0;
6144:   PetscSectionGetNumFields(section,&numFields);
6145:   for (p = pStart; p < pEnd; p++) {
6146:     PetscSectionGetDof(aSec,p,&dof);
6147:     if (!dof) continue;
6148:     PetscSectionGetOffset(aSec,p,&off);
6149:     if (numFields) {
6150:       for (f = 0; f < numFields; f++) {
6151:         annz = 0;
6152:         for (q = 0; q < dof; q++) {
6153:           a = anchors[off + q];
6154:           PetscSectionGetFieldDof(section,a,f,&aDof);
6155:           annz += aDof;
6156:         }
6157:         PetscSectionGetFieldDof(cSec,p,f,&dof);
6158:         PetscSectionGetFieldOffset(cSec,p,f,&off);
6159:         for (q = 0; q < dof; q++) {
6160:           i[off + q + 1] = i[off + q] + annz;
6161:         }
6162:       }
6163:     }
6164:     else {
6165:       annz = 0;
6166:       for (q = 0; q < dof; q++) {
6167:         a = anchors[off + q];
6168:         PetscSectionGetDof(section,a,&aDof);
6169:         annz += aDof;
6170:       }
6171:       PetscSectionGetDof(cSec,p,&dof);
6172:       PetscSectionGetOffset(cSec,p,&off);
6173:       for (q = 0; q < dof; q++) {
6174:         i[off + q + 1] = i[off + q] + annz;
6175:       }
6176:     }
6177:   }
6178:   nnz = i[m];
6179:   PetscMalloc1(nnz,&j);
6180:   offset = 0;
6181:   for (p = pStart; p < pEnd; p++) {
6182:     if (numFields) {
6183:       for (f = 0; f < numFields; f++) {
6184:         PetscSectionGetFieldDof(cSec,p,f,&dof);
6185:         for (q = 0; q < dof; q++) {
6186:           PetscInt rDof, rOff, r;
6187:           PetscSectionGetDof(aSec,p,&rDof);
6188:           PetscSectionGetOffset(aSec,p,&rOff);
6189:           for (r = 0; r < rDof; r++) {
6190:             PetscInt s;

6192:             a = anchors[rOff + r];

6194:             PetscSectionGetFieldDof(section,a,f,&aDof);
6195:             PetscSectionGetFieldOffset(section,a,f,&aOff);
6196:             for (s = 0; s < aDof; s++) {
6197:               j[offset++] = aOff + s;
6198:             }
6199:           }
6200:         }
6201:       }
6202:     }
6203:     else {
6204:       PetscSectionGetDof(cSec,p,&dof);
6205:       for (q = 0; q < dof; q++) {
6206:         PetscInt rDof, rOff, r;
6207:         PetscSectionGetDof(aSec,p,&rDof);
6208:         PetscSectionGetOffset(aSec,p,&rOff);
6209:         for (r = 0; r < rDof; r++) {
6210:           PetscInt s;

6212:           a = anchors[rOff + r];

6214:           PetscSectionGetDof(section,a,&aDof);
6215:           PetscSectionGetOffset(section,a,&aOff);
6216:           for (s = 0; s < aDof; s++) {
6217:             j[offset++] = aOff + s;
6218:           }
6219:         }
6220:       }
6221:     }
6222:   }
6223:   MatSeqAIJSetPreallocationCSR(*cMat,i,j,NULL);
6224:   PetscFree(i);
6225:   PetscFree(j);
6226:   ISRestoreIndices(aIS,&anchors);
6227:   return(0);
6228: }

6232: PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm)
6233: {
6234:   DM_Plex        *plex = (DM_Plex *)dm->data;
6235:   PetscSection   anchorSection, section, cSec;
6236:   Mat            cMat;

6241:   DMPlexGetAnchors(dm,&anchorSection,NULL);
6242:   if (anchorSection) {
6243:     PetscDS  ds;
6244:     PetscInt nf;

6246:     DMGetDefaultSection(dm,&section);
6247:     DMPlexCreateConstraintSection_Anchors(dm,section,&cSec);
6248:     DMPlexCreateConstraintMatrix_Anchors(dm,section,cSec,&cMat);
6249:     DMGetDS(dm,&ds);
6250:     PetscDSGetNumFields(ds,&nf);
6251:     if (nf && plex->computeanchormatrix) {(*plex->computeanchormatrix)(dm,section,cSec,cMat);}
6252:     DMSetDefaultConstraints(dm,cSec,cMat);
6253:     PetscSectionDestroy(&cSec);
6254:     MatDestroy(&cMat);
6255:   }
6256:   return(0);
6257: }