Actual source code: plex.c

petsc-master 2015-03-02
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_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM, DMPLEX_InterpolatorFEM, DMPLEX_InjectorFEM;

  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, 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:   PetscSectionGetChart(section, &pStart, &pEnd);
 26:   if (field >= 0) {
 27:     if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetFieldDof(section, vStart, field, &vdof);}
 28:     if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetFieldDof(section, cStart, field, &cdof);}
 29:   } else {
 30:     if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetDof(section, vStart, &vdof);}
 31:     if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetDof(section, cStart, &cdof);}
 32:   }
 33:   if (vdof) {
 34:     *sStart = vStart;
 35:     *sEnd   = vEnd;
 36:     if (vdof == dim) *ft = PETSC_VTK_POINT_VECTOR_FIELD;
 37:     else             *ft = PETSC_VTK_POINT_FIELD;
 38:   } else if (cdof) {
 39:     *sStart = cStart;
 40:     *sEnd   = cEnd;
 41:     if (cdof == dim) *ft = PETSC_VTK_CELL_VECTOR_FIELD;
 42:     else             *ft = PETSC_VTK_CELL_FIELD;
 43:   } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Could not classify input Vec for VTK");
 44:   return(0);
 45: }

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

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

 65:     DMGetNumFields(dm, &numFields);
 66:     if (numFields) {
 67:       PetscObject fe;

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

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

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

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

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

137: PetscErrorCode VecLoad_Plex_Local(Vec v, PetscViewer viewer)
138: {
139:   DM             dm;
140:   PetscBool      ishdf5;

144:   VecGetDM(v, &dm);
145:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
146:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
147:   if (ishdf5) {
148:     DM          dmBC;
149:     Vec         gv;
150:     const char *name;

152:     DMGetOutputDM(dm, &dmBC);
153:     DMGetGlobalVector(dmBC, &gv);
154:     PetscObjectGetName((PetscObject) v, &name);
155:     PetscObjectSetName((PetscObject) gv, name);
156:     VecLoad_Default(gv, viewer);
157:     DMGlobalToLocalBegin(dmBC, gv, INSERT_VALUES, v);
158:     DMGlobalToLocalEnd(dmBC, gv, INSERT_VALUES, v);
159:     DMRestoreGlobalVector(dmBC, &gv);
160:   } else {
161:     VecLoad_Default(v, viewer);
162:   }
163:   return(0);
164: }

168: PetscErrorCode VecLoad_Plex(Vec v, PetscViewer viewer)
169: {
170:   DM             dm;
171:   PetscBool      ishdf5;

175:   VecGetDM(v, &dm);
176:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
177:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
178:   if (ishdf5) {
179: #if defined(PETSC_HAVE_HDF5)
180:     VecLoad_Plex_HDF5(v, viewer);
181: #else
182:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
183: #endif
184:   } else {
185:     VecLoad_Default(v, viewer);
186:   }
187:   return(0);
188: }

192: PetscErrorCode DMPlexView_Ascii_Geometry(DM dm, PetscViewer viewer)
193: {
194:   PetscSection       coordSection;
195:   Vec                coordinates;
196:   DMLabel            depthLabel;
197:   const char        *name[4];
198:   const PetscScalar *a;
199:   PetscInt           dim, pStart, pEnd, cStart, cEnd, c;
200:   PetscErrorCode     ierr;

203:   DMGetDimension(dm, &dim);
204:   DMGetCoordinatesLocal(dm, &coordinates);
205:   DMGetCoordinateSection(dm, &coordSection);
206:   DMPlexGetDepthLabel(dm, &depthLabel);
207:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
208:   PetscSectionGetChart(coordSection, &pStart, &pEnd);
209:   VecGetArrayRead(coordinates, &a);
210:   name[0]     = "vertex";
211:   name[1]     = "edge";
212:   name[dim-1] = "face";
213:   name[dim]   = "cell";
214:   for (c = cStart; c < cEnd; ++c) {
215:     PetscInt *closure = NULL;
216:     PetscInt  closureSize, cl;

218:     PetscViewerASCIIPrintf(viewer, "Geometry for cell %d:\n", c);
219:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
220:     PetscViewerASCIIPushTab(viewer);
221:     for (cl = 0; cl < closureSize*2; cl += 2) {
222:       PetscInt point = closure[cl], depth, dof, off, d, p;

224:       if ((point < pStart) || (point >= pEnd)) continue;
225:       PetscSectionGetDof(coordSection, point, &dof);
226:       if (!dof) continue;
227:       DMLabelGetValue(depthLabel, point, &depth);
228:       PetscSectionGetOffset(coordSection, point, &off);
229:       PetscViewerASCIIPrintf(viewer, "%s %d coords:", name[depth], point);
230:       for (p = 0; p < dof/dim; ++p) {
231:         PetscViewerASCIIPrintf(viewer, " (");
232:         for (d = 0; d < dim; ++d) {
233:           if (d > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
234:           PetscViewerASCIIPrintf(viewer, "%g", PetscRealPart(a[off+p*dim+d]));
235:         }
236:         PetscViewerASCIIPrintf(viewer, ")");
237:       }
238:       PetscViewerASCIIPrintf(viewer, "\n");
239:     }
240:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
241:     PetscViewerASCIIPopTab(viewer);
242:   }
243:   VecRestoreArrayRead(coordinates, &a);
244:   return(0);
245: }

249: PetscErrorCode DMPlexView_Ascii(DM dm, PetscViewer viewer)
250: {
251:   DM_Plex          *mesh = (DM_Plex*) dm->data;
252:   DM                cdm;
253:   DMLabel           markers;
254:   PetscSection      coordSection;
255:   Vec               coordinates;
256:   PetscViewerFormat format;
257:   PetscErrorCode    ierr;

260:   DMGetCoordinateDM(dm, &cdm);
261:   DMGetDefaultSection(cdm, &coordSection);
262:   DMGetCoordinatesLocal(dm, &coordinates);
263:   PetscViewerGetFormat(viewer, &format);
264:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
265:     const char *name;
266:     PetscInt    maxConeSize, maxSupportSize;
267:     PetscInt    pStart, pEnd, p;
268:     PetscMPIInt rank, size;

270:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
271:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
272:     PetscObjectGetName((PetscObject) dm, &name);
273:     DMPlexGetChart(dm, &pStart, &pEnd);
274:     DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
275:     PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
276:     PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
277:     PetscViewerASCIISynchronizedPrintf(viewer, "Max sizes cone: %D support: %D\n", maxConeSize, maxSupportSize);
278:     PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
279:     PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
280:     for (p = pStart; p < pEnd; ++p) {
281:       PetscInt dof, off, s;

283:       PetscSectionGetDof(mesh->supportSection, p, &dof);
284:       PetscSectionGetOffset(mesh->supportSection, p, &off);
285:       for (s = off; s < off+dof; ++s) {
286:         PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D ----> %D\n", rank, p, mesh->supports[s]);
287:       }
288:     }
289:     PetscViewerFlush(viewer);
290:     PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
291:     for (p = pStart; p < pEnd; ++p) {
292:       PetscInt dof, off, c;

294:       PetscSectionGetDof(mesh->coneSection, p, &dof);
295:       PetscSectionGetOffset(mesh->coneSection, p, &off);
296:       for (c = off; c < off+dof; ++c) {
297:         PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D <---- %D (%D)\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]);
298:       }
299:     }
300:     PetscViewerFlush(viewer);
301:     PetscSectionGetChart(coordSection, &pStart, NULL);
302:     if (pStart >= 0) {PetscSectionVecView(coordSection, coordinates, viewer);}
303:     DMPlexGetLabel(dm, "marker", &markers);
304:     DMLabelView(markers,viewer);
305:     if (size > 1) {
306:       PetscSF sf;

308:       DMGetPointSF(dm, &sf);
309:       PetscSFView(sf, viewer);
310:     }
311:     PetscViewerFlush(viewer);
312:   } else if (format == PETSC_VIEWER_ASCII_LATEX) {
313:     const char  *name, *color;
314:     const char  *defcolors[3]  = {"gray", "orange", "green"};
315:     const char  *deflcolors[4] = {"blue", "cyan", "red", "magenta"};
316:     PetscReal    scale         = 2.0;
317:     PetscBool    useNumbers    = PETSC_TRUE, useLabels, useColors;
318:     double       tcoords[3];
319:     PetscScalar *coords;
320:     PetscInt     numLabels, l, numColors, numLColors, dim, depth, cStart, cEnd, c, vStart, vEnd, v, eStart = 0, eEnd = 0, e, p;
321:     PetscMPIInt  rank, size;
322:     char       **names, **colors, **lcolors;

324:     DMGetDimension(dm, &dim);
325:     DMPlexGetDepth(dm, &depth);
326:     DMPlexGetNumLabels(dm, &numLabels);
327:     numLabels  = PetscMax(numLabels, 10);
328:     numColors  = 10;
329:     numLColors = 10;
330:     PetscCalloc3(numLabels, &names, numColors, &colors, numLColors, &lcolors);
331:     PetscOptionsGetReal(((PetscObject) viewer)->prefix, "-dm_plex_view_scale", &scale, NULL);
332:     PetscOptionsGetBool(((PetscObject) viewer)->prefix, "-dm_plex_view_numbers", &useNumbers, NULL);
333:     PetscOptionsGetStringArray(((PetscObject) viewer)->prefix, "-dm_plex_view_labels", names, &numLabels, &useLabels);
334:     if (!useLabels) numLabels = 0;
335:     PetscOptionsGetStringArray(((PetscObject) viewer)->prefix, "-dm_plex_view_colors", colors, &numColors, &useColors);
336:     if (!useColors) {
337:       numColors = 3;
338:       for (c = 0; c < numColors; ++c) {PetscStrallocpy(defcolors[c], &colors[c]);}
339:     }
340:     PetscOptionsGetStringArray(((PetscObject) viewer)->prefix, "-dm_plex_view_lcolors", lcolors, &numLColors, &useColors);
341:     if (!useColors) {
342:       numLColors = 4;
343:       for (c = 0; c < numLColors; ++c) {PetscStrallocpy(deflcolors[c], &lcolors[c]);}
344:     }
345:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
346:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
347:     PetscObjectGetName((PetscObject) dm, &name);
348:     PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
349:     PetscViewerASCIIPrintf(viewer, "\
350: \\documentclass[tikz]{standalone}\n\n\
351: \\usepackage{pgflibraryshapes}\n\
352: \\usetikzlibrary{backgrounds}\n\
353: \\usetikzlibrary{arrows}\n\
354: \\begin{document}\n");
355:     if (size > 1) {
356:       PetscViewerASCIIPrintf(viewer, "%s for process ", name);
357:       for (p = 0; p < size; ++p) {
358:         if (p > 0 && p == size-1) {
359:           PetscViewerASCIIPrintf(viewer, ", and ", colors[p%numColors], p);
360:         } else if (p > 0) {
361:           PetscViewerASCIIPrintf(viewer, ", ", colors[p%numColors], p);
362:         }
363:         PetscViewerASCIIPrintf(viewer, "{\\textcolor{%s}%D}", colors[p%numColors], p);
364:       }
365:       PetscViewerASCIIPrintf(viewer, ".\n\n\n");
366:     }
367:     PetscViewerASCIIPrintf(viewer, "\\begin{tikzpicture}[scale = %g,font=\\fontsize{8}{8}\\selectfont]\n", 1.0);
368:     /* Plot vertices */
369:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
370:     VecGetArray(coordinates, &coords);
371:     for (v = vStart; v < vEnd; ++v) {
372:       PetscInt  off, dof, d;
373:       PetscBool isLabeled = PETSC_FALSE;

375:       PetscSectionGetDof(coordSection, v, &dof);
376:       PetscSectionGetOffset(coordSection, v, &off);
377:       PetscViewerASCIISynchronizedPrintf(viewer, "\\path (");
378:       for (d = 0; d < dof; ++d) {
379:         tcoords[d] = (double) (scale*PetscRealPart(coords[off+d]));
380:         tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
381:       }
382:       /* Rotate coordinates since PGF makes z point out of the page instead of up */
383:       if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
384:       for (d = 0; d < dof; ++d) {
385:         if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
386:         PetscViewerASCIISynchronizedPrintf(viewer, "%g", tcoords[d]);
387:       }
388:       color = colors[rank%numColors];
389:       for (l = 0; l < numLabels; ++l) {
390:         PetscInt val;
391:         DMPlexGetLabelValue(dm, names[l], v, &val);
392:         if (val >= 0) {color = lcolors[l%numLColors]; isLabeled = PETSC_TRUE; break;}
393:       }
394:       if (useNumbers) {
395:         PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [draw,shape=circle,color=%s] {%D};\n", v, rank, color, v);
396:       } else {
397:         PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [fill,inner sep=%dpt,shape=circle,color=%s] {};\n", v, rank, !isLabeled ? 1 : 2, color);
398:       }
399:     }
400:     VecRestoreArray(coordinates, &coords);
401:     PetscViewerFlush(viewer);
402:     /* Plot edges */
403:     if (depth > 1) {DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);}
404:     if (dim < 3 && useNumbers) {
405:       VecGetArray(coordinates, &coords);
406:       PetscViewerASCIIPrintf(viewer, "\\path\n");
407:       for (e = eStart; e < eEnd; ++e) {
408:         const PetscInt *cone;
409:         PetscInt        coneSize, offA, offB, dof, d;

411:         DMPlexGetConeSize(dm, e, &coneSize);
412:         if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %d cone should have two vertices, not %d", e, coneSize);
413:         DMPlexGetCone(dm, e, &cone);
414:         PetscSectionGetDof(coordSection, cone[0], &dof);
415:         PetscSectionGetOffset(coordSection, cone[0], &offA);
416:         PetscSectionGetOffset(coordSection, cone[1], &offB);
417:         PetscViewerASCIISynchronizedPrintf(viewer, "(");
418:         for (d = 0; d < dof; ++d) {
419:           tcoords[d] = (double) (scale*PetscRealPart(coords[offA+d]+coords[offB+d]));
420:           tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
421:         }
422:         /* Rotate coordinates since PGF makes z point out of the page instead of up */
423:         if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
424:         for (d = 0; d < dof; ++d) {
425:           if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
426:           PetscViewerASCIISynchronizedPrintf(viewer, "%g", tcoords[d]);
427:         }
428:         color = colors[rank%numColors];
429:         for (l = 0; l < numLabels; ++l) {
430:           PetscInt val;
431:           DMPlexGetLabelValue(dm, names[l], v, &val);
432:           if (val >= 0) {color = lcolors[l%numLColors]; break;}
433:         }
434:         PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [draw,shape=circle,color=%s] {%D} --\n", e, rank, color, e);
435:       }
436:       VecRestoreArray(coordinates, &coords);
437:       PetscViewerFlush(viewer);
438:       PetscViewerASCIIPrintf(viewer, "(0,0);\n");
439:     }
440:     /* Plot cells */
441:     if (dim == 3 || !useNumbers) {
442:       for (e = eStart; e < eEnd; ++e) {
443:         const PetscInt *cone;

445:         color = colors[rank%numColors];
446:         for (l = 0; l < numLabels; ++l) {
447:           PetscInt val;
448:           DMPlexGetLabelValue(dm, names[l], e, &val);
449:           if (val >= 0) {color = lcolors[l%numLColors]; break;}
450:         }
451:         DMPlexGetCone(dm, e, &cone);
452:         PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] (%D_%D) -- (%D_%D);\n", color, cone[0], rank, cone[1], rank);
453:       }
454:     } else {
455:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
456:       for (c = cStart; c < cEnd; ++c) {
457:         PetscInt *closure = NULL;
458:         PetscInt  closureSize, firstPoint = -1;

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

465:           if ((point < vStart) || (point >= vEnd)) continue;
466:           if (firstPoint >= 0) {PetscViewerASCIISynchronizedPrintf(viewer, " -- ");}
467:           PetscViewerASCIISynchronizedPrintf(viewer, "(%D_%D)", point, rank);
468:           if (firstPoint < 0) firstPoint = point;
469:         }
470:         /* Why doesn't this work? PetscViewerASCIISynchronizedPrintf(viewer, " -- cycle;\n"); */
471:         PetscViewerASCIISynchronizedPrintf(viewer, " -- (%D_%D);\n", firstPoint, rank);
472:         DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
473:       }
474:     }
475:     PetscViewerFlush(viewer);
476:     PetscViewerASCIIPrintf(viewer, "\\end{tikzpicture}\n");
477:     PetscViewerASCIIPrintf(viewer, "\\end{document}\n", name);
478:     for (l = 0; l < numLabels;  ++l) {PetscFree(names[l]);}
479:     for (c = 0; c < numColors;  ++c) {PetscFree(colors[c]);}
480:     for (c = 0; c < numLColors; ++c) {PetscFree(lcolors[c]);}
481:     PetscFree3(names, colors, lcolors);
482:   } else {
483:     MPI_Comm    comm;
484:     PetscInt   *sizes, *hybsizes;
485:     PetscInt    locDepth, depth, dim, d, pMax[4];
486:     PetscInt    pStart, pEnd, p;
487:     PetscInt    numLabels, l;
488:     const char *name;
489:     PetscMPIInt size;

491:     PetscObjectGetComm((PetscObject)dm,&comm);
492:     MPI_Comm_size(comm, &size);
493:     DMGetDimension(dm, &dim);
494:     PetscObjectGetName((PetscObject) dm, &name);
495:     if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dim);}
496:     else      {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dim);}
497:     DMPlexGetDepth(dm, &locDepth);
498:     MPI_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
499:     DMPlexGetHybridBounds(dm, &pMax[depth], depth > 0 ? &pMax[depth-1] : NULL, &pMax[1], &pMax[0]);
500:     PetscMalloc2(size,&sizes,size,&hybsizes);
501:     if (depth == 1) {
502:       DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd);
503:       pEnd = pEnd - pStart;
504:       MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
505:       PetscViewerASCIIPrintf(viewer, "  %D-cells:", 0);
506:       for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
507:       PetscViewerASCIIPrintf(viewer, "\n");
508:       DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd);
509:       pEnd = pEnd - pStart;
510:       MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
511:       PetscViewerASCIIPrintf(viewer, "  %D-cells:", dim);
512:       for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
513:       PetscViewerASCIIPrintf(viewer, "\n");
514:     } else {
515:       for (d = 0; d <= dim; d++) {
516:         DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
517:         pEnd    -= pStart;
518:         pMax[d] -= pStart;
519:         MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
520:         MPI_Gather(&pMax[d], 1, MPIU_INT, hybsizes, 1, MPIU_INT, 0, comm);
521:         PetscViewerASCIIPrintf(viewer, "  %D-cells:", d);
522:         for (p = 0; p < size; ++p) {
523:           if (hybsizes[p] >= 0) {PetscViewerASCIIPrintf(viewer, " %D (%D)", sizes[p], sizes[p] - hybsizes[p]);}
524:           else                  {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
525:         }
526:         PetscViewerASCIIPrintf(viewer, "\n");
527:       }
528:     }
529:     PetscFree2(sizes,hybsizes);
530:     DMPlexGetNumLabels(dm, &numLabels);
531:     if (numLabels) {PetscViewerASCIIPrintf(viewer, "Labels:\n");}
532:     for (l = 0; l < numLabels; ++l) {
533:       DMLabel         label;
534:       const char     *name;
535:       IS              valueIS;
536:       const PetscInt *values;
537:       PetscInt        numValues, v;

539:       DMPlexGetLabelName(dm, l, &name);
540:       DMPlexGetLabel(dm, name, &label);
541:       DMLabelGetNumValues(label, &numValues);
542:       PetscViewerASCIIPrintf(viewer, "  %s: %d strata of sizes (", name, numValues);
543:       DMLabelGetValueIS(label, &valueIS);
544:       ISGetIndices(valueIS, &values);
545:       for (v = 0; v < numValues; ++v) {
546:         PetscInt size;

548:         DMLabelGetStratumSize(label, values[v], &size);
549:         if (v > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
550:         PetscViewerASCIIPrintf(viewer, "%d", size);
551:       }
552:       PetscViewerASCIIPrintf(viewer, ")\n");
553:       ISRestoreIndices(valueIS, &values);
554:       ISDestroy(&valueIS);
555:     }
556:   }
557:   return(0);
558: }

562: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
563: {
564:   PetscBool      iascii, ishdf5, isvtk;

570:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
571:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,   &isvtk);
572:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,  &ishdf5);
573:   if (iascii) {
574:     DMPlexView_Ascii(dm, viewer);
575:   } else if (ishdf5) {
576: #if defined(PETSC_HAVE_HDF5)
577:     PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
578:     DMPlexView_HDF5(dm, viewer);
579:     PetscViewerPopFormat(viewer);
580: #else
581:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
582: #endif
583:   }
584:   else if (isvtk) {
585:     DMPlexVTKWriteAll((PetscObject) dm,viewer);
586:   }
587:   return(0);
588: }

592: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
593: {
594:   PetscBool      isbinary, ishdf5;

600:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
601:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,   &ishdf5);
602:   if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");}
603:   else if (ishdf5) {
604: #if defined(PETSC_HAVE_HDF5)
605:     DMPlexLoad_HDF5(dm, viewer);
606: #else
607:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
608: #endif
609:   }
610:   return(0);
611: }

615: static PetscErrorCode BoundaryDestroy(DMBoundary *boundary)
616: {
617:   DMBoundary     b, next;

621:   if (!boundary) return(0);
622:   b = *boundary;
623:   *boundary = NULL;
624:   for (; b; b = next) {
625:     next = b->next;
626:     PetscFree(b->ids);
627:     PetscFree(b->name);
628:     PetscFree(b->labelname);
629:     PetscFree(b);
630:   }
631:   return(0);
632: }

636: PetscErrorCode DMDestroy_Plex(DM dm)
637: {
638:   DM_Plex       *mesh = (DM_Plex*) dm->data;
639:   PlexLabel      next = mesh->labels;

643:   if (--mesh->refct > 0) return(0);
644:   PetscSectionDestroy(&mesh->coneSection);
645:   PetscFree(mesh->cones);
646:   PetscFree(mesh->coneOrientations);
647:   PetscSectionDestroy(&mesh->supportSection);
648:   PetscFree(mesh->supports);
649:   PetscFree(mesh->facesTmp);
650:   PetscFree(mesh->tetgenOpts);
651:   PetscFree(mesh->triangleOpts);
652:   PetscPartitionerDestroy(&mesh->partitioner);
653:   while (next) {
654:     PlexLabel tmp = next->next;

656:     DMLabelDestroy(&next->label);
657:     PetscFree(next);
658:     next = tmp;
659:   }
660:   DMDestroy(&mesh->coarseMesh);
661:   DMLabelDestroy(&mesh->subpointMap);
662:   ISDestroy(&mesh->globalVertexNumbers);
663:   ISDestroy(&mesh->globalCellNumbers);
664:   BoundaryDestroy(&mesh->boundary);
665:   PetscSectionDestroy(&mesh->anchorSection);
666:   ISDestroy(&mesh->anchorIS);
667:   PetscSectionDestroy(&mesh->parentSection);
668:   PetscFree(mesh->parents);
669:   PetscFree(mesh->childIDs);
670:   PetscSectionDestroy(&mesh->childSection);
671:   PetscFree(mesh->children);
672:   DMDestroy(&mesh->referenceTree);
673:   /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
674:   PetscFree(mesh);
675:   return(0);
676: }

680: PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
681: {
682:   PetscSection   section, sectionGlobal;
683:   PetscInt       bs = -1;
684:   PetscInt       localSize;
685:   PetscBool      isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock;
687:   MatType        mtype;
688:   ISLocalToGlobalMapping ltog;

691:   MatInitializePackage();
692:   mtype = dm->mattype;
693:   DMGetDefaultSection(dm, &section);
694:   DMGetDefaultGlobalSection(dm, &sectionGlobal);
695:   /* PetscSectionGetStorageSize(sectionGlobal, &localSize); */
696:   PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
697:   MatCreate(PetscObjectComm((PetscObject)dm), J);
698:   MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
699:   MatSetType(*J, mtype);
700:   MatSetFromOptions(*J);
701:   PetscStrcmp(mtype, MATSHELL, &isShell);
702:   PetscStrcmp(mtype, MATBAIJ, &isBlock);
703:   PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
704:   PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
705:   PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
706:   PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
707:   PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
708:   if (!isShell) {
709:     PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
710:     PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin;

712:     if (bs < 0) {
713:       if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
714:         PetscInt pStart, pEnd, p, dof, cdof;

716:         PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
717:         for (p = pStart; p < pEnd; ++p) {
718:           PetscSectionGetDof(sectionGlobal, p, &dof);
719:           PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
720:           if (dof-cdof) {
721:             if (bs < 0) {
722:               bs = dof-cdof;
723:             } else if (bs != dof-cdof) {
724:               /* Layout does not admit a pointwise block size */
725:               bs = 1;
726:               break;
727:             }
728:           }
729:         }
730:         /* Must have same blocksize on all procs (some might have no points) */
731:         bsLocal = bs;
732:         MPI_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
733:         bsLocal = bs < 0 ? bsMax : bs;
734:         MPI_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
735:         if (bsMin != bsMax) {
736:           bs = 1;
737:         } else {
738:           bs = bsMax;
739:         }
740:       } else {
741:         bs = 1;
742:       }
743:     }
744:     PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
745:     DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix);
746:     PetscFree4(dnz, onz, dnzu, onzu);

748:     /* Set localtoglobalmapping on the matrix for MatSetValuesLocal() to work */
749:     DMGetLocalToGlobalMapping(dm,&ltog);
750:     MatSetLocalToGlobalMapping(*J,ltog,ltog);
751:   }
752:   return(0);
753: }

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

760:   Not collective

762:   Input Parameter:
763: . mesh - The DMPlex

765:   Output Parameters:
766: + pStart - The first mesh point
767: - pEnd   - The upper bound for mesh points

769:   Level: beginner

771: .seealso: DMPlexCreate(), DMPlexSetChart()
772: @*/
773: PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
774: {
775:   DM_Plex       *mesh = (DM_Plex*) dm->data;

780:   PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
781:   return(0);
782: }

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

789:   Not collective

791:   Input Parameters:
792: + mesh - The DMPlex
793: . pStart - The first mesh point
794: - pEnd   - The upper bound for mesh points

796:   Output Parameters:

798:   Level: beginner

800: .seealso: DMPlexCreate(), DMPlexGetChart()
801: @*/
802: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
803: {
804:   DM_Plex       *mesh = (DM_Plex*) dm->data;

809:   PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
810:   PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
811:   return(0);
812: }

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

819:   Not collective

821:   Input Parameters:
822: + mesh - The DMPlex
823: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

825:   Output Parameter:
826: . size - The cone size for point p

828:   Level: beginner

830: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
831: @*/
832: PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size)
833: {
834:   DM_Plex       *mesh = (DM_Plex*) dm->data;

840:   PetscSectionGetDof(mesh->coneSection, p, size);
841:   return(0);
842: }

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

849:   Not collective

851:   Input Parameters:
852: + mesh - The DMPlex
853: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
854: - size - The cone size for point p

856:   Output Parameter:

858:   Note:
859:   This should be called after DMPlexSetChart().

861:   Level: beginner

863: .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart()
864: @*/
865: PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size)
866: {
867:   DM_Plex       *mesh = (DM_Plex*) dm->data;

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

874:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
875:   return(0);
876: }

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

883:   Not collective

885:   Input Parameters:
886: + mesh - The DMPlex
887: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
888: - size - The additional cone size for point p

890:   Output Parameter:

892:   Note:
893:   This should be called after DMPlexSetChart().

895:   Level: beginner

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

907:   PetscSectionAddDof(mesh->coneSection, p, size);
908:   PetscSectionGetDof(mesh->coneSection, p, &csize);

910:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
911:   return(0);
912: }

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

919:   Not collective

921:   Input Parameters:
922: + mesh - The DMPlex
923: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

928:   Level: beginner

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

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

936: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart()
937: @*/
938: PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
939: {
940:   DM_Plex       *mesh = (DM_Plex*) dm->data;
941:   PetscInt       off;

947:   PetscSectionGetOffset(mesh->coneSection, p, &off);
948:   *cone = &mesh->cones[off];
949:   return(0);
950: }

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

957:   Not collective

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

964:   Output Parameter:

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

969:   Level: beginner

971: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
972: @*/
973: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
974: {
975:   DM_Plex       *mesh = (DM_Plex*) dm->data;
976:   PetscInt       pStart, pEnd;
977:   PetscInt       dof, off, c;

982:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
983:   PetscSectionGetDof(mesh->coneSection, p, &dof);
985:   PetscSectionGetOffset(mesh->coneSection, p, &off);
986:   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);
987:   for (c = 0; c < dof; ++c) {
988:     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);
989:     mesh->cones[off+c] = cone[c];
990:   }
991:   return(0);
992: }

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

999:   Not collective

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

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

1011:   Level: beginner

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

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

1019: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
1020: @*/
1021: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
1022: {
1023:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1024:   PetscInt       off;

1029: #if defined(PETSC_USE_DEBUG)
1030:   {
1031:     PetscInt dof;
1032:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1034:   }
1035: #endif
1036:   PetscSectionGetOffset(mesh->coneSection, p, &off);

1038:   *coneOrientation = &mesh->coneOrientations[off];
1039:   return(0);
1040: }

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

1047:   Not collective

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

1057:   Output Parameter:

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

1062:   Level: beginner

1064: .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1065: @*/
1066: PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
1067: {
1068:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1069:   PetscInt       pStart, pEnd;
1070:   PetscInt       dof, off, c;

1075:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1076:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1078:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1079:   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);
1080:   for (c = 0; c < dof; ++c) {
1081:     PetscInt cdof, o = coneOrientation[c];

1083:     PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1084:     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);
1085:     mesh->coneOrientations[off+c] = o;
1086:   }
1087:   return(0);
1088: }

1092: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
1093: {
1094:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1095:   PetscInt       pStart, pEnd;
1096:   PetscInt       dof, off;

1101:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1102:   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);
1103:   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);
1104:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1105:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1106:   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);
1107:   mesh->cones[off+conePos] = conePoint;
1108:   return(0);
1109: }

1113: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1114: {
1115:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1116:   PetscInt       pStart, pEnd;
1117:   PetscInt       dof, off;

1122:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1123:   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);
1124:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1125:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1126:   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);
1127:   mesh->coneOrientations[off+conePos] = coneOrientation;
1128:   return(0);
1129: }

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

1136:   Not collective

1138:   Input Parameters:
1139: + mesh - The DMPlex
1140: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1142:   Output Parameter:
1143: . size - The support size for point p

1145:   Level: beginner

1147: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1148: @*/
1149: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1150: {
1151:   DM_Plex       *mesh = (DM_Plex*) dm->data;

1157:   PetscSectionGetDof(mesh->supportSection, p, size);
1158:   return(0);
1159: }

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

1166:   Not collective

1168:   Input Parameters:
1169: + mesh - The DMPlex
1170: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1171: - size - The support size for point p

1173:   Output Parameter:

1175:   Note:
1176:   This should be called after DMPlexSetChart().

1178:   Level: beginner

1180: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1181: @*/
1182: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1183: {
1184:   DM_Plex       *mesh = (DM_Plex*) dm->data;

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

1191:   mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1192:   return(0);
1193: }

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

1200:   Not collective

1202:   Input Parameters:
1203: + mesh - The DMPlex
1204: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1209:   Level: beginner

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

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

1217: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1218: @*/
1219: PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1220: {
1221:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1222:   PetscInt       off;

1228:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1229:   *support = &mesh->supports[off];
1230:   return(0);
1231: }

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

1238:   Not collective

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

1245:   Output Parameter:

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

1250:   Level: beginner

1252: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1253: @*/
1254: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1255: {
1256:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1257:   PetscInt       pStart, pEnd;
1258:   PetscInt       dof, off, c;

1263:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1264:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1266:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1267:   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);
1268:   for (c = 0; c < dof; ++c) {
1269:     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);
1270:     mesh->supports[off+c] = support[c];
1271:   }
1272:   return(0);
1273: }

1277: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1278: {
1279:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1280:   PetscInt       pStart, pEnd;
1281:   PetscInt       dof, off;

1286:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1287:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1288:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1289:   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);
1290:   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);
1291:   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);
1292:   mesh->supports[off+supportPos] = supportPoint;
1293:   return(0);
1294: }

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

1301:   Not collective

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

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

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

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

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

1322:   Level: beginner

1324: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1325: @*/
1326: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1327: {
1328:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1329:   PetscInt       *closure, *fifo;
1330:   const PetscInt *tmp = NULL, *tmpO = NULL;
1331:   PetscInt        tmpSize, t;
1332:   PetscInt        depth       = 0, maxSize;
1333:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1334:   PetscErrorCode  ierr;

1338:   DMPlexGetDepth(dm, &depth);
1339:   /* This is only 1-level */
1340:   if (useCone) {
1341:     DMPlexGetConeSize(dm, p, &tmpSize);
1342:     DMPlexGetCone(dm, p, &tmp);
1343:     DMPlexGetConeOrientation(dm, p, &tmpO);
1344:   } else {
1345:     DMPlexGetSupportSize(dm, p, &tmpSize);
1346:     DMPlexGetSupport(dm, p, &tmp);
1347:   }
1348:   if (depth == 1) {
1349:     if (*points) {
1350:       closure = *points;
1351:     } else {
1352:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1353:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1354:     }
1355:     closure[0] = p; closure[1] = 0;
1356:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1357:       closure[closureSize]   = tmp[t];
1358:       closure[closureSize+1] = tmpO ? tmpO[t] : 0;
1359:     }
1360:     if (numPoints) *numPoints = closureSize/2;
1361:     if (points)    *points    = closure;
1362:     return(0);
1363:   }
1364:   {
1365:     PetscInt c, coneSeries, s,supportSeries;

1367:     c = mesh->maxConeSize;
1368:     coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1369:     s = mesh->maxSupportSize;
1370:     supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1371:     maxSize = 2*PetscMax(coneSeries,supportSeries);
1372:   }
1373:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1374:   if (*points) {
1375:     closure = *points;
1376:   } else {
1377:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1378:   }
1379:   closure[0] = p; closure[1] = 0;
1380:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1381:     const PetscInt cp = tmp[t];
1382:     const PetscInt co = tmpO ? tmpO[t] : 0;

1384:     closure[closureSize]   = cp;
1385:     closure[closureSize+1] = co;
1386:     fifo[fifoSize]         = cp;
1387:     fifo[fifoSize+1]       = co;
1388:   }
1389:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1390:   while (fifoSize - fifoStart) {
1391:     const PetscInt q   = fifo[fifoStart];
1392:     const PetscInt o   = fifo[fifoStart+1];
1393:     const PetscInt rev = o >= 0 ? 0 : 1;
1394:     const PetscInt off = rev ? -(o+1) : o;

1396:     if (useCone) {
1397:       DMPlexGetConeSize(dm, q, &tmpSize);
1398:       DMPlexGetCone(dm, q, &tmp);
1399:       DMPlexGetConeOrientation(dm, q, &tmpO);
1400:     } else {
1401:       DMPlexGetSupportSize(dm, q, &tmpSize);
1402:       DMPlexGetSupport(dm, q, &tmp);
1403:       tmpO = NULL;
1404:     }
1405:     for (t = 0; t < tmpSize; ++t) {
1406:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1407:       const PetscInt cp = tmp[i];
1408:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1409:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1410:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1411:       PetscInt       co = tmpO ? tmpO[i] : 0;
1412:       PetscInt       c;

1414:       if (rev) {
1415:         PetscInt childSize, coff;
1416:         DMPlexGetConeSize(dm, cp, &childSize);
1417:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1418:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1419:       }
1420:       /* Check for duplicate */
1421:       for (c = 0; c < closureSize; c += 2) {
1422:         if (closure[c] == cp) break;
1423:       }
1424:       if (c == closureSize) {
1425:         closure[closureSize]   = cp;
1426:         closure[closureSize+1] = co;
1427:         fifo[fifoSize]         = cp;
1428:         fifo[fifoSize+1]       = co;
1429:         closureSize           += 2;
1430:         fifoSize              += 2;
1431:       }
1432:     }
1433:     fifoStart += 2;
1434:   }
1435:   if (numPoints) *numPoints = closureSize/2;
1436:   if (points)    *points    = closure;
1437:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1438:   return(0);
1439: }

1443: /*@C
1444:   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

1446:   Not collective

1448:   Input Parameters:
1449: + mesh - The DMPlex
1450: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1451: . orientation - The orientation of the point
1452: . useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
1453: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used

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

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

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

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

1468:   Level: beginner

1470: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1471: @*/
1472: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1473: {
1474:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1475:   PetscInt       *closure, *fifo;
1476:   const PetscInt *tmp = NULL, *tmpO = NULL;
1477:   PetscInt        tmpSize, t;
1478:   PetscInt        depth       = 0, maxSize;
1479:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1480:   PetscErrorCode  ierr;

1484:   DMPlexGetDepth(dm, &depth);
1485:   /* This is only 1-level */
1486:   if (useCone) {
1487:     DMPlexGetConeSize(dm, p, &tmpSize);
1488:     DMPlexGetCone(dm, p, &tmp);
1489:     DMPlexGetConeOrientation(dm, p, &tmpO);
1490:   } else {
1491:     DMPlexGetSupportSize(dm, p, &tmpSize);
1492:     DMPlexGetSupport(dm, p, &tmp);
1493:   }
1494:   if (depth == 1) {
1495:     if (*points) {
1496:       closure = *points;
1497:     } else {
1498:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1499:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1500:     }
1501:     closure[0] = p; closure[1] = ornt;
1502:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1503:       const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1504:       closure[closureSize]   = tmp[i];
1505:       closure[closureSize+1] = tmpO ? tmpO[i] : 0;
1506:     }
1507:     if (numPoints) *numPoints = closureSize/2;
1508:     if (points)    *points    = closure;
1509:     return(0);
1510:   }
1511:   {
1512:     PetscInt c, coneSeries, s,supportSeries;

1514:     c = mesh->maxConeSize;
1515:     coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1516:     s = mesh->maxSupportSize;
1517:     supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1518:     maxSize = 2*PetscMax(coneSeries,supportSeries);
1519:   }
1520:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1521:   if (*points) {
1522:     closure = *points;
1523:   } else {
1524:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1525:   }
1526:   closure[0] = p; closure[1] = ornt;
1527:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1528:     const PetscInt i  = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1529:     const PetscInt cp = tmp[i];
1530:     PetscInt       co = tmpO ? tmpO[i] : 0;

1532:     if (ornt < 0) {
1533:       PetscInt childSize, coff;
1534:       DMPlexGetConeSize(dm, cp, &childSize);
1535:       coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1536:       co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1537:     }
1538:     closure[closureSize]   = cp;
1539:     closure[closureSize+1] = co;
1540:     fifo[fifoSize]         = cp;
1541:     fifo[fifoSize+1]       = co;
1542:   }
1543:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1544:   while (fifoSize - fifoStart) {
1545:     const PetscInt q   = fifo[fifoStart];
1546:     const PetscInt o   = fifo[fifoStart+1];
1547:     const PetscInt rev = o >= 0 ? 0 : 1;
1548:     const PetscInt off = rev ? -(o+1) : o;

1550:     if (useCone) {
1551:       DMPlexGetConeSize(dm, q, &tmpSize);
1552:       DMPlexGetCone(dm, q, &tmp);
1553:       DMPlexGetConeOrientation(dm, q, &tmpO);
1554:     } else {
1555:       DMPlexGetSupportSize(dm, q, &tmpSize);
1556:       DMPlexGetSupport(dm, q, &tmp);
1557:       tmpO = NULL;
1558:     }
1559:     for (t = 0; t < tmpSize; ++t) {
1560:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1561:       const PetscInt cp = tmp[i];
1562:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1563:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1564:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1565:       PetscInt       co = tmpO ? tmpO[i] : 0;
1566:       PetscInt       c;

1568:       if (rev) {
1569:         PetscInt childSize, coff;
1570:         DMPlexGetConeSize(dm, cp, &childSize);
1571:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1572:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1573:       }
1574:       /* Check for duplicate */
1575:       for (c = 0; c < closureSize; c += 2) {
1576:         if (closure[c] == cp) break;
1577:       }
1578:       if (c == closureSize) {
1579:         closure[closureSize]   = cp;
1580:         closure[closureSize+1] = co;
1581:         fifo[fifoSize]         = cp;
1582:         fifo[fifoSize+1]       = co;
1583:         closureSize           += 2;
1584:         fifoSize              += 2;
1585:       }
1586:     }
1587:     fifoStart += 2;
1588:   }
1589:   if (numPoints) *numPoints = closureSize/2;
1590:   if (points)    *points    = closure;
1591:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1592:   return(0);
1593: }

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

1600:   Not collective

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

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

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

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

1618:   Level: beginner

1620: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1621: @*/
1622: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1623: {

1630:   DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1631:   if (numPoints) *numPoints = 0;
1632:   return(0);
1633: }

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

1640:   Not collective

1642:   Input Parameter:
1643: . mesh - The DMPlex

1645:   Output Parameters:
1646: + maxConeSize - The maximum number of in-edges
1647: - maxSupportSize - The maximum number of out-edges

1649:   Level: beginner

1651: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1652: @*/
1653: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1654: {
1655:   DM_Plex *mesh = (DM_Plex*) dm->data;

1659:   if (maxConeSize)    *maxConeSize    = mesh->maxConeSize;
1660:   if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1661:   return(0);
1662: }

1666: PetscErrorCode DMSetUp_Plex(DM dm)
1667: {
1668:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1669:   PetscInt       size;

1674:   PetscSectionSetUp(mesh->coneSection);
1675:   PetscSectionGetStorageSize(mesh->coneSection, &size);
1676:   PetscMalloc1(size, &mesh->cones);
1677:   PetscCalloc1(size, &mesh->coneOrientations);
1678:   if (mesh->maxSupportSize) {
1679:     PetscSectionSetUp(mesh->supportSection);
1680:     PetscSectionGetStorageSize(mesh->supportSection, &size);
1681:     PetscMalloc1(size, &mesh->supports);
1682:   }
1683:   return(0);
1684: }

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

1693:   if (subdm) {DMClone(dm, subdm);}
1694:   DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1695:   return(0);
1696: }

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

1703:   Not collective

1705:   Input Parameter:
1706: . mesh - The DMPlex

1708:   Output Parameter:

1710:   Note:
1711:   This should be called after all calls to DMPlexSetCone()

1713:   Level: beginner

1715: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1716: @*/
1717: PetscErrorCode DMPlexSymmetrize(DM dm)
1718: {
1719:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1720:   PetscInt      *offsets;
1721:   PetscInt       supportSize;
1722:   PetscInt       pStart, pEnd, p;

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

1733:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1734:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1735:     for (c = off; c < off+dof; ++c) {
1736:       PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1737:     }
1738:   }
1739:   for (p = pStart; p < pEnd; ++p) {
1740:     PetscInt dof;

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

1744:     mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1745:   }
1746:   PetscSectionSetUp(mesh->supportSection);
1747:   /* Calculate supports */
1748:   PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1749:   PetscMalloc1(supportSize, &mesh->supports);
1750:   PetscCalloc1(pEnd - pStart, &offsets);
1751:   for (p = pStart; p < pEnd; ++p) {
1752:     PetscInt dof, off, c;

1754:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1755:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1756:     for (c = off; c < off+dof; ++c) {
1757:       const PetscInt q = mesh->cones[c];
1758:       PetscInt       offS;

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

1762:       mesh->supports[offS+offsets[q]] = p;
1763:       ++offsets[q];
1764:     }
1765:   }
1766:   PetscFree(offsets);
1767:   return(0);
1768: }

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

1778:   Not collective

1780:   Input Parameter:
1781: . mesh - The DMPlex

1783:   Output Parameter:

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

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

1793:   Level: beginner

1795: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1796: @*/
1797: PetscErrorCode DMPlexStratify(DM dm)
1798: {
1799:   DMLabel        label;
1800:   PetscInt       pStart, pEnd, p;
1801:   PetscInt       numRoots = 0, numLeaves = 0;

1806:   PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1807:   /* Calculate depth */
1808:   DMPlexGetChart(dm, &pStart, &pEnd);
1809:   DMPlexCreateLabel(dm, "depth");
1810:   DMPlexGetDepthLabel(dm, &label);
1811:   /* Initialize roots and count leaves */
1812:   for (p = pStart; p < pEnd; ++p) {
1813:     PetscInt coneSize, supportSize;

1815:     DMPlexGetConeSize(dm, p, &coneSize);
1816:     DMPlexGetSupportSize(dm, p, &supportSize);
1817:     if (!coneSize && supportSize) {
1818:       ++numRoots;
1819:       DMLabelSetValue(label, p, 0);
1820:     } else if (!supportSize && coneSize) {
1821:       ++numLeaves;
1822:     } else if (!supportSize && !coneSize) {
1823:       /* Isolated points */
1824:       DMLabelSetValue(label, p, 0);
1825:     }
1826:   }
1827:   if (numRoots + numLeaves == (pEnd - pStart)) {
1828:     for (p = pStart; p < pEnd; ++p) {
1829:       PetscInt coneSize, supportSize;

1831:       DMPlexGetConeSize(dm, p, &coneSize);
1832:       DMPlexGetSupportSize(dm, p, &supportSize);
1833:       if (!supportSize && coneSize) {
1834:         DMLabelSetValue(label, p, 1);
1835:       }
1836:     }
1837:   } else {
1838:     IS       pointIS;
1839:     PetscInt numPoints = 0, level = 0;

1841:     DMLabelGetStratumIS(label, level, &pointIS);
1842:     if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1843:     while (numPoints) {
1844:       const PetscInt *points;
1845:       const PetscInt  newLevel = level+1;

1847:       ISGetIndices(pointIS, &points);
1848:       for (p = 0; p < numPoints; ++p) {
1849:         const PetscInt  point = points[p];
1850:         const PetscInt *support;
1851:         PetscInt        supportSize, s;

1853:         DMPlexGetSupportSize(dm, point, &supportSize);
1854:         DMPlexGetSupport(dm, point, &support);
1855:         for (s = 0; s < supportSize; ++s) {
1856:           DMLabelSetValue(label, support[s], newLevel);
1857:         }
1858:       }
1859:       ISRestoreIndices(pointIS, &points);
1860:       ++level;
1861:       ISDestroy(&pointIS);
1862:       DMLabelGetStratumIS(label, level, &pointIS);
1863:       if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1864:       else         {numPoints = 0;}
1865:     }
1866:     ISDestroy(&pointIS);
1867:   }
1868:   PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1869:   return(0);
1870: }

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

1877:   Not Collective

1879:   Input Parameters:
1880: + dm - The DMPlex object
1881: . numPoints - The number of input points for the join
1882: - points - The input points

1884:   Output Parameters:
1885: + numCoveredPoints - The number of points in the join
1886: - coveredPoints - The points in the join

1888:   Level: intermediate

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

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

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

1898: .keywords: mesh
1899: .seealso: DMPlexRestoreJoin(), DMPlexGetMeet()
1900: @*/
1901: PetscErrorCode DMPlexGetJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1902: {
1903:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1904:   PetscInt      *join[2];
1905:   PetscInt       joinSize, i = 0;
1906:   PetscInt       dof, off, p, c, m;

1914:   DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[0]);
1915:   DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1]);
1916:   /* Copy in support of first point */
1917:   PetscSectionGetDof(mesh->supportSection, points[0], &dof);
1918:   PetscSectionGetOffset(mesh->supportSection, points[0], &off);
1919:   for (joinSize = 0; joinSize < dof; ++joinSize) {
1920:     join[i][joinSize] = mesh->supports[off+joinSize];
1921:   }
1922:   /* Check each successive support */
1923:   for (p = 1; p < numPoints; ++p) {
1924:     PetscInt newJoinSize = 0;

1926:     PetscSectionGetDof(mesh->supportSection, points[p], &dof);
1927:     PetscSectionGetOffset(mesh->supportSection, points[p], &off);
1928:     for (c = 0; c < dof; ++c) {
1929:       const PetscInt point = mesh->supports[off+c];

1931:       for (m = 0; m < joinSize; ++m) {
1932:         if (point == join[i][m]) {
1933:           join[1-i][newJoinSize++] = point;
1934:           break;
1935:         }
1936:       }
1937:     }
1938:     joinSize = newJoinSize;
1939:     i        = 1-i;
1940:   }
1941:   *numCoveredPoints = joinSize;
1942:   *coveredPoints    = join[i];
1943:   DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
1944:   return(0);
1945: }

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

1952:   Not Collective

1954:   Input Parameters:
1955: + dm - The DMPlex object
1956: . numPoints - The number of input points for the join
1957: - points - The input points

1959:   Output Parameters:
1960: + numCoveredPoints - The number of points in the join
1961: - coveredPoints - The points in the join

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

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

1969:   Level: intermediate

1971: .keywords: mesh
1972: .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet()
1973: @*/
1974: PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1975: {

1983:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
1984:   if (numCoveredPoints) *numCoveredPoints = 0;
1985:   return(0);
1986: }

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

1993:   Not Collective

1995:   Input Parameters:
1996: + dm - The DMPlex object
1997: . numPoints - The number of input points for the join
1998: - points - The input points

2000:   Output Parameters:
2001: + numCoveredPoints - The number of points in the join
2002: - coveredPoints - The points in the join

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

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

2010:   Level: intermediate

2012: .keywords: mesh
2013: .seealso: DMPlexGetJoin(), DMPlexRestoreJoin(), DMPlexGetMeet()
2014: @*/
2015: PetscErrorCode DMPlexGetFullJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2016: {
2017:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2018:   PetscInt      *offsets, **closures;
2019:   PetscInt      *join[2];
2020:   PetscInt       depth = 0, maxSize, joinSize = 0, i = 0;
2021:   PetscInt       p, d, c, m, ms;


2030:   DMPlexGetDepth(dm, &depth);
2031:   PetscCalloc1(numPoints, &closures);
2032:   DMGetWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2033:   ms      = mesh->maxSupportSize;
2034:   maxSize = (ms > 1) ? ((PetscPowInt(ms,depth+1)-1)/(ms-1)) : depth + 1;
2035:   DMGetWorkArray(dm, maxSize, PETSC_INT, &join[0]);
2036:   DMGetWorkArray(dm, maxSize, PETSC_INT, &join[1]);

2038:   for (p = 0; p < numPoints; ++p) {
2039:     PetscInt closureSize;

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

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

2047:       DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
2048:       for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) {
2049:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2050:           offsets[p*(depth+2)+d+1] = i;
2051:           break;
2052:         }
2053:       }
2054:       if (i == closureSize) offsets[p*(depth+2)+d+1] = i;
2055:     }
2056:     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);
2057:   }
2058:   for (d = 0; d < depth+1; ++d) {
2059:     PetscInt dof;

2061:     /* Copy in support of first point */
2062:     dof = offsets[d+1] - offsets[d];
2063:     for (joinSize = 0; joinSize < dof; ++joinSize) {
2064:       join[i][joinSize] = closures[0][(offsets[d]+joinSize)*2];
2065:     }
2066:     /* Check each successive cone */
2067:     for (p = 1; p < numPoints && joinSize; ++p) {
2068:       PetscInt newJoinSize = 0;

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

2074:         for (m = 0; m < joinSize; ++m) {
2075:           if (point == join[i][m]) {
2076:             join[1-i][newJoinSize++] = point;
2077:             break;
2078:           }
2079:         }
2080:       }
2081:       joinSize = newJoinSize;
2082:       i        = 1-i;
2083:     }
2084:     if (joinSize) break;
2085:   }
2086:   *numCoveredPoints = joinSize;
2087:   *coveredPoints    = join[i];
2088:   for (p = 0; p < numPoints; ++p) {
2089:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, NULL, &closures[p]);
2090:   }
2091:   PetscFree(closures);
2092:   DMRestoreWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2093:   DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2094:   return(0);
2095: }

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

2102:   Not Collective

2104:   Input Parameters:
2105: + dm - The DMPlex object
2106: . numPoints - The number of input points for the meet
2107: - points - The input points

2109:   Output Parameters:
2110: + numCoveredPoints - The number of points in the meet
2111: - coveredPoints - The points in the meet

2113:   Level: intermediate

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

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

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

2123: .keywords: mesh
2124: .seealso: DMPlexRestoreMeet(), DMPlexGetJoin()
2125: @*/
2126: PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2127: {
2128:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2129:   PetscInt      *meet[2];
2130:   PetscInt       meetSize, i = 0;
2131:   PetscInt       dof, off, p, c, m;

2139:   DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[0]);
2140:   DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1]);
2141:   /* Copy in cone of first point */
2142:   PetscSectionGetDof(mesh->coneSection, points[0], &dof);
2143:   PetscSectionGetOffset(mesh->coneSection, points[0], &off);
2144:   for (meetSize = 0; meetSize < dof; ++meetSize) {
2145:     meet[i][meetSize] = mesh->cones[off+meetSize];
2146:   }
2147:   /* Check each successive cone */
2148:   for (p = 1; p < numPoints; ++p) {
2149:     PetscInt newMeetSize = 0;

2151:     PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2152:     PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2153:     for (c = 0; c < dof; ++c) {
2154:       const PetscInt point = mesh->cones[off+c];

2156:       for (m = 0; m < meetSize; ++m) {
2157:         if (point == meet[i][m]) {
2158:           meet[1-i][newMeetSize++] = point;
2159:           break;
2160:         }
2161:       }
2162:     }
2163:     meetSize = newMeetSize;
2164:     i        = 1-i;
2165:   }
2166:   *numCoveringPoints = meetSize;
2167:   *coveringPoints    = meet[i];
2168:   DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2169:   return(0);
2170: }

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

2177:   Not Collective

2179:   Input Parameters:
2180: + dm - The DMPlex object
2181: . numPoints - The number of input points for the meet
2182: - points - The input points

2184:   Output Parameters:
2185: + numCoveredPoints - The number of points in the meet
2186: - coveredPoints - The points in the meet

2188:   Level: intermediate

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

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

2196: .keywords: mesh
2197: .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin()
2198: @*/
2199: PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2200: {

2208:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2209:   if (numCoveredPoints) *numCoveredPoints = 0;
2210:   return(0);
2211: }

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

2218:   Not Collective

2220:   Input Parameters:
2221: + dm - The DMPlex object
2222: . numPoints - The number of input points for the meet
2223: - points - The input points

2225:   Output Parameters:
2226: + numCoveredPoints - The number of points in the meet
2227: - coveredPoints - The points in the meet

2229:   Level: intermediate

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

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

2237: .keywords: mesh
2238: .seealso: DMPlexGetMeet(), DMPlexRestoreMeet(), DMPlexGetJoin()
2239: @*/
2240: PetscErrorCode DMPlexGetFullMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2241: {
2242:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2243:   PetscInt      *offsets, **closures;
2244:   PetscInt      *meet[2];
2245:   PetscInt       height = 0, maxSize, meetSize = 0, i = 0;
2246:   PetscInt       p, h, c, m, mc;


2255:   DMPlexGetDepth(dm, &height);
2256:   PetscMalloc1(numPoints, &closures);
2257:   DMGetWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2258:   mc      = mesh->maxConeSize;
2259:   maxSize = (mc > 1) ? ((PetscPowInt(mc,height+1)-1)/(mc-1)) : height + 1;
2260:   DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[0]);
2261:   DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[1]);

2263:   for (p = 0; p < numPoints; ++p) {
2264:     PetscInt closureSize;

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

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

2272:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
2273:       for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) {
2274:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2275:           offsets[p*(height+2)+h+1] = i;
2276:           break;
2277:         }
2278:       }
2279:       if (i == closureSize) offsets[p*(height+2)+h+1] = i;
2280:     }
2281:     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);
2282:   }
2283:   for (h = 0; h < height+1; ++h) {
2284:     PetscInt dof;

2286:     /* Copy in cone of first point */
2287:     dof = offsets[h+1] - offsets[h];
2288:     for (meetSize = 0; meetSize < dof; ++meetSize) {
2289:       meet[i][meetSize] = closures[0][(offsets[h]+meetSize)*2];
2290:     }
2291:     /* Check each successive cone */
2292:     for (p = 1; p < numPoints && meetSize; ++p) {
2293:       PetscInt newMeetSize = 0;

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

2299:         for (m = 0; m < meetSize; ++m) {
2300:           if (point == meet[i][m]) {
2301:             meet[1-i][newMeetSize++] = point;
2302:             break;
2303:           }
2304:         }
2305:       }
2306:       meetSize = newMeetSize;
2307:       i        = 1-i;
2308:     }
2309:     if (meetSize) break;
2310:   }
2311:   *numCoveredPoints = meetSize;
2312:   *coveredPoints    = meet[i];
2313:   for (p = 0; p < numPoints; ++p) {
2314:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, NULL, &closures[p]);
2315:   }
2316:   PetscFree(closures);
2317:   DMRestoreWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2318:   DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2319:   return(0);
2320: }

2324: /*@C
2325:   DMPlexEqual - Determine if two DMs have the same topology

2327:   Not Collective

2329:   Input Parameters:
2330: + dmA - A DMPlex object
2331: - dmB - A DMPlex object

2333:   Output Parameters:
2334: . equal - PETSC_TRUE if the topologies are identical

2336:   Level: intermediate

2338:   Notes:
2339:   We are not solving graph isomorphism, so we do not permutation.

2341: .keywords: mesh
2342: .seealso: DMPlexGetCone()
2343: @*/
2344: PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal)
2345: {
2346:   PetscInt       depth, depthB, pStart, pEnd, pStartB, pEndB, p;


2354:   *equal = PETSC_FALSE;
2355:   DMPlexGetDepth(dmA, &depth);
2356:   DMPlexGetDepth(dmB, &depthB);
2357:   if (depth != depthB) return(0);
2358:   DMPlexGetChart(dmA, &pStart,  &pEnd);
2359:   DMPlexGetChart(dmB, &pStartB, &pEndB);
2360:   if ((pStart != pStartB) || (pEnd != pEndB)) return(0);
2361:   for (p = pStart; p < pEnd; ++p) {
2362:     const PetscInt *cone, *coneB, *ornt, *orntB, *support, *supportB;
2363:     PetscInt        coneSize, coneSizeB, c, supportSize, supportSizeB, s;

2365:     DMPlexGetConeSize(dmA, p, &coneSize);
2366:     DMPlexGetCone(dmA, p, &cone);
2367:     DMPlexGetConeOrientation(dmA, p, &ornt);
2368:     DMPlexGetConeSize(dmB, p, &coneSizeB);
2369:     DMPlexGetCone(dmB, p, &coneB);
2370:     DMPlexGetConeOrientation(dmB, p, &orntB);
2371:     if (coneSize != coneSizeB) return(0);
2372:     for (c = 0; c < coneSize; ++c) {
2373:       if (cone[c] != coneB[c]) return(0);
2374:       if (ornt[c] != orntB[c]) return(0);
2375:     }
2376:     DMPlexGetSupportSize(dmA, p, &supportSize);
2377:     DMPlexGetSupport(dmA, p, &support);
2378:     DMPlexGetSupportSize(dmB, p, &supportSizeB);
2379:     DMPlexGetSupport(dmB, p, &supportB);
2380:     if (supportSize != supportSizeB) return(0);
2381:     for (s = 0; s < supportSize; ++s) {
2382:       if (support[s] != supportB[s]) return(0);
2383:     }
2384:   }
2385:   *equal = PETSC_TRUE;
2386:   return(0);
2387: }

2391: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2392: {
2393:   MPI_Comm       comm;

2397:   PetscObjectGetComm((PetscObject)dm,&comm);
2399:   switch (cellDim) {
2400:   case 0:
2401:     *numFaceVertices = 0;
2402:     break;
2403:   case 1:
2404:     *numFaceVertices = 1;
2405:     break;
2406:   case 2:
2407:     switch (numCorners) {
2408:     case 3: /* triangle */
2409:       *numFaceVertices = 2; /* Edge has 2 vertices */
2410:       break;
2411:     case 4: /* quadrilateral */
2412:       *numFaceVertices = 2; /* Edge has 2 vertices */
2413:       break;
2414:     case 6: /* quadratic triangle, tri and quad cohesive Lagrange cells */
2415:       *numFaceVertices = 3; /* Edge has 3 vertices */
2416:       break;
2417:     case 9: /* quadratic quadrilateral, quadratic quad cohesive Lagrange cells */
2418:       *numFaceVertices = 3; /* Edge has 3 vertices */
2419:       break;
2420:     default:
2421:       SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim);
2422:     }
2423:     break;
2424:   case 3:
2425:     switch (numCorners) {
2426:     case 4: /* tetradehdron */
2427:       *numFaceVertices = 3; /* Face has 3 vertices */
2428:       break;
2429:     case 6: /* tet cohesive cells */
2430:       *numFaceVertices = 4; /* Face has 4 vertices */
2431:       break;
2432:     case 8: /* hexahedron */
2433:       *numFaceVertices = 4; /* Face has 4 vertices */
2434:       break;
2435:     case 9: /* tet cohesive Lagrange cells */
2436:       *numFaceVertices = 6; /* Face has 6 vertices */
2437:       break;
2438:     case 10: /* quadratic tetrahedron */
2439:       *numFaceVertices = 6; /* Face has 6 vertices */
2440:       break;
2441:     case 12: /* hex cohesive Lagrange cells */
2442:       *numFaceVertices = 6; /* Face has 6 vertices */
2443:       break;
2444:     case 18: /* quadratic tet cohesive Lagrange cells */
2445:       *numFaceVertices = 6; /* Face has 6 vertices */
2446:       break;
2447:     case 27: /* quadratic hexahedron, quadratic hex cohesive Lagrange cells */
2448:       *numFaceVertices = 9; /* Face has 9 vertices */
2449:       break;
2450:     default:
2451:       SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim);
2452:     }
2453:     break;
2454:   default:
2455:     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid cell dimension %d", cellDim);
2456:   }
2457:   return(0);
2458: }

2462: PetscErrorCode DMPlexLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2463: {
2464:   PetscInt d;

2467:   if (!dm->maxCell) {
2468:     for (d = 0; d < dim; ++d) out[d] = in[d];
2469:   } else {
2470:     for (d = 0; d < dim; ++d) {
2471:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2472:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2473:       } else {
2474:         out[d] = in[d];
2475:       }
2476:     }
2477:   }
2478:   return(0);
2479: }

2483: PetscErrorCode DMPlexLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2484: {
2485:   PetscInt d;

2488:   if (!dm->maxCell) {
2489:     for (d = 0; d < dim; ++d) out[d] += in[d];
2490:   } else {
2491:     for (d = 0; d < dim; ++d) {
2492:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2493:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2494:       } else {
2495:         out[d] += in[d];
2496:       }
2497:     }
2498:   }
2499:   return(0);
2500: }

2504: PetscErrorCode DMPlexLocalizeCoordinates(DM dm)
2505: {
2506:   PetscSection   coordSection, cSection;
2507:   Vec            coordinates,  cVec;
2508:   PetscScalar   *coords, *coords2, *anchor;
2509:   PetscInt       Nc, cStart, cEnd, c, vStart, vEnd, v, dof, d, off, off2, bs, coordSize;

2514:   if (!dm->maxCell) return(0);
2515:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2516:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2517:   DMGetCoordinatesLocal(dm, &coordinates);
2518:   DMGetCoordinateSection(dm, &coordSection);
2519:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
2520:   PetscSectionSetNumFields(cSection, 1);
2521:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
2522:   PetscSectionSetFieldComponents(cSection, 0, Nc);
2523:   PetscSectionSetChart(cSection, cStart, vEnd);
2524:   for (v = vStart; v < vEnd; ++v) {
2525:     PetscSectionGetDof(coordSection, v, &dof);
2526:     PetscSectionSetDof(cSection,     v,  dof);
2527:     PetscSectionSetFieldDof(cSection, v, 0, dof);
2528:   }
2529:   for (c = cStart; c < cEnd; ++c) {
2530:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, NULL);
2531:     PetscSectionSetDof(cSection, c, dof);
2532:     PetscSectionSetFieldDof(cSection, c, 0, dof);
2533:   }
2534:   PetscSectionSetUp(cSection);
2535:   PetscSectionGetStorageSize(cSection, &coordSize);
2536:   VecCreate(PetscObjectComm((PetscObject) dm), &cVec);
2537:   VecGetBlockSize(coordinates, &bs);
2538:   VecSetBlockSize(cVec,         bs);
2539:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
2540:   VecSetType(cVec,VECSTANDARD);
2541:   VecGetArray(coordinates, &coords);
2542:   VecGetArray(cVec,        &coords2);
2543:   for (v = vStart; v < vEnd; ++v) {
2544:     PetscSectionGetDof(coordSection, v, &dof);
2545:     PetscSectionGetOffset(coordSection, v, &off);
2546:     PetscSectionGetOffset(cSection,     v, &off2);
2547:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
2548:   }
2549:   DMGetWorkArray(dm, 3, PETSC_SCALAR, &anchor);
2550:   for (c = cStart; c < cEnd; ++c) {
2551:     PetscScalar *cellCoords = NULL;
2552:     PetscInt     b;

2554:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2555:     PetscSectionGetOffset(cSection, c, &off2);
2556:     for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
2557:     for (d = 0; d < dof/bs; ++d) {DMPlexLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
2558:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2559:   }
2560:   DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);
2561:   VecRestoreArray(coordinates, &coords);
2562:   VecRestoreArray(cVec,        &coords2);
2563:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
2564:   DMSetCoordinatesLocal(dm, cVec);
2565:   VecDestroy(&cVec);
2566:   PetscSectionDestroy(&cSection);
2567:   return(0);
2568: }

2572: /*@
2573:   DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point

2575:   Not Collective

2577:   Input Parameter:
2578: . dm    - The DMPlex object

2580:   Output Parameter:
2581: . depthLabel - The DMLabel recording point depth

2583:   Level: developer

2585: .keywords: mesh, points
2586: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2587: @*/
2588: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
2589: {
2590:   DM_Plex       *mesh = (DM_Plex*) dm->data;

2596:   if (!mesh->depthLabel) {DMPlexGetLabel(dm, "depth", &mesh->depthLabel);}
2597:   *depthLabel = mesh->depthLabel;
2598:   return(0);
2599: }

2603: /*@
2604:   DMPlexGetDepth - Get the depth of the DAG representing this mesh

2606:   Not Collective

2608:   Input Parameter:
2609: . dm    - The DMPlex object

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

2614:   Level: developer

2616: .keywords: mesh, points
2617: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2618: @*/
2619: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
2620: {
2621:   DMLabel        label;
2622:   PetscInt       d = 0;

2628:   DMPlexGetDepthLabel(dm, &label);
2629:   if (label) {DMLabelGetNumValues(label, &d);}
2630:   *depth = d-1;
2631:   return(0);
2632: }

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

2639:   Not Collective

2641:   Input Parameters:
2642: + dm           - The DMPlex object
2643: - stratumValue - The requested depth

2645:   Output Parameters:
2646: + start - The first point at this depth
2647: - end   - One beyond the last point at this depth

2649:   Level: developer

2651: .keywords: mesh, points
2652: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
2653: @*/
2654: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2655: {
2656:   DMLabel        label;
2657:   PetscInt       pStart, pEnd;

2664:   DMPlexGetChart(dm, &pStart, &pEnd);
2665:   if (pStart == pEnd) return(0);
2666:   if (stratumValue < 0) {
2667:     if (start) *start = pStart;
2668:     if (end)   *end   = pEnd;
2669:     return(0);
2670:   }
2671:   DMPlexGetDepthLabel(dm, &label);
2672:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2673:   DMLabelGetStratumBounds(label, stratumValue, start, end);
2674:   return(0);
2675: }

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

2682:   Not Collective

2684:   Input Parameters:
2685: + dm           - The DMPlex object
2686: - stratumValue - The requested height

2688:   Output Parameters:
2689: + start - The first point at this height
2690: - end   - One beyond the last point at this height

2692:   Level: developer

2694: .keywords: mesh, points
2695: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
2696: @*/
2697: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2698: {
2699:   DMLabel        label;
2700:   PetscInt       depth, pStart, pEnd;

2707:   DMPlexGetChart(dm, &pStart, &pEnd);
2708:   if (pStart == pEnd) return(0);
2709:   if (stratumValue < 0) {
2710:     if (start) *start = pStart;
2711:     if (end)   *end   = pEnd;
2712:     return(0);
2713:   }
2714:   DMPlexGetDepthLabel(dm, &label);
2715:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2716:   DMLabelGetNumValues(label, &depth);
2717:   DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
2718:   return(0);
2719: }

2723: /* Set the number of dof on each point and separate by fields */
2724: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
2725: {
2726:   PetscInt      *pMax;
2727:   PetscInt       depth, pStart = 0, pEnd = 0;
2728:   PetscInt       Nf, p, d, dep, f;
2729:   PetscBool     *isFE;

2733:   PetscMalloc1(numFields, &isFE);
2734:   DMGetNumFields(dm, &Nf);
2735:   for (f = 0; f < numFields; ++f) {
2736:     PetscObject  obj;
2737:     PetscClassId id;

2739:     isFE[f] = PETSC_FALSE;
2740:     if (f >= Nf) continue;
2741:     DMGetField(dm, f, &obj);
2742:     PetscObjectGetClassId(obj, &id);
2743:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
2744:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
2745:   }
2746:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
2747:   if (numFields > 0) {
2748:     PetscSectionSetNumFields(*section, numFields);
2749:     if (numComp) {
2750:       for (f = 0; f < numFields; ++f) {
2751:         PetscSectionSetFieldComponents(*section, f, numComp[f]);
2752:       }
2753:     }
2754:   }
2755:   DMPlexGetChart(dm, &pStart, &pEnd);
2756:   PetscSectionSetChart(*section, pStart, pEnd);
2757:   DMPlexGetDepth(dm, &depth);
2758:   PetscMalloc1(depth+1,&pMax);
2759:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2760:   for (dep = 0; dep <= depth; ++dep) {
2761:     d    = dim == depth ? dep : (!dep ? 0 : dim);
2762:     DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
2763:     pMax[dep] = pMax[dep] < 0 ? pEnd : pMax[dep];
2764:     for (p = pStart; p < pEnd; ++p) {
2765:       PetscInt tot = 0;

2767:       for (f = 0; f < numFields; ++f) {
2768:         if (isFE[f] && p >= pMax[dep]) continue;
2769:         PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
2770:         tot += numDof[f*(dim+1)+d];
2771:       }
2772:       PetscSectionSetDof(*section, p, tot);
2773:     }
2774:   }
2775:   PetscFree(pMax);
2776:   PetscFree(isFE);
2777:   return(0);
2778: }

2782: /* Set the number of dof on each point and separate by fields
2783:    If constDof is PETSC_DETERMINE, constrain every dof on the point
2784: */
2785: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC,const PetscInt bcField[],const IS bcPoints[], PetscInt constDof, PetscSection section)
2786: {
2787:   PetscInt       numFields;
2788:   PetscInt       bc;
2789:   PetscSection   aSec;

2793:   PetscSectionGetNumFields(section, &numFields);
2794:   for (bc = 0; bc < numBC; ++bc) {
2795:     PetscInt        field = 0;
2796:     const PetscInt *idx;
2797:     PetscInt        n, i;

2799:     if (numFields) field = bcField[bc];
2800:     ISGetLocalSize(bcPoints[bc], &n);
2801:     ISGetIndices(bcPoints[bc], &idx);
2802:     for (i = 0; i < n; ++i) {
2803:       const PetscInt p        = idx[i];
2804:       PetscInt       numConst = constDof;

2806:       /* Constrain every dof on the point */
2807:       if (numConst < 0) {
2808:         if (numFields) {
2809:           PetscSectionGetFieldDof(section, p, field, &numConst);
2810:         } else {
2811:           PetscSectionGetDof(section, p, &numConst);
2812:         }
2813:       }
2814:       if (numFields) {
2815:         PetscSectionAddFieldConstraintDof(section, p, field, numConst);
2816:       }
2817:       PetscSectionAddConstraintDof(section, p, numConst);
2818:     }
2819:     ISRestoreIndices(bcPoints[bc], &idx);
2820:   }
2821:   DMPlexGetAnchors(dm, &aSec, NULL);
2822:   if (aSec) {
2823:     PetscInt aStart, aEnd, a;

2825:     PetscSectionGetChart(aSec, &aStart, &aEnd);
2826:     for (a = aStart; a < aEnd; a++) {
2827:       PetscInt dof;

2829:       PetscSectionGetDof(aSec, a, &dof);
2830:       if (dof) {
2831:         /* if there are point-to-point constraints, then all dofs are constrained */
2832:         PetscSectionGetDof(section, a, &dof);
2833:         PetscSectionSetConstraintDof(section, a, dof);
2834:         if (numFields) {
2835:           PetscInt f;

2837:           for (f = 0; f < numFields; f++) {
2838:             PetscSectionGetFieldDof(section, a, f, &dof);
2839:             PetscSectionSetFieldConstraintDof(section, a, f, dof);
2840:           }
2841:         }
2842:       }
2843:     }
2844:   }
2845:   return(0);
2846: }

2850: /* Set the constrained indices on each point and separate by fields */
2851: PetscErrorCode DMPlexCreateSectionBCIndicesAll(DM dm, PetscSection section)
2852: {
2853:   PetscInt      *maxConstraints;
2854:   PetscInt       numFields, f, pStart = 0, pEnd = 0, p;

2858:   PetscSectionGetNumFields(section, &numFields);
2859:   PetscSectionGetChart(section, &pStart, &pEnd);
2860:   PetscMalloc1(numFields+1, &maxConstraints);
2861:   for (f = 0; f <= numFields; ++f) maxConstraints[f] = 0;
2862:   for (p = pStart; p < pEnd; ++p) {
2863:     PetscInt cdof;

2865:     if (numFields) {
2866:       for (f = 0; f < numFields; ++f) {
2867:         PetscSectionGetFieldConstraintDof(section, p, f, &cdof);
2868:         maxConstraints[f] = PetscMax(maxConstraints[f], cdof);
2869:       }
2870:     } else {
2871:       PetscSectionGetConstraintDof(section, p, &cdof);
2872:       maxConstraints[0] = PetscMax(maxConstraints[0], cdof);
2873:     }
2874:   }
2875:   for (f = 0; f < numFields; ++f) {
2876:     maxConstraints[numFields] += maxConstraints[f];
2877:   }
2878:   if (maxConstraints[numFields]) {
2879:     PetscInt *indices;

2881:     PetscMalloc1(maxConstraints[numFields], &indices);
2882:     for (p = pStart; p < pEnd; ++p) {
2883:       PetscInt cdof, d;

2885:       PetscSectionGetConstraintDof(section, p, &cdof);
2886:       if (cdof) {
2887:         if (cdof > maxConstraints[numFields]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, "Likely memory corruption, point %D cDof %D > maxConstraints %D", p, cdof, maxConstraints[numFields]);
2888:         if (numFields) {
2889:           PetscInt numConst = 0, foff = 0;

2891:           for (f = 0; f < numFields; ++f) {
2892:             PetscInt cfdof, fdof;

2894:             PetscSectionGetFieldDof(section, p, f, &fdof);
2895:             PetscSectionGetFieldConstraintDof(section, p, f, &cfdof);
2896:             /* Change constraint numbering from absolute local dof number to field relative local dof number */
2897:             for (d = 0; d < cfdof; ++d) indices[numConst+d] = d;
2898:             PetscSectionSetFieldConstraintIndices(section, p, f, &indices[numConst]);
2899:             for (d = 0; d < cfdof; ++d) indices[numConst+d] += foff;
2900:             numConst += cfdof;
2901:             foff     += fdof;
2902:           }
2903:           if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
2904:         } else {
2905:           for (d = 0; d < cdof; ++d) indices[d] = d;
2906:         }
2907:         PetscSectionSetConstraintIndices(section, p, indices);
2908:       }
2909:     }
2910:     PetscFree(indices);
2911:   }
2912:   PetscFree(maxConstraints);
2913:   return(0);
2914: }

2918: /* Set the constrained field indices on each point */
2919: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt field, IS bcPoints, IS constraintIndices, PetscSection section)
2920: {
2921:   const PetscInt *points, *indices;
2922:   PetscInt        numFields, maxDof, numPoints, p, numConstraints;
2923:   PetscErrorCode  ierr;

2926:   PetscSectionGetNumFields(section, &numFields);
2927:   if ((field < 0) || (field >= numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, numFields);

2929:   ISGetLocalSize(bcPoints, &numPoints);
2930:   ISGetIndices(bcPoints, &points);
2931:   if (!constraintIndices) {
2932:     PetscInt *idx, i;

2934:     PetscSectionGetMaxDof(section, &maxDof);
2935:     PetscMalloc1(maxDof, &idx);
2936:     for (i = 0; i < maxDof; ++i) idx[i] = i;
2937:     for (p = 0; p < numPoints; ++p) {
2938:       PetscSectionSetFieldConstraintIndices(section, points[p], field, idx);
2939:     }
2940:     PetscFree(idx);
2941:   } else {
2942:     ISGetLocalSize(constraintIndices, &numConstraints);
2943:     ISGetIndices(constraintIndices, &indices);
2944:     for (p = 0; p < numPoints; ++p) {
2945:       PetscInt fcdof;

2947:       PetscSectionGetFieldConstraintDof(section, points[p], field, &fcdof);
2948:       if (fcdof != numConstraints) SETERRQ4(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Section point %d field %d has %d constraints, but yo ugave %d indices", p, field, fcdof, numConstraints);
2949:       PetscSectionSetFieldConstraintIndices(section, points[p], field, indices);
2950:     }
2951:     ISRestoreIndices(constraintIndices, &indices);
2952:   }
2953:   ISRestoreIndices(bcPoints, &points);
2954:   return(0);
2955: }

2959: /* Set the constrained indices on each point and separate by fields */
2960: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
2961: {
2962:   PetscInt      *indices;
2963:   PetscInt       numFields, maxDof, f, pStart = 0, pEnd = 0, p;

2967:   PetscSectionGetMaxDof(section, &maxDof);
2968:   PetscMalloc1(maxDof, &indices);
2969:   PetscSectionGetNumFields(section, &numFields);
2970:   if (!numFields) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This function only works after users have set field constraint indices.");
2971:   PetscSectionGetChart(section, &pStart, &pEnd);
2972:   for (p = pStart; p < pEnd; ++p) {
2973:     PetscInt cdof, d;

2975:     PetscSectionGetConstraintDof(section, p, &cdof);
2976:     if (cdof) {
2977:       PetscInt numConst = 0, foff = 0;

2979:       for (f = 0; f < numFields; ++f) {
2980:         const PetscInt *fcind;
2981:         PetscInt        fdof, fcdof;

2983:         PetscSectionGetFieldDof(section, p, f, &fdof);
2984:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
2985:         if (fcdof) {PetscSectionGetFieldConstraintIndices(section, p, f, &fcind);}
2986:         /* Change constraint numbering from field relative local dof number to absolute local dof number */
2987:         for (d = 0; d < fcdof; ++d) indices[numConst+d] = fcind[d]+foff;
2988:         foff     += fdof;
2989:         numConst += fcdof;
2990:       }
2991:       if (cdof != numConst) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
2992:       PetscSectionSetConstraintIndices(section, p, indices);
2993:     }
2994:   }
2995:   PetscFree(indices);
2996:   return(0);
2997: }

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

3004:   Not Collective

3006:   Input Parameters:
3007: + dm        - The DMPlex object
3008: . dim       - The spatial dimension of the problem
3009: . numFields - The number of fields in the problem
3010: . numComp   - An array of size numFields that holds the number of components for each field
3011: . numDof    - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
3012: . numBC     - The number of boundary conditions
3013: . bcField   - An array of size numBC giving the field number for each boundry condition
3014: . bcPoints  - An array of size numBC giving an IS holding the sieve points to which each boundary condition applies
3015: - perm      - Optional permutation of the chart, or NULL

3017:   Output Parameter:
3018: . section - The PetscSection object

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

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

3025:   Level: developer

3027:   Fortran Notes:
3028:   A Fortran 90 version is available as DMPlexCreateSectionF90()

3030: .keywords: mesh, elements
3031: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
3032: @*/
3033: PetscErrorCode DMPlexCreateSection(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscInt numBC,const PetscInt bcField[],const IS bcPoints[], IS perm, PetscSection *section)
3034: {
3035:   PetscSection   aSec;

3039:   DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
3040:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcPoints, PETSC_DETERMINE, *section);
3041:   if (perm) {PetscSectionSetPermutation(*section, perm);}
3042:   PetscSectionSetUp(*section);
3043:   DMPlexGetAnchors(dm,&aSec,NULL);
3044:   if (numBC || aSec) {DMPlexCreateSectionBCIndicesAll(dm, *section);}
3045:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
3046:   return(0);
3047: }

3051: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
3052: {
3053:   PetscSection   section;

3057:   DMClone(dm, cdm);
3058:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
3059:   DMSetDefaultSection(*cdm, section);
3060:   PetscSectionDestroy(&section);
3061:   return(0);
3062: }

3066: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
3067: {
3068:   DM_Plex *mesh = (DM_Plex*) dm->data;

3072:   if (section) *section = mesh->coneSection;
3073:   return(0);
3074: }

3078: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
3079: {
3080:   DM_Plex *mesh = (DM_Plex*) dm->data;

3084:   if (section) *section = mesh->supportSection;
3085:   return(0);
3086: }

3090: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
3091: {
3092:   DM_Plex *mesh = (DM_Plex*) dm->data;

3096:   if (cones) *cones = mesh->cones;
3097:   return(0);
3098: }

3102: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
3103: {
3104:   DM_Plex *mesh = (DM_Plex*) dm->data;

3108:   if (coneOrientations) *coneOrientations = mesh->coneOrientations;
3109:   return(0);
3110: }

3112: /******************************** FEM Support **********************************/

3116: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3117: {
3118:   PetscScalar    *array, *vArray;
3119:   const PetscInt *cone, *coneO;
3120:   PetscInt        pStart, pEnd, p, numPoints, size = 0, offset = 0;
3121:   PetscErrorCode  ierr;

3124:   PetscSectionGetChart(section, &pStart, &pEnd);
3125:   DMPlexGetConeSize(dm, point, &numPoints);
3126:   DMPlexGetCone(dm, point, &cone);
3127:   DMPlexGetConeOrientation(dm, point, &coneO);
3128:   if (!values || !*values) {
3129:     if ((point >= pStart) && (point < pEnd)) {
3130:       PetscInt dof;

3132:       PetscSectionGetDof(section, point, &dof);
3133:       size += dof;
3134:     }
3135:     for (p = 0; p < numPoints; ++p) {
3136:       const PetscInt cp = cone[p];
3137:       PetscInt       dof;

3139:       if ((cp < pStart) || (cp >= pEnd)) continue;
3140:       PetscSectionGetDof(section, cp, &dof);
3141:       size += dof;
3142:     }
3143:     if (!values) {
3144:       if (csize) *csize = size;
3145:       return(0);
3146:     }
3147:     DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
3148:   } else {
3149:     array = *values;
3150:   }
3151:   size = 0;
3152:   VecGetArray(v, &vArray);
3153:   if ((point >= pStart) && (point < pEnd)) {
3154:     PetscInt     dof, off, d;
3155:     PetscScalar *varr;

3157:     PetscSectionGetDof(section, point, &dof);
3158:     PetscSectionGetOffset(section, point, &off);
3159:     varr = &vArray[off];
3160:     for (d = 0; d < dof; ++d, ++offset) {
3161:       array[offset] = varr[d];
3162:     }
3163:     size += dof;
3164:   }
3165:   for (p = 0; p < numPoints; ++p) {
3166:     const PetscInt cp = cone[p];
3167:     PetscInt       o  = coneO[p];
3168:     PetscInt       dof, off, d;
3169:     PetscScalar   *varr;

3171:     if ((cp < pStart) || (cp >= pEnd)) continue;
3172:     PetscSectionGetDof(section, cp, &dof);
3173:     PetscSectionGetOffset(section, cp, &off);
3174:     varr = &vArray[off];
3175:     if (o >= 0) {
3176:       for (d = 0; d < dof; ++d, ++offset) {
3177:         array[offset] = varr[d];
3178:       }
3179:     } else {
3180:       for (d = dof-1; d >= 0; --d, ++offset) {
3181:         array[offset] = varr[d];
3182:       }
3183:     }
3184:     size += dof;
3185:   }
3186:   VecRestoreArray(v, &vArray);
3187:   if (!*values) {
3188:     if (csize) *csize = size;
3189:     *values = array;
3190:   } else {
3191:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3192:     *csize = size;
3193:   }
3194:   return(0);
3195: }

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

3205:   *size = 0;
3206:   for (p = 0; p < numPoints*2; p += 2) {
3207:     const PetscInt point = points[p];
3208:     const PetscInt o     = points[p+1];
3209:     PetscInt       dof, off, d;
3210:     const PetscScalar *varr;

3212:     PetscSectionGetDof(section, point, &dof);
3213:     PetscSectionGetOffset(section, point, &off);
3214:     varr = &vArray[off];
3215:     if (o >= 0) {
3216:       for (d = 0; d < dof; ++d, ++offset)    array[offset] = varr[d];
3217:     } else {
3218:       for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
3219:     }
3220:   }
3221:   *size = offset;
3222:   return(0);
3223: }

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

3233:   *size = 0;
3234:   for (f = 0; f < numFields; ++f) {
3235:     PetscInt fcomp, p;

3237:     PetscSectionGetFieldComponents(section, f, &fcomp);
3238:     for (p = 0; p < numPoints*2; p += 2) {
3239:       const PetscInt point = points[p];
3240:       const PetscInt o     = points[p+1];
3241:       PetscInt       fdof, foff, d, c;
3242:       const PetscScalar *varr;

3244:       PetscSectionGetFieldDof(section, point, f, &fdof);
3245:       PetscSectionGetFieldOffset(section, point, f, &foff);
3246:       varr = &vArray[foff];
3247:       if (o >= 0) {
3248:         for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
3249:       } else {
3250:         for (d = fdof/fcomp-1; d >= 0; --d) {
3251:           for (c = 0; c < fcomp; ++c, ++offset) {
3252:             array[offset] = varr[d*fcomp+c];
3253:           }
3254:         }
3255:       }
3256:     }
3257:   }
3258:   *size = offset;
3259:   return(0);
3260: }

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

3267:   Not collective

3269:   Input Parameters:
3270: + dm - The DM
3271: . section - The section describing the layout in v, or NULL to use the default section
3272: . v - The local vector
3273: - point - The sieve point in the DM

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

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

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

3285:   Level: intermediate

3287: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3288: @*/
3289: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3290: {
3291:   PetscSection    clSection;
3292:   IS              clPoints;
3293:   PetscScalar    *array, *vArray;
3294:   PetscInt       *points = NULL;
3295:   const PetscInt *clp;
3296:   PetscInt        depth, numFields, numPoints, size;
3297:   PetscErrorCode  ierr;

3301:   if (!section) {DMGetDefaultSection(dm, &section);}
3304:   DMPlexGetDepth(dm, &depth);
3305:   PetscSectionGetNumFields(section, &numFields);
3306:   if (depth == 1 && numFields < 2) {
3307:     DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
3308:     return(0);
3309:   }
3310:   /* Get points */
3311:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3312:   if (!clPoints) {
3313:     PetscInt pStart, pEnd, p, q;

3315:     PetscSectionGetChart(section, &pStart, &pEnd);
3316:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3317:     /* Compress out points not in the section */
3318:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3319:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3320:         points[q*2]   = points[p];
3321:         points[q*2+1] = points[p+1];
3322:         ++q;
3323:       }
3324:     }
3325:     numPoints = q;
3326:   } else {
3327:     PetscInt dof, off;

3329:     PetscSectionGetDof(clSection, point, &dof);
3330:     PetscSectionGetOffset(clSection, point, &off);
3331:     ISGetIndices(clPoints, &clp);
3332:     numPoints = dof/2;
3333:     points    = (PetscInt *) &clp[off];
3334:   }
3335:   /* Get array */
3336:   if (!values || !*values) {
3337:     PetscInt asize = 0, dof, p;

3339:     for (p = 0; p < numPoints*2; p += 2) {
3340:       PetscSectionGetDof(section, points[p], &dof);
3341:       asize += dof;
3342:     }
3343:     if (!values) {
3344:       if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3345:       else           {ISRestoreIndices(clPoints, &clp);}
3346:       if (csize) *csize = asize;
3347:       return(0);
3348:     }
3349:     DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
3350:   } else {
3351:     array = *values;
3352:   }
3353:   VecGetArray(v, &vArray);
3354:   /* Get values */
3355:   if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);}
3356:   else               {DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);}
3357:   /* Cleanup points */
3358:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3359:   else           {ISRestoreIndices(clPoints, &clp);}
3360:   /* Cleanup array */
3361:   VecRestoreArray(v, &vArray);
3362:   if (!*values) {
3363:     if (csize) *csize = size;
3364:     *values = array;
3365:   } else {
3366:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3367:     *csize = size;
3368:   }
3369:   return(0);
3370: }

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

3377:   Not collective

3379:   Input Parameters:
3380: + dm - The DM
3381: . section - The section describing the layout in v, or NULL to use the default section
3382: . v - The local vector
3383: . point - The sieve point in the DM
3384: . csize - The number of values in the closure, or NULL
3385: - values - The array of values, which is a borrowed array and should not be freed

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

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

3393:   Level: intermediate

3395: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3396: @*/
3397: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3398: {
3399:   PetscInt       size = 0;

3403:   /* Should work without recalculating size */
3404:   DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
3405:   return(0);
3406: }

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

3413: 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[])
3414: {
3415:   PetscInt        cdof;   /* The number of constraints on this point */
3416:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3417:   PetscScalar    *a;
3418:   PetscInt        off, cind = 0, k;
3419:   PetscErrorCode  ierr;

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

3454: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3455: {
3456:   PetscInt        cdof;   /* The number of constraints on this point */
3457:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3458:   PetscScalar    *a;
3459:   PetscInt        off, cind = 0, k;
3460:   PetscErrorCode  ierr;

3463:   PetscSectionGetConstraintDof(section, point, &cdof);
3464:   PetscSectionGetOffset(section, point, &off);
3465:   a    = &array[off];
3466:   if (cdof) {
3467:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3468:     if (orientation >= 0) {
3469:       for (k = 0; k < dof; ++k) {
3470:         if ((cind < cdof) && (k == cdofs[cind])) {
3471:           fuse(&a[k], values[k]);
3472:           ++cind;
3473:         }
3474:       }
3475:     } else {
3476:       for (k = 0; k < dof; ++k) {
3477:         if ((cind < cdof) && (k == cdofs[cind])) {
3478:           fuse(&a[k], values[dof-k-1]);
3479:           ++cind;
3480:         }
3481:       }
3482:     }
3483:   }
3484:   return(0);
3485: }

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

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

3534: 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[])
3535: {
3536:   PetscScalar    *a;
3537:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3538:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3539:   PetscInt        cind = 0, k, c;
3540:   PetscErrorCode  ierr;

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

3573: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3574: {
3575:   PetscScalar    *array;
3576:   const PetscInt *cone, *coneO;
3577:   PetscInt        pStart, pEnd, p, numPoints, off, dof;
3578:   PetscErrorCode  ierr;

3581:   PetscSectionGetChart(section, &pStart, &pEnd);
3582:   DMPlexGetConeSize(dm, point, &numPoints);
3583:   DMPlexGetCone(dm, point, &cone);
3584:   DMPlexGetConeOrientation(dm, point, &coneO);
3585:   VecGetArray(v, &array);
3586:   for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
3587:     const PetscInt cp = !p ? point : cone[p-1];
3588:     const PetscInt o  = !p ? 0     : coneO[p-1];

3590:     if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
3591:     PetscSectionGetDof(section, cp, &dof);
3592:     /* ADD_VALUES */
3593:     {
3594:       const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3595:       PetscScalar    *a;
3596:       PetscInt        cdof, coff, cind = 0, k;

3598:       PetscSectionGetConstraintDof(section, cp, &cdof);
3599:       PetscSectionGetOffset(section, cp, &coff);
3600:       a    = &array[coff];
3601:       if (!cdof) {
3602:         if (o >= 0) {
3603:           for (k = 0; k < dof; ++k) {
3604:             a[k] += values[off+k];
3605:           }
3606:         } else {
3607:           for (k = 0; k < dof; ++k) {
3608:             a[k] += values[off+dof-k-1];
3609:           }
3610:         }
3611:       } else {
3612:         PetscSectionGetConstraintIndices(section, cp, &cdofs);
3613:         if (o >= 0) {
3614:           for (k = 0; k < dof; ++k) {
3615:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3616:             a[k] += values[off+k];
3617:           }
3618:         } else {
3619:           for (k = 0; k < dof; ++k) {
3620:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3621:             a[k] += values[off+dof-k-1];
3622:           }
3623:         }
3624:       }
3625:     }
3626:   }
3627:   VecRestoreArray(v, &array);
3628:   return(0);
3629: }

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

3636:   Not collective

3638:   Input Parameters:
3639: + dm - The DM
3640: . section - The section describing the layout in v, or NULL to use the default section
3641: . v - The local vector
3642: . point - The sieve point in the DM
3643: . values - The array of values
3644: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

3649:   Level: intermediate

3651: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
3652: @*/
3653: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3654: {
3655:   PetscSection    clSection;
3656:   IS              clPoints;
3657:   PetscScalar    *array;
3658:   PetscInt       *points = NULL;
3659:   const PetscInt *clp;
3660:   PetscInt        depth, numFields, numPoints, p;
3661:   PetscErrorCode  ierr;

3665:   if (!section) {DMGetDefaultSection(dm, &section);}
3668:   DMPlexGetDepth(dm, &depth);
3669:   PetscSectionGetNumFields(section, &numFields);
3670:   if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
3671:     DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
3672:     return(0);
3673:   }
3674:   /* Get points */
3675:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3676:   if (!clPoints) {
3677:     PetscInt pStart, pEnd, q;

3679:     PetscSectionGetChart(section, &pStart, &pEnd);
3680:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3681:     /* Compress out points not in the section */
3682:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3683:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3684:         points[q*2]   = points[p];
3685:         points[q*2+1] = points[p+1];
3686:         ++q;
3687:       }
3688:     }
3689:     numPoints = q;
3690:   } else {
3691:     PetscInt dof, off;

3693:     PetscSectionGetDof(clSection, point, &dof);
3694:     PetscSectionGetOffset(clSection, point, &off);
3695:     ISGetIndices(clPoints, &clp);
3696:     numPoints = dof/2;
3697:     points    = (PetscInt *) &clp[off];
3698:   }
3699:   /* Get array */
3700:   VecGetArray(v, &array);
3701:   /* Get values */
3702:   if (numFields > 0) {
3703:     PetscInt offset = 0, fcomp, f;
3704:     for (f = 0; f < numFields; ++f) {
3705:       PetscSectionGetFieldComponents(section, f, &fcomp);
3706:       switch (mode) {
3707:       case INSERT_VALUES:
3708:         for (p = 0; p < numPoints*2; p += 2) {
3709:           const PetscInt point = points[p];
3710:           const PetscInt o     = points[p+1];
3711:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3712:         } break;
3713:       case INSERT_ALL_VALUES:
3714:         for (p = 0; p < numPoints*2; p += 2) {
3715:           const PetscInt point = points[p];
3716:           const PetscInt o     = points[p+1];
3717:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3718:         } break;
3719:       case INSERT_BC_VALUES:
3720:         for (p = 0; p < numPoints*2; p += 2) {
3721:           const PetscInt point = points[p];
3722:           const PetscInt o     = points[p+1];
3723:           updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3724:         } break;
3725:       case ADD_VALUES:
3726:         for (p = 0; p < numPoints*2; p += 2) {
3727:           const PetscInt point = points[p];
3728:           const PetscInt o     = points[p+1];
3729:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3730:         } break;
3731:       case ADD_ALL_VALUES:
3732:         for (p = 0; p < numPoints*2; p += 2) {
3733:           const PetscInt point = points[p];
3734:           const PetscInt o     = points[p+1];
3735:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3736:         } break;
3737:       default:
3738:         SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3739:       }
3740:     }
3741:   } else {
3742:     PetscInt dof, off;

3744:     switch (mode) {
3745:     case INSERT_VALUES:
3746:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3747:         PetscInt o = points[p+1];
3748:         PetscSectionGetDof(section, points[p], &dof);
3749:         updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
3750:       } break;
3751:     case INSERT_ALL_VALUES:
3752:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3753:         PetscInt o = points[p+1];
3754:         PetscSectionGetDof(section, points[p], &dof);
3755:         updatePoint_private(section, points[p], dof, insert, PETSC_TRUE,  o, &values[off], array);
3756:       } break;
3757:     case INSERT_BC_VALUES:
3758:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3759:         PetscInt o = points[p+1];
3760:         PetscSectionGetDof(section, points[p], &dof);
3761:         updatePointBC_private(section, points[p], dof, insert,  o, &values[off], array);
3762:       } break;
3763:     case ADD_VALUES:
3764:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3765:         PetscInt o = points[p+1];
3766:         PetscSectionGetDof(section, points[p], &dof);
3767:         updatePoint_private(section, points[p], dof, add,    PETSC_FALSE, o, &values[off], array);
3768:       } break;
3769:     case ADD_ALL_VALUES:
3770:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3771:         PetscInt o = points[p+1];
3772:         PetscSectionGetDof(section, points[p], &dof);
3773:         updatePoint_private(section, points[p], dof, add,    PETSC_TRUE,  o, &values[off], array);
3774:       } break;
3775:     default:
3776:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3777:     }
3778:   }
3779:   /* Cleanup points */
3780:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3781:   else           {ISRestoreIndices(clPoints, &clp);}
3782:   /* Cleanup array */
3783:   VecRestoreArray(v, &array);
3784:   return(0);
3785: }

3789: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
3790: {
3791:   PetscSection    clSection;
3792:   IS              clPoints;
3793:   PetscScalar    *array;
3794:   PetscInt       *points = NULL;
3795:   const PetscInt *clp;
3796:   PetscInt        numFields, numPoints, p;
3797:   PetscInt        offset = 0, fcomp, f;
3798:   PetscErrorCode  ierr;

3802:   if (!section) {DMGetDefaultSection(dm, &section);}
3805:   PetscSectionGetNumFields(section, &numFields);
3806:   /* Get points */
3807:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3808:   if (!clPoints) {
3809:     PetscInt pStart, pEnd, q;

3811:     PetscSectionGetChart(section, &pStart, &pEnd);
3812:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3813:     /* Compress out points not in the section */
3814:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3815:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3816:         points[q*2]   = points[p];
3817:         points[q*2+1] = points[p+1];
3818:         ++q;
3819:       }
3820:     }
3821:     numPoints = q;
3822:   } else {
3823:     PetscInt dof, off;

3825:     PetscSectionGetDof(clSection, point, &dof);
3826:     PetscSectionGetOffset(clSection, point, &off);
3827:     ISGetIndices(clPoints, &clp);
3828:     numPoints = dof/2;
3829:     points    = (PetscInt *) &clp[off];
3830:   }
3831:   /* Get array */
3832:   VecGetArray(v, &array);
3833:   /* Get values */
3834:   for (f = 0; f < numFields; ++f) {
3835:     PetscSectionGetFieldComponents(section, f, &fcomp);
3836:     if (!fieldActive[f]) {
3837:       for (p = 0; p < numPoints*2; p += 2) {
3838:         PetscInt fdof;
3839:         PetscSectionGetFieldDof(section, points[p], f, &fdof);
3840:         offset += fdof;
3841:       }
3842:       continue;
3843:     }
3844:     switch (mode) {
3845:     case INSERT_VALUES:
3846:       for (p = 0; p < numPoints*2; p += 2) {
3847:         const PetscInt point = points[p];
3848:         const PetscInt o     = points[p+1];
3849:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3850:       } break;
3851:     case INSERT_ALL_VALUES:
3852:       for (p = 0; p < numPoints*2; p += 2) {
3853:         const PetscInt point = points[p];
3854:         const PetscInt o     = points[p+1];
3855:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3856:         } break;
3857:     case INSERT_BC_VALUES:
3858:       for (p = 0; p < numPoints*2; p += 2) {
3859:         const PetscInt point = points[p];
3860:         const PetscInt o     = points[p+1];
3861:         updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3862:       } break;
3863:     case ADD_VALUES:
3864:       for (p = 0; p < numPoints*2; p += 2) {
3865:         const PetscInt point = points[p];
3866:         const PetscInt o     = points[p+1];
3867:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3868:       } break;
3869:     case ADD_ALL_VALUES:
3870:       for (p = 0; p < numPoints*2; p += 2) {
3871:         const PetscInt point = points[p];
3872:         const PetscInt o     = points[p+1];
3873:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3874:       } break;
3875:     default:
3876:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3877:     }
3878:   }
3879:   /* Cleanup points */
3880:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3881:   else           {ISRestoreIndices(clPoints, &clp);}
3882:   /* Cleanup array */
3883:   VecRestoreArray(v, &array);
3884:   return(0);
3885: }

3889: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
3890: {
3891:   PetscMPIInt    rank;
3892:   PetscInt       i, j;

3896:   MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
3897:   PetscViewerASCIIPrintf(viewer, "[%D]mat for sieve point %D\n", rank, point);
3898:   for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%D]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
3899:   for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%D]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
3900:   numCIndices = numCIndices ? numCIndices : numRIndices;
3901:   for (i = 0; i < numRIndices; i++) {
3902:     PetscViewerASCIIPrintf(viewer, "[%D]", rank);
3903:     for (j = 0; j < numCIndices; j++) {
3904: #if defined(PETSC_USE_COMPLEX)
3905:       PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
3906: #else
3907:       PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
3908: #endif
3909:     }
3910:     PetscViewerASCIIPrintf(viewer, "\n");
3911:   }
3912:   return(0);
3913: }

3917: /* . off - The global offset of this point */
3918: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
3919: {
3920:   PetscInt        dof;    /* The number of unknowns on this point */
3921:   PetscInt        cdof;   /* The number of constraints on this point */
3922:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3923:   PetscInt        cind = 0, k;
3924:   PetscErrorCode  ierr;

3927:   PetscSectionGetDof(section, point, &dof);
3928:   PetscSectionGetConstraintDof(section, point, &cdof);
3929:   if (!cdof || setBC) {
3930:     if (orientation >= 0) {
3931:       for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
3932:     } else {
3933:       for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
3934:     }
3935:   } else {
3936:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3937:     if (orientation >= 0) {
3938:       for (k = 0; k < dof; ++k) {
3939:         if ((cind < cdof) && (k == cdofs[cind])) {
3940:           /* Insert check for returning constrained indices */
3941:           indices[*loff+k] = -(off+k+1);
3942:           ++cind;
3943:         } else {
3944:           indices[*loff+k] = off+k-cind;
3945:         }
3946:       }
3947:     } else {
3948:       for (k = 0; k < dof; ++k) {
3949:         if ((cind < cdof) && (k == cdofs[cind])) {
3950:           /* Insert check for returning constrained indices */
3951:           indices[*loff+dof-k-1] = -(off+k+1);
3952:           ++cind;
3953:         } else {
3954:           indices[*loff+dof-k-1] = off+k-cind;
3955:         }
3956:       }
3957:     }
3958:   }
3959:   *loff += dof;
3960:   return(0);
3961: }

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

3972:   PetscSectionGetNumFields(section, &numFields);
3973:   for (f = 0, foff = 0; f < numFields; ++f) {
3974:     PetscInt        fdof, fcomp, cfdof;
3975:     const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3976:     PetscInt        cind = 0, k, c;

3978:     PetscSectionGetFieldComponents(section, f, &fcomp);
3979:     PetscSectionGetFieldDof(section, point, f, &fdof);
3980:     PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
3981:     if (!cfdof || setBC) {
3982:       if (orientation >= 0) {
3983:         for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
3984:       } else {
3985:         for (k = fdof/fcomp-1; k >= 0; --k) {
3986:           for (c = 0; c < fcomp; ++c) {
3987:             indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
3988:           }
3989:         }
3990:       }
3991:     } else {
3992:       PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3993:       if (orientation >= 0) {
3994:         for (k = 0; k < fdof; ++k) {
3995:           if ((cind < cfdof) && (k == fcdofs[cind])) {
3996:             indices[foffs[f]+k] = -(off+foff+k+1);
3997:             ++cind;
3998:           } else {
3999:             indices[foffs[f]+k] = off+foff+k-cind;
4000:           }
4001:         }
4002:       } else {
4003:         for (k = fdof/fcomp-1; k >= 0; --k) {
4004:           for (c = 0; c < fcomp; ++c) {
4005:             if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
4006:               indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
4007:               ++cind;
4008:             } else {
4009:               indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
4010:             }
4011:           }
4012:         }
4013:       }
4014:     }
4015:     foff     += fdof - cfdof;
4016:     foffs[f] += fdof;
4017:   }
4018:   return(0);
4019: }

4023: 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[])
4024: {
4025:   Mat             cMat;
4026:   PetscSection    aSec, cSec;
4027:   IS              aIS;
4028:   PetscInt        aStart = -1, aEnd = -1;
4029:   const PetscInt  *anchors;
4030:   PetscInt        numFields, f, p, q, newP = 0;
4031:   PetscInt        newNumPoints = 0, newNumIndices = 0;
4032:   PetscInt        *newPoints, *indices, *newIndices;
4033:   PetscInt        maxAnchor, maxDof;
4034:   PetscInt        newOffsets[32];
4035:   PetscInt        *pointMatOffsets[32];
4036:   PetscInt        *newPointOffsets[32];
4037:   PetscScalar     *pointMat[32];
4038:   PetscScalar     *newValues,*tmpValues;
4039:   PetscBool       anyConstrained = PETSC_FALSE;
4040:   PetscErrorCode  ierr;

4045:   PetscSectionGetNumFields(section, &numFields);

4047:   DMPlexGetAnchors(dm,&aSec,&aIS);
4048:   /* if there are point-to-point constraints */
4049:   if (aSec) {
4050:     PetscMemzero(newOffsets, 32 * sizeof(PetscInt));
4051:     ISGetIndices(aIS,&anchors);
4052:     PetscSectionGetChart(aSec,&aStart,&aEnd);
4053:     /* figure out how many points are going to be in the new element matrix
4054:      * (we allow double counting, because it's all just going to be summed
4055:      * into the global matrix anyway) */
4056:     for (p = 0; p < 2*numPoints; p+=2) {
4057:       PetscInt b    = points[p];
4058:       PetscInt bDof = 0;

4060:       if (b >= aStart && b < aEnd) {
4061:         PetscSectionGetDof(aSec,b,&bDof);
4062:       }
4063:       if (bDof) {
4064:         /* this point is constrained */
4065:         /* it is going to be replaced by its anchors */
4066:         PetscInt bOff, q;

4068:         anyConstrained = PETSC_TRUE;
4069:         newNumPoints  += bDof;
4070:         PetscSectionGetOffset(aSec,b,&bOff);
4071:         for (q = 0; q < bDof; q++) {
4072:           PetscInt a = anchors[bOff + q];
4073:           PetscInt aDof;

4075:           PetscSectionGetDof(section,a,&aDof);
4076:           newNumIndices += aDof;
4077:           for (f = 0; f < numFields; ++f) {
4078:             PetscInt fDof;

4080:             PetscSectionGetFieldDof(section, a, f, &fDof);
4081:             newOffsets[f+1] += fDof;
4082:           }
4083:         }
4084:       }
4085:       else {
4086:         /* this point is not constrained */
4087:         newNumPoints++;
4088:         PetscSectionGetDof(section,b,&bDof);
4089:         newNumIndices += bDof;
4090:         for (f = 0; f < numFields; ++f) {
4091:           PetscInt fDof;

4093:           PetscSectionGetFieldDof(section, b, f, &fDof);
4094:           newOffsets[f+1] += fDof;
4095:         }
4096:       }
4097:     }
4098:   }
4099:   if (!anyConstrained) {
4100:     *outNumPoints  = 0;
4101:     *outNumIndices = 0;
4102:     *outPoints     = NULL;
4103:     *outValues     = NULL;
4104:     if (aSec) {
4105:       ISRestoreIndices(aIS,&anchors);
4106:     }
4107:     return(0);
4108:   }

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

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

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

4116:   /* output arrays */
4117:   DMGetWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4118:   DMGetWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);

4120:   /* workspaces */
4121:   DMGetWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4122:   if (numFields) {
4123:     for (f = 0; f < numFields; f++) {
4124:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4125:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4126:     }
4127:   }
4128:   else {
4129:     DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4130:     DMGetWorkArray(dm,numPoints,PETSC_INT,&newPointOffsets[0]);
4131:   }

4133:   /* get workspaces for the point-to-point matrices */
4134:   if (numFields) {
4135:     for (p = 0; p < numPoints; p++) {
4136:       PetscInt b    = points[2*p];
4137:       PetscInt bDof = 0;

4139:       if (b >= aStart && b < aEnd) {
4140:         PetscSectionGetDof(aSec, b, &bDof);
4141:       }
4142:       if (bDof) {
4143:         for (f = 0; f < numFields; f++) {
4144:           PetscInt fDof, q, bOff, allFDof = 0;

4146:           PetscSectionGetFieldDof(section, b, f, &fDof);
4147:           PetscSectionGetOffset(aSec, b, &bOff);
4148:           for (q = 0; q < bDof; q++) {
4149:             PetscInt a = anchors[bOff + q];
4150:             PetscInt aFDof;

4152:             PetscSectionGetFieldDof(section, a, f, &aFDof);
4153:             allFDof += aFDof;
4154:           }
4155:           newPointOffsets[f][p+1] = allFDof;
4156:           pointMatOffsets[f][p+1] = fDof * allFDof;
4157:         }
4158:       }
4159:       else {
4160:         for (f = 0; f < numFields; f++) {
4161:           PetscInt fDof;

4163:           PetscSectionGetFieldDof(section, b, f, &fDof);
4164:           newPointOffsets[f][p+1] = fDof;
4165:           pointMatOffsets[f][p+1] = 0;
4166:         }
4167:       }
4168:     }
4169:     for (f = 0; f < numFields; f++) {
4170:       newPointOffsets[f][0] = 0;
4171:       pointMatOffsets[f][0] = 0;
4172:       for (p = 0; p < numPoints; p++) {
4173:         newPointOffsets[f][p+1] += newPointOffsets[f][p];
4174:         pointMatOffsets[f][p+1] += pointMatOffsets[f][p];
4175:       }
4176:       DMGetWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4177:     }
4178:   }
4179:   else {
4180:     for (p = 0; p < numPoints; p++) {
4181:       PetscInt b    = points[2*p];
4182:       PetscInt bDof = 0;

4184:       if (b >= aStart && b < aEnd) {
4185:         PetscSectionGetDof(aSec, b, &bDof);
4186:       }
4187:       if (bDof) {
4188:         PetscInt dof, bOff, q, allDof = 0;

4190:         PetscSectionGetDof(section, b, &dof);
4191:         PetscSectionGetOffset(aSec, b, &bOff);
4192:         for (q = 0; q < bDof; q++) {
4193:           PetscInt a = anchors[bOff + q], aDof;

4195:           PetscSectionGetDof(section, a, &aDof);
4196:           allDof += aDof;
4197:         }
4198:         newPointOffsets[0][p+1] = allDof;
4199:         pointMatOffsets[0][p+1] = dof * allDof;
4200:       }
4201:       else {
4202:         PetscInt dof;

4204:         PetscSectionGetDof(section, b, &dof);
4205:         newPointOffsets[0][p+1] = dof;
4206:         pointMatOffsets[0][p+1] = 0;
4207:       }
4208:     }
4209:     newPointOffsets[0][0] = 0;
4210:     pointMatOffsets[0][0] = 0;
4211:     for (p = 0; p < numPoints; p++) {
4212:       newPointOffsets[0][p+1] += newPointOffsets[0][p];
4213:       pointMatOffsets[0][p+1] += pointMatOffsets[0][p];
4214:     }
4215:     DMGetWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4216:   }

4218:   /* get the point-to-point matrices; construct newPoints */
4219:   PetscSectionGetMaxDof(aSec, &maxAnchor);
4220:   PetscSectionGetMaxDof(section, &maxDof);
4221:   DMGetWorkArray(dm,maxDof,PETSC_INT,&indices);
4222:   DMGetWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4223:   if (numFields) {
4224:     for (p = 0, newP = 0; p < numPoints; p++) {
4225:       PetscInt b    = points[2*p];
4226:       PetscInt o    = points[2*p+1];
4227:       PetscInt bDof = 0;

4229:       if (b >= aStart && b < aEnd) {
4230:         PetscSectionGetDof(aSec, b, &bDof);
4231:       }
4232:       if (bDof) {
4233:         PetscInt fStart[32], fEnd[32], fAnchorStart[32], fAnchorEnd[32], bOff, q;

4235:         fStart[0] = 0;
4236:         fEnd[0]   = 0;
4237:         for (f = 0; f < numFields; f++) {
4238:           PetscInt fDof;

4240:           PetscSectionGetFieldDof(cSec, b, f, &fDof);
4241:           fStart[f+1] = fStart[f] + fDof;
4242:           fEnd[f+1]   = fStart[f+1];
4243:         }
4244:         PetscSectionGetOffset(cSec, b, &bOff);
4245:         indicesPointFields_private(cSec, b, bOff, fEnd, PETSC_TRUE, o, indices);

4247:         fAnchorStart[0] = 0;
4248:         fAnchorEnd[0]   = 0;
4249:         for (f = 0; f < numFields; f++) {
4250:           PetscInt fDof = newPointOffsets[f][p + 1] - newPointOffsets[f][p];

4252:           fAnchorStart[f+1] = fAnchorStart[f] + fDof;
4253:           fAnchorEnd[f+1]   = fAnchorStart[f + 1];
4254:         }
4255:         PetscSectionGetOffset (aSec, b, &bOff);
4256:         for (q = 0; q < bDof; q++) {
4257:           PetscInt a = anchors[bOff + q], aOff;

4259:           /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4260:           newPoints[2*(newP + q)]     = a;
4261:           newPoints[2*(newP + q) + 1] = 0;
4262:           PetscSectionGetOffset(section, a, &aOff);
4263:           indicesPointFields_private(section, a, aOff, fAnchorEnd, PETSC_TRUE, 0, newIndices);
4264:         }
4265:         newP += bDof;

4267:         /* get the point-to-point submatrix */
4268:         for (f = 0; f < numFields; f++) {
4269:           MatGetValues(cMat,fEnd[f]-fStart[f],indices + fStart[f],fAnchorEnd[f] - fAnchorStart[f],newIndices + fAnchorStart[f],pointMat[f] + pointMatOffsets[f][p]);
4270:         }
4271:       }
4272:       else {
4273:         newPoints[2 * newP]     = b;
4274:         newPoints[2 * newP + 1] = o;
4275:         newP++;
4276:       }
4277:     }
4278:   } else {
4279:     for (p = 0; p < numPoints; p++) {
4280:       PetscInt b    = points[2*p];
4281:       PetscInt o    = points[2*p+1];
4282:       PetscInt bDof = 0;

4284:       if (b >= aStart && b < aEnd) {
4285:         PetscSectionGetDof(aSec, b, &bDof);
4286:       }
4287:       if (bDof) {
4288:         PetscInt bEnd = 0, bAnchorEnd = 0, bOff;

4290:         PetscSectionGetOffset(cSec, b, &bOff);
4291:         indicesPoint_private(cSec, b, bOff, &bEnd, PETSC_TRUE, o, indices);

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

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

4299:           newPoints[2*(newP + q)]     = a;
4300:           newPoints[2*(newP + q) + 1] = 0;
4301:           PetscSectionGetOffset(section, a, &aOff);
4302:           indicesPoint_private(section, a, aOff, &bAnchorEnd, PETSC_TRUE, 0, newIndices);
4303:         }
4304:         newP += bDof;

4306:         /* get the point-to-point submatrix */
4307:         MatGetValues(cMat,bEnd,indices,bAnchorEnd,newIndices,pointMat[0] + pointMatOffsets[0][p]);
4308:       }
4309:       else {
4310:         newPoints[2 * newP]     = b;
4311:         newPoints[2 * newP + 1] = o;
4312:         newP++;
4313:       }
4314:     }
4315:   }

4317:   PetscMemzero(tmpValues,newNumIndices*numIndices*sizeof(*tmpValues));
4318:   /* multiply constraints on the right */
4319:   if (numFields) {
4320:     for (f = 0; f < numFields; f++) {
4321:       PetscInt oldOff = offsets[f];

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

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

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

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

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

4387:   PetscMemzero(newValues,newNumIndices*newNumIndices*sizeof(*newValues));
4388:   /* multiply constraints transpose on the left */
4389:   if (numFields) {
4390:     for (f = 0; f < numFields; f++) {
4391:       PetscInt oldOff = offsets[f];

4393:       for (p = 0; p < numPoints; p++) {
4394:         PetscInt rStart = newPointOffsets[f][p];
4395:         PetscInt b      = points[2 * p];
4396:         PetscInt c, r, k;
4397:         PetscInt dof;

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

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

4427:     for (p = 0; p < numPoints; p++) {
4428:       PetscInt rStart = newPointOffsets[0][p];
4429:       PetscInt b      = points[2 * p];
4430:       PetscInt c, r, k;
4431:       PetscInt dof;

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

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

4458:   /* clean up */
4459:   DMRestoreWorkArray(dm,maxDof,PETSC_INT,&indices);
4460:   DMRestoreWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4461:   if (numFields) {
4462:     for (f = 0; f < numFields; f++) {
4463:       DMRestoreWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4464:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4465:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4466:     }
4467:   }
4468:   else {
4469:     DMRestoreWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4470:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4471:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[0]);
4472:   }
4473:   DMRestoreWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4474:   ISRestoreIndices(aIS,&anchors);

4476:   /* output */
4477:   *outNumPoints  = newNumPoints;
4478:   *outNumIndices = newNumIndices;
4479:   *outPoints     = newPoints;
4480:   *outValues     = newValues;
4481:   for (f = 0; f < numFields; f++) {
4482:     offsets[f] = newOffsets[f];
4483:   }
4484:   return(0);
4485: }

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

4492:   Not collective

4494:   Input Parameters:
4495: + dm - The DM
4496: . section - The section describing the layout in v, or NULL to use the default section
4497: . globalSection - The section describing the layout in v, or NULL to use the default global section
4498: . A - The matrix
4499: . point - The sieve point in the DM
4500: . values - The array of values
4501: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

4506:   Level: intermediate

4508: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
4509: @*/
4510: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4511: {
4512:   DM_Plex        *mesh   = (DM_Plex*) dm->data;
4513:   PetscSection    clSection;
4514:   IS              clPoints;
4515:   PetscInt       *points = NULL, *newPoints;
4516:   const PetscInt *clp;
4517:   PetscInt       *indices;
4518:   PetscInt        offsets[32];
4519:   PetscInt        numFields, numPoints, newNumPoints, numIndices, newNumIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
4520:   PetscScalar    *newValues;
4521:   PetscErrorCode  ierr;

4525:   if (!section) {DMGetDefaultSection(dm, &section);}
4527:   if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
4530:   PetscSectionGetNumFields(section, &numFields);
4531:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4532:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
4533:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4534:   if (!clPoints) {
4535:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4536:     /* Compress out points not in the section */
4537:     PetscSectionGetChart(section, &pStart, &pEnd);
4538:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4539:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4540:         points[q*2]   = points[p];
4541:         points[q*2+1] = points[p+1];
4542:         ++q;
4543:       }
4544:     }
4545:     numPoints = q;
4546:   } else {
4547:     PetscInt dof, off;

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

4558:     PetscSectionGetDof(section, points[p], &dof);
4559:     for (f = 0; f < numFields; ++f) {
4560:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
4561:       offsets[f+1] += fdof;
4562:     }
4563:     numIndices += dof;
4564:   }
4565:   for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];

4567:   if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
4568:   DMPlexAnchorsModifyMat(dm,section,numPoints,numIndices,points,values,&newNumPoints,&newNumIndices,&newPoints,&newValues,offsets);
4569:   if (newNumPoints) {
4570:     if (!clPoints) {
4571:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4572:     } else {
4573:       ISRestoreIndices(clPoints, &clp);
4574:     }
4575:     numPoints  = newNumPoints;
4576:     numIndices = newNumIndices;
4577:     points     = newPoints;
4578:     values     = newValues;
4579:   }
4580:   DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
4581:   if (numFields) {
4582:     for (p = 0; p < numPoints*2; p += 2) {
4583:       PetscInt o = points[p+1];
4584:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4585:       indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
4586:     }
4587:   } else {
4588:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
4589:       PetscInt o = points[p+1];
4590:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4591:       indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
4592:     }
4593:   }
4594:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
4595:   MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
4596:   if (ierr) {
4597:     PetscMPIInt    rank;
4598:     PetscErrorCode ierr2;

4600:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4601:     ierr2 = (*PetscErrorPrintf)("[%D]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4602:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
4603:     ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
4604: 
4605:   }
4606:   if (newNumPoints) {
4607:     DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4608:     DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4609:   }
4610:   else {
4611:     if (!clPoints) {
4612:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4613:     } else {
4614:       ISRestoreIndices(clPoints, &clp);
4615:     }
4616:   }
4617:   DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
4618:   return(0);
4619: }

4623: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4624: {
4625:   DM_Plex        *mesh   = (DM_Plex*) dmf->data;
4626:   PetscInt       *fpoints = NULL, *ftotpoints = NULL;
4627:   PetscInt       *cpoints = NULL;
4628:   PetscInt       *findices, *cindices;
4629:   PetscInt        foffsets[32], coffsets[32];
4630:   CellRefiner     cellRefiner;
4631:   PetscInt        numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4632:   PetscErrorCode  ierr;

4637:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4639:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
4641:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4643:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4646:   PetscSectionGetNumFields(fsection, &numFields);
4647:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4648:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4649:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4650:   /* Column indices */
4651:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4652:   maxFPoints = numCPoints;
4653:   /* Compress out points not in the section */
4654:   /*   TODO: Squeeze out points with 0 dof as well */
4655:   PetscSectionGetChart(csection, &pStart, &pEnd);
4656:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4657:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4658:       cpoints[q*2]   = cpoints[p];
4659:       cpoints[q*2+1] = cpoints[p+1];
4660:       ++q;
4661:     }
4662:   }
4663:   numCPoints = q;
4664:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4665:     PetscInt fdof;

4667:     PetscSectionGetDof(csection, cpoints[p], &dof);
4668:     if (!dof) continue;
4669:     for (f = 0; f < numFields; ++f) {
4670:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4671:       coffsets[f+1] += fdof;
4672:     }
4673:     numCIndices += dof;
4674:   }
4675:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4676:   /* Row indices */
4677:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4678:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4679:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4680:   for (r = 0, q = 0; r < numSubcells; ++r) {
4681:     /* TODO Map from coarse to fine cells */
4682:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4683:     /* Compress out points not in the section */
4684:     PetscSectionGetChart(fsection, &pStart, &pEnd);
4685:     for (p = 0; p < numFPoints*2; p += 2) {
4686:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4687:         PetscSectionGetDof(fsection, fpoints[p], &dof);
4688:         if (!dof) continue;
4689:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4690:         if (s < q) continue;
4691:         ftotpoints[q*2]   = fpoints[p];
4692:         ftotpoints[q*2+1] = fpoints[p+1];
4693:         ++q;
4694:       }
4695:     }
4696:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4697:   }
4698:   numFPoints = q;
4699:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4700:     PetscInt fdof;

4702:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4703:     if (!dof) continue;
4704:     for (f = 0; f < numFields; ++f) {
4705:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4706:       foffsets[f+1] += fdof;
4707:     }
4708:     numFIndices += dof;
4709:   }
4710:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

4712:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4713:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4714:   DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4715:   DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4716:   if (numFields) {
4717:     for (p = 0; p < numFPoints*2; p += 2) {
4718:       PetscInt o = ftotpoints[p+1];
4719:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4720:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4721:     }
4722:     for (p = 0; p < numCPoints*2; p += 2) {
4723:       PetscInt o = cpoints[p+1];
4724:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4725:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4726:     }
4727:   } else {
4728:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4729:       PetscInt o = ftotpoints[p+1];
4730:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4731:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4732:     }
4733:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4734:       PetscInt o = cpoints[p+1];
4735:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4736:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4737:     }
4738:   }
4739:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
4740:   MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
4741:   if (ierr) {
4742:     PetscMPIInt    rank;
4743:     PetscErrorCode ierr2;

4745:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4746:     ierr2 = (*PetscErrorPrintf)("[%D]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4747:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
4748:     ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
4749:     ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
4750: 
4751:   }
4752:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4753:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4754:   DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4755:   DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4756:   return(0);
4757: }

4761: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
4762: {
4763:   PetscInt      *fpoints = NULL, *ftotpoints = NULL;
4764:   PetscInt      *cpoints = NULL;
4765:   PetscInt       foffsets[32], coffsets[32];
4766:   CellRefiner    cellRefiner;
4767:   PetscInt       numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;

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

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

4837:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4838:     if (!dof) continue;
4839:     for (f = 0; f < numFields; ++f) {
4840:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4841:       foffsets[f+1] += fdof;
4842:     }
4843:     numFIndices += dof;
4844:   }
4845:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

4847:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4848:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4849:   if (numFields) {
4850:     for (p = 0; p < numFPoints*2; p += 2) {
4851:       PetscInt o = ftotpoints[p+1];
4852:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4853:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4854:     }
4855:     for (p = 0; p < numCPoints*2; p += 2) {
4856:       PetscInt o = cpoints[p+1];
4857:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4858:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4859:     }
4860:   } else {
4861:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4862:       PetscInt o = ftotpoints[p+1];
4863:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4864:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4865:     }
4866:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4867:       PetscInt o = cpoints[p+1];
4868:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4869:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4870:     }
4871:   }
4872:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4873:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4874:   return(0);
4875: }

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

4882:   Input Parameter:
4883: . dm - The DMPlex object

4885:   Output Parameters:
4886: + cMax - The first hybrid cell
4887: . fMax - The first hybrid face
4888: . eMax - The first hybrid edge
4889: - vMax - The first hybrid vertex

4891:   Level: developer

4893: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
4894: @*/
4895: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
4896: {
4897:   DM_Plex       *mesh = (DM_Plex*) dm->data;
4898:   PetscInt       dim;

4903:   DMGetDimension(dm, &dim);
4904:   if (cMax) *cMax = mesh->hybridPointMax[dim];
4905:   if (fMax) *fMax = mesh->hybridPointMax[dim-1];
4906:   if (eMax) *eMax = mesh->hybridPointMax[1];
4907:   if (vMax) *vMax = mesh->hybridPointMax[0];
4908:   return(0);
4909: }

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

4916:   Input Parameters:
4917: . dm   - The DMPlex object
4918: . cMax - The first hybrid cell
4919: . fMax - The first hybrid face
4920: . eMax - The first hybrid edge
4921: - vMax - The first hybrid vertex

4923:   Level: developer

4925: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
4926: @*/
4927: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
4928: {
4929:   DM_Plex       *mesh = (DM_Plex*) dm->data;
4930:   PetscInt       dim;

4935:   DMGetDimension(dm, &dim);
4936:   if (cMax >= 0) mesh->hybridPointMax[dim]   = cMax;
4937:   if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
4938:   if (eMax >= 0) mesh->hybridPointMax[1]     = eMax;
4939:   if (vMax >= 0) mesh->hybridPointMax[0]     = vMax;
4940:   return(0);
4941: }

4945: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
4946: {
4947:   DM_Plex *mesh = (DM_Plex*) dm->data;

4952:   *cellHeight = mesh->vtkCellHeight;
4953:   return(0);
4954: }

4958: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
4959: {
4960:   DM_Plex *mesh = (DM_Plex*) dm->data;

4964:   mesh->vtkCellHeight = cellHeight;
4965:   return(0);
4966: }

4970: /* We can easily have a form that takes an IS instead */
4971: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
4972: {
4973:   PetscSection   section, globalSection;
4974:   PetscInt      *numbers, p;

4978:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
4979:   PetscSectionSetChart(section, pStart, pEnd);
4980:   for (p = pStart; p < pEnd; ++p) {
4981:     PetscSectionSetDof(section, p, 1);
4982:   }
4983:   PetscSectionSetUp(section);
4984:   PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, &globalSection);
4985:   PetscMalloc1(pEnd - pStart, &numbers);
4986:   for (p = pStart; p < pEnd; ++p) {
4987:     PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
4988:     if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
4989:     else                       numbers[p-pStart] += shift;
4990:   }
4991:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
4992:   if (globalSize) {
4993:     PetscLayout layout;
4994:     PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
4995:     PetscLayoutGetSize(layout, globalSize);
4996:     PetscLayoutDestroy(&layout);
4997:   }
4998:   PetscSectionDestroy(&section);
4999:   PetscSectionDestroy(&globalSection);
5000:   return(0);
5001: }

5005: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5006: {
5007:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5008:   PetscInt       cellHeight, cStart, cEnd, cMax;

5013:   if (!mesh->globalCellNumbers) {
5014:     DMPlexGetVTKCellHeight(dm, &cellHeight);
5015:     DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5016:     DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5017:     if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
5018:     DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, &mesh->globalCellNumbers);
5019:   }
5020:   *globalCellNumbers = mesh->globalCellNumbers;
5021:   return(0);
5022: }

5026: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5027: {
5028:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5029:   PetscInt       vStart, vEnd, vMax;

5034:   if (!mesh->globalVertexNumbers) {
5035:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5036:     DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5037:     if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
5038:     DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, &mesh->globalVertexNumbers);
5039:   }
5040:   *globalVertexNumbers = mesh->globalVertexNumbers;
5041:   return(0);
5042: }

5046: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5047: {
5048:   IS             nums[4];
5049:   PetscInt       depths[4];
5050:   PetscInt       depth, d, shift = 0;

5055:   DMPlexGetDepth(dm, &depth);
5056:   /* For unstratified meshes use dim instead of depth */
5057:   if (depth < 0) {DMGetDimension(dm, &depth);}
5058:   depths[0] = depth; depths[1] = 0;
5059:   for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5060:   for (d = 0; d <= depth; ++d) {
5061:     PetscInt pStart, pEnd, gsize;

5063:     DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5064:     DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5065:     shift += gsize;
5066:   }
5067:   ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5068:   for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5069:   return(0);
5070: }


5075: /*@C
5076:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
5077:   the local section and an SF describing the section point overlap.

5079:   Input Parameters:
5080:   + s - The PetscSection for the local field layout
5081:   . sf - The SF describing parallel layout of the section points
5082:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
5083:   . label - The label specifying the points
5084:   - labelValue - The label stratum specifying the points

5086:   Output Parameter:
5087:   . gsection - The PetscSection for the global field layout

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

5091:   Level: developer

5093: .seealso: PetscSectionCreate()
5094: @*/
5095: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
5096: {
5097:   PetscInt      *neg = NULL, *tmpOff = NULL;
5098:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

5102:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
5103:   PetscSectionGetChart(s, &pStart, &pEnd);
5104:   PetscSectionSetChart(*gsection, pStart, pEnd);
5105:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
5106:   if (nroots >= 0) {
5107:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
5108:     PetscCalloc1(nroots, &neg);
5109:     if (nroots > pEnd-pStart) {
5110:       PetscCalloc1(nroots, &tmpOff);
5111:     } else {
5112:       tmpOff = &(*gsection)->atlasDof[-pStart];
5113:     }
5114:   }
5115:   /* Mark ghost points with negative dof */
5116:   for (p = pStart; p < pEnd; ++p) {
5117:     PetscInt value;

5119:     DMLabelGetValue(label, p, &value);
5120:     if (value != labelValue) continue;
5121:     PetscSectionGetDof(s, p, &dof);
5122:     PetscSectionSetDof(*gsection, p, dof);
5123:     PetscSectionGetConstraintDof(s, p, &cdof);
5124:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
5125:     if (neg) neg[p] = -(dof+1);
5126:   }
5127:   PetscSectionSetUpBC(*gsection);
5128:   if (nroots >= 0) {
5129:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5130:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5131:     if (nroots > pEnd-pStart) {
5132:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
5133:     }
5134:   }
5135:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
5136:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5137:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
5138:     (*gsection)->atlasOff[p] = off;
5139:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
5140:   }
5141:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
5142:   globalOff -= off;
5143:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5144:     (*gsection)->atlasOff[p] += globalOff;
5145:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
5146:   }
5147:   /* Put in negative offsets for ghost points */
5148:   if (nroots >= 0) {
5149:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5150:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5151:     if (nroots > pEnd-pStart) {
5152:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
5153:     }
5154:   }
5155:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
5156:   PetscFree(neg);
5157:   return(0);
5158: }

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

5165:   Input Parameters:
5166:   + dm - The DMPlex object

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

5170:   Level: developer

5172: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5173: @*/
5174: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5175: {
5176:   PetscSection    coneSection, supportSection;
5177:   const PetscInt *cone, *support;
5178:   PetscInt        coneSize, c, supportSize, s;
5179:   PetscInt        pStart, pEnd, p, csize, ssize;
5180:   PetscErrorCode  ierr;

5184:   DMPlexGetConeSection(dm, &coneSection);
5185:   DMPlexGetSupportSection(dm, &supportSection);
5186:   /* Check that point p is found in the support of its cone points, and vice versa */
5187:   DMPlexGetChart(dm, &pStart, &pEnd);
5188:   for (p = pStart; p < pEnd; ++p) {
5189:     DMPlexGetConeSize(dm, p, &coneSize);
5190:     DMPlexGetCone(dm, p, &cone);
5191:     for (c = 0; c < coneSize; ++c) {
5192:       PetscBool dup = PETSC_FALSE;
5193:       PetscInt  d;
5194:       for (d = c-1; d >= 0; --d) {
5195:         if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5196:       }
5197:       DMPlexGetSupportSize(dm, cone[c], &supportSize);
5198:       DMPlexGetSupport(dm, cone[c], &support);
5199:       for (s = 0; s < supportSize; ++s) {
5200:         if (support[s] == p) break;
5201:       }
5202:       if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5203:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5204:         for (s = 0; s < coneSize; ++s) {
5205:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
5206:         }
5207:         PetscPrintf(PETSC_COMM_SELF, "\n");
5208:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
5209:         for (s = 0; s < supportSize; ++s) {
5210:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
5211:         }
5212:         PetscPrintf(PETSC_COMM_SELF, "\n");
5213:         if (dup) {
5214:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
5215:         } else {
5216:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
5217:         }
5218:       }
5219:     }
5220:     DMPlexGetSupportSize(dm, p, &supportSize);
5221:     DMPlexGetSupport(dm, p, &support);
5222:     for (s = 0; s < supportSize; ++s) {
5223:       DMPlexGetConeSize(dm, support[s], &coneSize);
5224:       DMPlexGetCone(dm, support[s], &cone);
5225:       for (c = 0; c < coneSize; ++c) {
5226:         if (cone[c] == p) break;
5227:       }
5228:       if (c >= coneSize) {
5229:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
5230:         for (c = 0; c < supportSize; ++c) {
5231:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
5232:         }
5233:         PetscPrintf(PETSC_COMM_SELF, "\n");
5234:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
5235:         for (c = 0; c < coneSize; ++c) {
5236:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
5237:         }
5238:         PetscPrintf(PETSC_COMM_SELF, "\n");
5239:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
5240:       }
5241:     }
5242:   }
5243:   PetscSectionGetStorageSize(coneSection, &csize);
5244:   PetscSectionGetStorageSize(supportSection, &ssize);
5245:   if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
5246:   return(0);
5247: }

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

5254:   Input Parameters:
5255: + dm - The DMPlex object
5256: . isSimplex - Are the cells simplices or tensor products
5257: - cellHeight - Normally 0

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

5261:   Level: developer

5263: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
5264: @*/
5265: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5266: {
5267:   PetscInt       dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;

5272:   DMGetDimension(dm, &dim);
5273:   switch (dim) {
5274:   case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
5275:   case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
5276:   case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
5277:   default:
5278:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
5279:   }
5280:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5281:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5282:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5283:   cMax = cMax >= 0 ? cMax : cEnd;
5284:   for (c = cStart; c < cMax; ++c) {
5285:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5287:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5288:     for (cl = 0; cl < closureSize*2; cl += 2) {
5289:       const PetscInt p = closure[cl];
5290:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5291:     }
5292:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5293:     if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has  %d vertices != %d", c, coneSize, numCorners);
5294:   }
5295:   for (c = cMax; c < cEnd; ++c) {
5296:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5298:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5299:     for (cl = 0; cl < closureSize*2; cl += 2) {
5300:       const PetscInt p = closure[cl];
5301:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5302:     }
5303:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5304:     if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has  %d vertices > %d", c, coneSize, numHybridCorners);
5305:   }
5306:   return(0);
5307: }

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

5314:   Input Parameters:
5315: + dm - The DMPlex object
5316: . isSimplex - Are the cells simplices or tensor products
5317: - cellHeight - Normally 0

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

5321:   Level: developer

5323: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
5324: @*/
5325: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5326: {
5327:   PetscInt       pMax[4];
5328:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, h;

5333:   DMGetDimension(dm, &dim);
5334:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5335:   DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
5336:   for (h = cellHeight; h < dim; ++h) {
5337:     DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
5338:     for (c = cStart; c < cEnd; ++c) {
5339:       const PetscInt *cone, *ornt, *faces;
5340:       PetscInt        numFaces, faceSize, coneSize,f;
5341:       PetscInt       *closure = NULL, closureSize, cl, numCorners = 0;

5343:       if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
5344:       DMPlexGetConeSize(dm, c, &coneSize);
5345:       DMPlexGetCone(dm, c, &cone);
5346:       DMPlexGetConeOrientation(dm, c, &ornt);
5347:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5348:       for (cl = 0; cl < closureSize*2; cl += 2) {
5349:         const PetscInt p = closure[cl];
5350:         if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
5351:       }
5352:       DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
5353:       if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
5354:       for (f = 0; f < numFaces; ++f) {
5355:         PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;

5357:         DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
5358:         for (cl = 0; cl < fclosureSize*2; cl += 2) {
5359:           const PetscInt p = fclosure[cl];
5360:           if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
5361:         }
5362:         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);
5363:         for (v = 0; v < fnumCorners; ++v) {
5364:           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]);
5365:         }
5366:         DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
5367:       }
5368:       DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
5369:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5370:     }
5371:   }
5372:   return(0);
5373: }

5377: /* Pointwise interpolation
5378:      Just code FEM for now
5379:      u^f = I u^c
5380:      sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
5381:      u^f_i = sum_j psi^f_i I phi^c_j u^c_j
5382:      I_{ij} = psi^f_i phi^c_j
5383: */
5384: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
5385: {
5386:   PetscSection   gsc, gsf;
5387:   PetscInt       m, n;
5388:   void          *ctx;

5392:   /*
5393:   Loop over coarse cells
5394:     Loop over coarse basis functions
5395:       Loop over fine cells in coarse cell
5396:         Loop over fine dual basis functions
5397:           Evaluate coarse basis on fine dual basis quad points
5398:           Sum
5399:           Update local element matrix
5400:     Accumulate to interpolation matrix

5402:    Can extend PetscFEIntegrateJacobian_Basic() to do a specialized cell loop
5403:   */
5404:   DMGetDefaultGlobalSection(dmFine, &gsf);
5405:   PetscSectionGetConstrainedStorageSize(gsf, &m);
5406:   DMGetDefaultGlobalSection(dmCoarse, &gsc);
5407:   PetscSectionGetConstrainedStorageSize(gsc, &n);
5408:   /* We need to preallocate properly */
5409:   MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
5410:   MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
5411:   MatSetType(*interpolation, dmCoarse->mattype);
5412:   DMGetApplicationContext(dmFine, &ctx);
5413:   DMPlexComputeInterpolatorFEM(dmCoarse, dmFine, *interpolation, ctx);
5414:   /* Use naive scaling */
5415:   DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
5416:   return(0);
5417: }

5421: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat)
5422: {
5424:   VecScatter     ctx;

5427:   DMPlexComputeInjectorFEM(dmCoarse, dmFine, &ctx, NULL);
5428:   MatCreateScatter(PetscObjectComm((PetscObject)ctx), ctx, mat);
5429:   VecScatterDestroy(&ctx);
5430:   return(0);
5431: }

5435: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
5436: {
5437:   PetscSection   section;
5438:   IS            *bcPoints;
5439:   PetscBool     *isFE;
5440:   PetscInt      *bcFields, *numComp, *numDof;
5441:   PetscInt       depth, dim, numBd, numBC = 0, numFields, bd, bc = 0, f;
5442:   PetscInt       cStart, cEnd, cEndInterior;

5446:   DMGetNumFields(dm, &numFields);
5447:   /* FE and FV boundary conditions are handled slightly differently */
5448:   PetscMalloc1(numFields, &isFE);
5449:   for (f = 0; f < numFields; ++f) {
5450:     PetscObject  obj;
5451:     PetscClassId id;

5453:     DMGetField(dm, f, &obj);
5454:     PetscObjectGetClassId(obj, &id);
5455:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
5456:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
5457:     else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
5458:   }
5459:   /* Allocate boundary point storage for FEM boundaries */
5460:   DMPlexGetDepth(dm, &depth);
5461:   DMGetDimension(dm, &dim);
5462:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
5463:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
5464:   DMPlexGetNumBoundary(dm, &numBd);
5465:   for (bd = 0; bd < numBd; ++bd) {
5466:     PetscInt  field;
5467:     PetscBool isEssential;

5469:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, NULL, &field, NULL, NULL, NULL, NULL);
5470:     if (isFE[field] && isEssential) ++numBC;
5471:   }
5472:   /* Add ghost cell boundaries for FVM */
5473:   for (f = 0; f < numFields; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
5474:   PetscMalloc2(numBC,&bcFields,numBC,&bcPoints);
5475:   /* Constrain ghost cells for FV */
5476:   for (f = 0; f < numFields; ++f) {
5477:     PetscInt *newidx, c;

5479:     if (isFE[f] || cEndInterior < 0) continue;
5480:     PetscMalloc1(cEnd-cEndInterior,&newidx);
5481:     for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
5482:     bcFields[bc] = f;
5483:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5484:   }
5485:   /* Handle FEM Dirichlet boundaries */
5486:   for (bd = 0; bd < numBd; ++bd) {
5487:     const char     *bdLabel;
5488:     DMLabel         label;
5489:     const PetscInt *values;
5490:     PetscInt        bd2, field, numValues;
5491:     PetscBool       isEssential, duplicate = PETSC_FALSE;

5493:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);
5494:     if (!isFE[field]) continue;
5495:     DMPlexGetLabel(dm, bdLabel, &label);
5496:     /* Only want to modify label once */
5497:     for (bd2 = 0; bd2 < bd; ++bd2) {
5498:       const char *bdname;
5499:       DMPlexGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL);
5500:       PetscStrcmp(bdname, bdLabel, &duplicate);
5501:       if (duplicate) break;
5502:     }
5503:     if (!duplicate && (isFE[field])) {
5504:       DMPlexLabelComplete(dm, label);
5505:       DMPlexLabelAddCells(dm, label);
5506:     }
5507:     /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
5508:     if (isEssential) {
5509:       PetscInt       *newidx;
5510:       PetscInt        n, newn = 0, p, v;

5512:       bcFields[bc] = field;
5513:       for (v = 0; v < numValues; ++v) {
5514:         IS              tmp;
5515:         const PetscInt *idx;

5517:         DMPlexGetStratumIS(dm, bdLabel, values[v], &tmp);
5518:         if (!tmp) continue;
5519:         ISGetLocalSize(tmp, &n);
5520:         ISGetIndices(tmp, &idx);
5521:         if (isFE[field]) {
5522:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
5523:         } else {
5524:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
5525:         }
5526:         ISRestoreIndices(tmp, &idx);
5527:         ISDestroy(&tmp);
5528:       }
5529:       PetscMalloc1(newn,&newidx);
5530:       newn = 0;
5531:       for (v = 0; v < numValues; ++v) {
5532:         IS              tmp;
5533:         const PetscInt *idx;

5535:         DMPlexGetStratumIS(dm, bdLabel, values[v], &tmp);
5536:         if (!tmp) continue;
5537:         ISGetLocalSize(tmp, &n);
5538:         ISGetIndices(tmp, &idx);
5539:         if (isFE[field]) {
5540:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
5541:         } else {
5542:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
5543:         }
5544:         ISRestoreIndices(tmp, &idx);
5545:         ISDestroy(&tmp);
5546:       }
5547:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5548:     }
5549:   }
5550:   /* Handle discretization */
5551:   PetscCalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
5552:   for (f = 0; f < numFields; ++f) {
5553:     PetscObject obj;

5555:     DMGetField(dm, f, &obj);
5556:     if (isFE[f]) {
5557:       PetscFE         fe = (PetscFE) obj;
5558:       const PetscInt *numFieldDof;
5559:       PetscInt        d;

5561:       PetscFEGetNumComponents(fe, &numComp[f]);
5562:       PetscFEGetNumDof(fe, &numFieldDof);
5563:       for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
5564:     } else {
5565:       PetscFV fv = (PetscFV) obj;

5567:       PetscFVGetNumComponents(fv, &numComp[f]);
5568:       numDof[f*(dim+1)+dim] = numComp[f];
5569:     }
5570:   }
5571:   for (f = 0; f < numFields; ++f) {
5572:     PetscInt d;
5573:     for (d = 1; d < dim; ++d) {
5574:       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.");
5575:     }
5576:   }
5577:   DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, NULL, &section);
5578:   for (f = 0; f < numFields; ++f) {
5579:     PetscFE     fe;
5580:     const char *name;

5582:     DMGetField(dm, f, (PetscObject *) &fe);
5583:     PetscObjectGetName((PetscObject) fe, &name);
5584:     PetscSectionSetFieldName(section, f, name);
5585:   }
5586:   DMSetDefaultSection(dm, section);
5587:   PetscSectionDestroy(&section);
5588:   for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);}
5589:   PetscFree2(bcFields,bcPoints);
5590:   PetscFree2(numComp,numDof);
5591:   PetscFree(isFE);
5592:   return(0);
5593: }

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

5600:   Input Parameter:
5601: . dm - The DMPlex object

5603:   Output Parameter:
5604: . cdm - The coarse DM

5606:   Level: intermediate

5608: .seealso: DMPlexSetCoarseDM()
5609: @*/
5610: PetscErrorCode DMPlexGetCoarseDM(DM dm, DM *cdm)
5611: {
5615:   *cdm = ((DM_Plex *) dm->data)->coarseMesh;
5616:   return(0);
5617: }

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

5624:   Input Parameters:
5625: + dm - The DMPlex object
5626: - cdm - The coarse DM

5628:   Level: intermediate

5630: .seealso: DMPlexGetCoarseDM()
5631: @*/
5632: PetscErrorCode DMPlexSetCoarseDM(DM dm, DM cdm)
5633: {
5634:   DM_Plex       *mesh;

5640:   mesh = (DM_Plex *) dm->data;
5641:   DMDestroy(&mesh->coarseMesh);
5642:   mesh->coarseMesh = cdm;
5643:   PetscObjectReference((PetscObject) mesh->coarseMesh);
5644:   return(0);
5645: }

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

5654:   not collective

5656:   Input Parameters:
5657: . dm - The DMPlex object

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


5664:   Level: intermediate

5666: .seealso: DMPlexSetAnchors(), DMGetConstraints(), DMSetConstraints()
5667: @*/
5668: PetscErrorCode DMPlexGetAnchors(DM dm, PetscSection *anchorSection, IS *anchorIS)
5669: {
5670:   DM_Plex *plex = (DM_Plex *)dm->data;

5675:   if (!plex->anchorSection && !plex->anchorIS && plex->createanchors) {(*plex->createanchors)(dm);}
5676:   if (anchorSection) *anchorSection = plex->anchorSection;
5677:   if (anchorIS) *anchorIS = plex->anchorIS;
5678:   return(0);
5679: }

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

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

5691:   collective on dm

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

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

5700:   Level: intermediate

5702: .seealso: DMPlexGetAnchors(), DMGetConstraints(), DMSetConstraints()
5703: @*/
5704: PetscErrorCode DMPlexSetAnchors(DM dm, PetscSection anchorSection, IS anchorIS)
5705: {
5706:   DM_Plex *plex = (DM_Plex *)dm->data;
5707:   PetscMPIInt result;

5712:   if (anchorSection) {
5714:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorSection),&result);
5715:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor section must have local communicator");
5716:   }
5717:   if (anchorIS) {
5719:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorIS),&result);
5720:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor IS must have local communicator");
5721:   }

5723:   PetscObjectReference((PetscObject)anchorSection);
5724:   PetscSectionDestroy(&plex->anchorSection);
5725:   plex->anchorSection = anchorSection;

5727:   PetscObjectReference((PetscObject)anchorIS);
5728:   ISDestroy(&plex->anchorIS);
5729:   plex->anchorIS = anchorIS;

5731: #if defined(PETSC_USE_DEBUG)
5732:   if (anchorIS && anchorSection) {
5733:     PetscInt size, a, pStart, pEnd;
5734:     const PetscInt *anchors;

5736:     PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5737:     ISGetLocalSize(anchorIS,&size);
5738:     ISGetIndices(anchorIS,&anchors);
5739:     for (a = 0; a < size; a++) {
5740:       PetscInt p;

5742:       p = anchors[a];
5743:       if (p >= pStart && p < pEnd) {
5744:         PetscInt dof;

5746:         PetscSectionGetDof(anchorSection,p,&dof);
5747:         if (dof) {
5748:           PetscErrorCode ierr2;

5750:           ierr2 = ISRestoreIndices(anchorIS,&anchors);CHKERRQ(ierr2);
5751:           SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point %d cannot be constrained and an anchor",p);
5752:         }
5753:       }
5754:     }
5755:     ISRestoreIndices(anchorIS,&anchors);
5756:   }
5757: #endif
5758:   /* reset the generic constraints */
5759:   DMSetDefaultConstraints(dm,NULL,NULL);
5760:   return(0);
5761: }

5765: static PetscErrorCode DMPlexCreateConstraintSection_Anchors(DM dm, PetscSection section, PetscSection *cSec)
5766: {
5767:   PetscSection anchorSection;
5768:   PetscInt pStart, pEnd, p, dof, numFields, f;

5773:   DMPlexGetAnchors(dm,&anchorSection,NULL);
5774:   PetscSectionCreate(PETSC_COMM_SELF,cSec);
5775:   PetscSectionGetNumFields(section,&numFields);
5776:   PetscSectionSetNumFields(*cSec,numFields);
5777:   PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5778:   PetscSectionSetChart(*cSec,pStart,pEnd);
5779:   for (p = pStart; p < pEnd; p++) {
5780:     PetscSectionGetDof(anchorSection,p,&dof);
5781:     if (dof) {
5782:       PetscSectionGetDof(section,p,&dof);
5783:       PetscSectionSetDof(*cSec,p,dof);
5784:       for (f = 0; f < numFields; f++) {
5785:         PetscSectionGetFieldDof(section,p,f,&dof);
5786:         PetscSectionSetFieldDof(*cSec,p,f,dof);
5787:       }
5788:     }
5789:   }
5790:   PetscSectionSetUp(*cSec);
5791:   return(0);
5792: }

5796: static PetscErrorCode DMPlexCreateConstraintMatrix_Anchors(DM dm, PetscSection section, PetscSection cSec, Mat *cMat)
5797: {
5798:   PetscSection aSec;
5799:   PetscInt pStart, pEnd, p, dof, aDof, aOff, off, nnz, annz, m, n, q, a, offset, *i, *j;
5800:   const PetscInt *anchors;
5801:   PetscInt numFields, f;
5802:   IS aIS;

5807:   PetscSectionGetStorageSize(cSec, &m);
5808:   PetscSectionGetStorageSize(section, &n);
5809:   MatCreate(PETSC_COMM_SELF,cMat);
5810:   MatSetSizes(*cMat,m,n,m,n);
5811:   MatSetType(*cMat,MATSEQAIJ);
5812:   DMPlexGetAnchors(dm,&aSec,&aIS);
5813:   ISGetIndices(aIS,&anchors);
5814:   PetscSectionGetChart(aSec,&pStart,&pEnd);
5815:   PetscMalloc1(m+1,&i);
5816:   i[0] = 0;
5817:   PetscSectionGetNumFields(section,&numFields);
5818:   for (p = pStart; p < pEnd; p++) {
5819:     PetscSectionGetDof(aSec,p,&dof);
5820:     if (!dof) continue;
5821:     PetscSectionGetOffset(aSec,p,&off);
5822:     if (numFields) {
5823:       for (f = 0; f < numFields; f++) {
5824:         annz = 0;
5825:         for (q = 0; q < dof; q++) {
5826:           a = anchors[off + q];
5827:           PetscSectionGetFieldDof(section,a,f,&aDof);
5828:           annz += aDof;
5829:         }
5830:         PetscSectionGetFieldDof(cSec,p,f,&dof);
5831:         PetscSectionGetFieldOffset(cSec,p,f,&off);
5832:         for (q = 0; q < dof; q++) {
5833:           i[off + q + 1] = i[off + q] + annz;
5834:         }
5835:       }
5836:     }
5837:     else {
5838:       annz = 0;
5839:       for (q = 0; q < dof; q++) {
5840:         a = anchors[off + q];
5841:         PetscSectionGetDof(section,a,&aDof);
5842:         annz += aDof;
5843:       }
5844:       PetscSectionGetDof(cSec,p,&dof);
5845:       PetscSectionGetOffset(cSec,p,&off);
5846:       for (q = 0; q < dof; q++) {
5847:         i[off + q + 1] = i[off + q] + annz;
5848:       }
5849:     }
5850:   }
5851:   nnz = i[m];
5852:   PetscMalloc1(nnz,&j);
5853:   offset = 0;
5854:   for (p = pStart; p < pEnd; p++) {
5855:     if (numFields) {
5856:       for (f = 0; f < numFields; f++) {
5857:         PetscSectionGetFieldDof(cSec,p,f,&dof);
5858:         for (q = 0; q < dof; q++) {
5859:           PetscInt rDof, rOff, r;
5860:           PetscSectionGetDof(aSec,p,&rDof);
5861:           PetscSectionGetOffset(aSec,p,&rOff);
5862:           for (r = 0; r < rDof; r++) {
5863:             PetscInt s;

5865:             a = anchors[rOff + r];

5867:             PetscSectionGetFieldDof(section,a,f,&aDof);
5868:             PetscSectionGetFieldOffset(section,a,f,&aOff);
5869:             for (s = 0; s < aDof; s++) {
5870:               j[offset++] = aOff + s;
5871:             }
5872:           }
5873:         }
5874:       }
5875:     }
5876:     else {
5877:       PetscSectionGetDof(cSec,p,&dof);
5878:       for (q = 0; q < dof; q++) {
5879:         PetscInt rDof, rOff, r;
5880:         PetscSectionGetDof(aSec,p,&rDof);
5881:         PetscSectionGetOffset(aSec,p,&rOff);
5882:         for (r = 0; r < rDof; r++) {
5883:           PetscInt s;

5885:           a = anchors[rOff + r];

5887:           PetscSectionGetDof(section,a,&aDof);
5888:           PetscSectionGetOffset(section,a,&aOff);
5889:           for (s = 0; s < aDof; s++) {
5890:             j[offset++] = aOff + s;
5891:           }
5892:         }
5893:       }
5894:     }
5895:   }
5896:   MatSeqAIJSetPreallocationCSR(*cMat,i,j,NULL);
5897:   PetscFree(i);
5898:   PetscFree(j);
5899:   ISRestoreIndices(aIS,&anchors);
5900:   return(0);
5901: }

5905: PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm)
5906: {
5907:   DM_Plex        *plex = (DM_Plex *)dm->data;
5908:   PetscSection   anchorSection, section, cSec;
5909:   Mat            cMat;

5914:   DMPlexGetAnchors(dm,&anchorSection,NULL);
5915:   if (anchorSection) {
5916:     PetscDS  ds;
5917:     PetscInt nf;

5919:     DMGetDefaultSection(dm,&section);
5920:     DMPlexCreateConstraintSection_Anchors(dm,section,&cSec);
5921:     DMPlexCreateConstraintMatrix_Anchors(dm,section,cSec,&cMat);
5922:     DMGetDS(dm,&ds);
5923:     PetscDSGetNumFields(ds,&nf);
5924:     if (nf && plex->computeanchormatrix) {(*plex->computeanchormatrix)(dm,section,cSec,cMat);}
5925:     DMSetDefaultConstraints(dm,cSec,cMat);
5926:     PetscSectionDestroy(&cSec);
5927:     MatDestroy(&cMat);
5928:   }
5929:   return(0);
5930: }