Actual source code: plex.c

petsc-dev 2014-08-19
Report Typos and Errors
  1: #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <../src/sys/utils/hash.h>
  3: #include <petsc-private/isimpl.h>
  4: #include <petscsf.h>

  6: /* Logging support */
  7: PetscLogEvent DMPLEX_Interpolate, DMPLEX_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeField, DMPLEX_DistributeData, 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:   DMPlexGetDimension(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) {DMPlexInsertBoundaryValuesFEM(dm, v);}
 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(DM dm, PetscViewer viewer)
193: {
194:   DM_Plex          *mesh = (DM_Plex*) dm->data;
195:   DM                cdm;
196:   DMLabel           markers;
197:   PetscSection      coordSection;
198:   Vec               coordinates;
199:   PetscViewerFormat format;
200:   PetscErrorCode    ierr;

203:   DMGetCoordinateDM(dm, &cdm);
204:   DMGetDefaultSection(cdm, &coordSection);
205:   DMGetCoordinatesLocal(dm, &coordinates);
206:   PetscViewerGetFormat(viewer, &format);
207:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
208:     const char *name;
209:     PetscInt    maxConeSize, maxSupportSize;
210:     PetscInt    pStart, pEnd, p;
211:     PetscMPIInt rank, size;

213:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
214:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
215:     PetscObjectGetName((PetscObject) dm, &name);
216:     DMPlexGetChart(dm, &pStart, &pEnd);
217:     DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
218:     PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
219:     PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
220:     PetscViewerASCIISynchronizedPrintf(viewer, "Max sizes cone: %D support: %D\n", maxConeSize, maxSupportSize);
221:     PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
222:     PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
223:     for (p = pStart; p < pEnd; ++p) {
224:       PetscInt dof, off, s;

226:       PetscSectionGetDof(mesh->supportSection, p, &dof);
227:       PetscSectionGetOffset(mesh->supportSection, p, &off);
228:       for (s = off; s < off+dof; ++s) {
229:         PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D ----> %D\n", rank, p, mesh->supports[s]);
230:       }
231:     }
232:     PetscViewerFlush(viewer);
233:     PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
234:     for (p = pStart; p < pEnd; ++p) {
235:       PetscInt dof, off, c;

237:       PetscSectionGetDof(mesh->coneSection, p, &dof);
238:       PetscSectionGetOffset(mesh->coneSection, p, &off);
239:       for (c = off; c < off+dof; ++c) {
240:         PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D <---- %D (%D)\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]);
241:       }
242:     }
243:     PetscViewerFlush(viewer);
244:     PetscSectionGetChart(coordSection, &pStart, NULL);
245:     if (pStart >= 0) {PetscSectionVecView(coordSection, coordinates, viewer);}
246:     DMPlexGetLabel(dm, "marker", &markers);
247:     DMLabelView(markers,viewer);
248:     if (size > 1) {
249:       PetscSF sf;

251:       DMGetPointSF(dm, &sf);
252:       PetscSFView(sf, viewer);
253:     }
254:     PetscViewerFlush(viewer);
255:   } else if (format == PETSC_VIEWER_ASCII_LATEX) {
256:     const char  *name;
257:     const char  *colors[3] = {"red", "blue", "green"};
258:     const int    numColors  = 3;
259:     PetscReal    scale      = 2.0;
260:     PetscScalar *coords;
261:     PetscInt     depth, cStart, cEnd, c, vStart, vEnd, v, eStart = 0, eEnd = 0, e, p;
262:     PetscMPIInt  rank, size;

264:     DMPlexGetDepth(dm, &depth);
265:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
266:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
267:     PetscObjectGetName((PetscObject) dm, &name);
268:     PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
269:     PetscViewerASCIIPrintf(viewer, "\
270: \\documentclass[crop,multi=false]{standalone}\n\n\
271: \\usepackage{tikz}\n\
272: \\usepackage{pgflibraryshapes}\n\
273: \\usetikzlibrary{backgrounds}\n\
274: \\usetikzlibrary{arrows}\n\
275: \\begin{document}\n\
276: \\section{%s}\n\
277: \\begin{center}\n", name, 8.0/scale);
278:     PetscViewerASCIIPrintf(viewer, "Mesh for process ");
279:     for (p = 0; p < size; ++p) {
280:       if (p > 0 && p == size-1) {
281:         PetscViewerASCIIPrintf(viewer, ", and ", colors[p%numColors], p);
282:       } else if (p > 0) {
283:         PetscViewerASCIIPrintf(viewer, ", ", colors[p%numColors], p);
284:       }
285:       PetscViewerASCIIPrintf(viewer, "{\\textcolor{%s}%D}", colors[p%numColors], p);
286:     }
287:     PetscViewerASCIIPrintf(viewer, ".\n\n\n\
288: \\begin{tikzpicture}[scale = %g,font=\\fontsize{8}{8}\\selectfont]\n");
289:     /* Plot vertices */
290:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
291:     VecGetArray(coordinates, &coords);
292:     PetscViewerASCIIPrintf(viewer, "\\path\n");
293:     for (v = vStart; v < vEnd; ++v) {
294:       PetscInt off, dof, d;

296:       PetscSectionGetDof(coordSection, v, &dof);
297:       PetscSectionGetOffset(coordSection, v, &off);
298:       PetscViewerASCIISynchronizedPrintf(viewer, "(");
299:       for (d = 0; d < dof; ++d) {
300:         if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
301:         PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)(scale*PetscRealPart(coords[off+d])));
302:       }
303:       PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [draw,shape=circle,color=%s] {%D} --\n", v, rank, colors[rank%numColors], v);
304:     }
305:     VecRestoreArray(coordinates, &coords);
306:     PetscViewerFlush(viewer);
307:     PetscViewerASCIIPrintf(viewer, "(0,0);\n");
308:     /* Plot edges */
309:     VecGetArray(coordinates, &coords);
310:     PetscViewerASCIIPrintf(viewer, "\\path\n");
311:     if (depth > 1) {DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);}
312:     for (e = eStart; e < eEnd; ++e) {
313:       const PetscInt *cone;
314:       PetscInt        coneSize, offA, offB, dof, d;

316:       DMPlexGetConeSize(dm, e, &coneSize);
317:       if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %d cone should have two vertices, not %d", e, coneSize);
318:       DMPlexGetCone(dm, e, &cone);
319:       PetscSectionGetDof(coordSection, cone[0], &dof);
320:       PetscSectionGetOffset(coordSection, cone[0], &offA);
321:       PetscSectionGetOffset(coordSection, cone[1], &offB);
322:       PetscViewerASCIISynchronizedPrintf(viewer, "(");
323:       for (d = 0; d < dof; ++d) {
324:         if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
325:         PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)(scale*0.5*PetscRealPart(coords[offA+d]+coords[offB+d])));
326:       }
327:       PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [draw,shape=circle,color=%s] {%D} --\n", e, rank, colors[rank%numColors], e);
328:     }
329:     VecRestoreArray(coordinates, &coords);
330:     PetscViewerFlush(viewer);
331:     PetscViewerASCIIPrintf(viewer, "(0,0);\n");
332:     /* Plot cells */
333:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
334:     for (c = cStart; c < cEnd; ++c) {
335:       PetscInt *closure = NULL;
336:       PetscInt  closureSize, firstPoint = -1;

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

343:         if ((point < vStart) || (point >= vEnd)) continue;
344:         if (firstPoint >= 0) {PetscViewerASCIISynchronizedPrintf(viewer, " -- ");}
345:         PetscViewerASCIISynchronizedPrintf(viewer, "(%D_%D)", point, rank);
346:         if (firstPoint < 0) firstPoint = point;
347:       }
348:       /* Why doesn't this work? PetscViewerASCIISynchronizedPrintf(viewer, " -- cycle;\n"); */
349:       PetscViewerASCIISynchronizedPrintf(viewer, " -- (%D_%D);\n", firstPoint, rank);
350:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
351:     }
352:     PetscViewerFlush(viewer);
353:     PetscViewerASCIIPrintf(viewer, "\\end{tikzpicture}\n\\end{center}\n");
354:     PetscViewerASCIIPrintf(viewer, "\\end{document}\n", name);
355:   } else {
356:     MPI_Comm    comm;
357:     PetscInt   *sizes, *hybsizes;
358:     PetscInt    locDepth, depth, dim, d, pMax[4];
359:     PetscInt    pStart, pEnd, p;
360:     PetscInt    numLabels, l;
361:     const char *name;
362:     PetscMPIInt size;

364:     PetscObjectGetComm((PetscObject)dm,&comm);
365:     MPI_Comm_size(comm, &size);
366:     DMPlexGetDimension(dm, &dim);
367:     PetscObjectGetName((PetscObject) dm, &name);
368:     if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dim);}
369:     else      {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dim);}
370:     DMPlexGetDepth(dm, &locDepth);
371:     MPI_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
372:     DMPlexGetHybridBounds(dm, &pMax[depth], depth > 0 ? &pMax[depth-1] : NULL, &pMax[1], &pMax[0]);
373:     PetscMalloc2(size,&sizes,size,&hybsizes);
374:     if (depth == 1) {
375:       DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd);
376:       pEnd = pEnd - pStart;
377:       MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
378:       PetscViewerASCIIPrintf(viewer, "  %D-cells:", 0);
379:       for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
380:       PetscViewerASCIIPrintf(viewer, "\n");
381:       DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd);
382:       pEnd = pEnd - pStart;
383:       MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
384:       PetscViewerASCIIPrintf(viewer, "  %D-cells:", dim);
385:       for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
386:       PetscViewerASCIIPrintf(viewer, "\n");
387:     } else {
388:       for (d = 0; d <= dim; d++) {
389:         DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
390:         pEnd    -= pStart;
391:         pMax[d] -= pStart;
392:         MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
393:         MPI_Gather(&pMax[d], 1, MPIU_INT, hybsizes, 1, MPIU_INT, 0, comm);
394:         PetscViewerASCIIPrintf(viewer, "  %D-cells:", d);
395:         for (p = 0; p < size; ++p) {
396:           if (hybsizes[p] >= 0) {PetscViewerASCIIPrintf(viewer, " %D (%D)", sizes[p], sizes[p] - hybsizes[p]);}
397:           else                  {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
398:         }
399:         PetscViewerASCIIPrintf(viewer, "\n");
400:       }
401:     }
402:     PetscFree2(sizes,hybsizes);
403:     DMPlexGetNumLabels(dm, &numLabels);
404:     if (numLabels) {PetscViewerASCIIPrintf(viewer, "Labels:\n");}
405:     for (l = 0; l < numLabels; ++l) {
406:       DMLabel         label;
407:       const char     *name;
408:       IS              valueIS;
409:       const PetscInt *values;
410:       PetscInt        numValues, v;

412:       DMPlexGetLabelName(dm, l, &name);
413:       DMPlexGetLabel(dm, name, &label);
414:       DMLabelGetNumValues(label, &numValues);
415:       PetscViewerASCIIPrintf(viewer, "  %s: %d strata of sizes (", name, numValues);
416:       DMLabelGetValueIS(label, &valueIS);
417:       ISGetIndices(valueIS, &values);
418:       for (v = 0; v < numValues; ++v) {
419:         PetscInt size;

421:         DMLabelGetStratumSize(label, values[v], &size);
422:         if (v > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
423:         PetscViewerASCIIPrintf(viewer, "%d", size);
424:       }
425:       PetscViewerASCIIPrintf(viewer, ")\n");
426:       ISRestoreIndices(valueIS, &values);
427:       ISDestroy(&valueIS);
428:     }
429:   }
430:   return(0);
431: }

435: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
436: {
437:   PetscBool      iascii, ishdf5;

443:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
444:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,  &ishdf5);
445:   if (iascii) {
446:     DMPlexView_Ascii(dm, viewer);
447:   } else if (ishdf5) {
448: #if defined(PETSC_HAVE_HDF5)
449:     PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
450:     DMPlexView_HDF5(dm, viewer);
451:     PetscViewerPopFormat(viewer);
452: #else
453:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
454: #endif
455:   }
456:   return(0);
457: }

461: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
462: {
463:   PetscBool      isbinary, ishdf5;

469:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
470:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,   &ishdf5);
471:   if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");}
472:   else if (ishdf5) {
473: #if defined(PETSC_HAVE_HDF5)
474:     DMPlexLoad_HDF5(dm, viewer);
475: #else
476:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
477: #endif
478:   }
479:   return(0);
480: }

484: static PetscErrorCode BoundaryDestroy(DMBoundary *boundary)
485: {
486:   DMBoundary     b, next;

490:   if (!boundary) return(0);
491:   b = *boundary;
492:   *boundary = NULL;
493:   for (; b; b = next) {
494:     next = b->next;
495:     PetscFree(b->ids);
496:     PetscFree(b->name);
497:     PetscFree(b->labelname);
498:     PetscFree(b);
499:   }
500:   return(0);
501: }

505: PetscErrorCode DMDestroy_Plex(DM dm)
506: {
507:   DM_Plex       *mesh = (DM_Plex*) dm->data;
508:   DMLabel        next  = mesh->labels;

512:   if (--mesh->refct > 0) return(0);
513:   PetscSectionDestroy(&mesh->coneSection);
514:   PetscFree(mesh->cones);
515:   PetscFree(mesh->coneOrientations);
516:   PetscSectionDestroy(&mesh->supportSection);
517:   PetscFree(mesh->supports);
518:   PetscFree(mesh->facesTmp);
519:   while (next) {
520:     DMLabel tmp = next->next;

522:     DMLabelDestroy(&next);
523:     next = tmp;
524:   }
525:   DMDestroy(&mesh->coarseMesh);
526:   DMLabelDestroy(&mesh->subpointMap);
527:   ISDestroy(&mesh->globalVertexNumbers);
528:   ISDestroy(&mesh->globalCellNumbers);
529:   BoundaryDestroy(&mesh->boundary);
530:   /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
531:   PetscFree(mesh);
532:   return(0);
533: }

537: PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
538: {
539:   PetscSection   section, sectionGlobal;
540:   PetscInt       bs = -1;
541:   PetscInt       localSize;
542:   PetscBool      isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock;
544:   MatType        mtype;

547:   MatInitializePackage();
548:   mtype = dm->mattype;
549:   DMGetDefaultSection(dm, &section);
550:   DMGetDefaultGlobalSection(dm, &sectionGlobal);
551:   /* PetscSectionGetStorageSize(sectionGlobal, &localSize); */
552:   PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
553:   MatCreate(PetscObjectComm((PetscObject)dm), J);
554:   MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
555:   MatSetType(*J, mtype);
556:   MatSetFromOptions(*J);
557:   PetscStrcmp(mtype, MATSHELL, &isShell);
558:   PetscStrcmp(mtype, MATBAIJ, &isBlock);
559:   PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
560:   PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
561:   PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
562:   PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
563:   PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
564:   if (!isShell) {
565:     PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
566:     PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin;

568:     if (bs < 0) {
569:       if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
570:         PetscInt pStart, pEnd, p, dof, cdof;

572:         PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
573:         for (p = pStart; p < pEnd; ++p) {
574:           PetscSectionGetDof(sectionGlobal, p, &dof);
575:           PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
576:           if (dof-cdof) {
577:             if (bs < 0) {
578:               bs = dof-cdof;
579:             } else if (bs != dof-cdof) {
580:               /* Layout does not admit a pointwise block size */
581:               bs = 1;
582:               break;
583:             }
584:           }
585:         }
586:         /* Must have same blocksize on all procs (some might have no points) */
587:         bsLocal = bs;
588:         MPI_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
589:         bsLocal = bs < 0 ? bsMax : bs;
590:         MPI_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
591:         if (bsMin != bsMax) {
592:           bs = 1;
593:         } else {
594:           bs = bsMax;
595:         }
596:       } else {
597:         bs = 1;
598:       }
599:     }
600:     PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
601:     DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix);
602:     PetscFree4(dnz, onz, dnzu, onzu);
603:   }
604:   return(0);
605: }

609: /*@
610:   DMPlexGetDimension - Return the topological mesh dimension

612:   Not collective

614:   Input Parameter:
615: . mesh - The DMPlex

617:   Output Parameter:
618: . dim - The topological mesh dimension

620:   Level: beginner

622: .seealso: DMPlexCreate()
623: @*/
624: PetscErrorCode DMPlexGetDimension(DM dm, PetscInt *dim)
625: {
626:   DM_Plex *mesh = (DM_Plex*) dm->data;

631:   *dim = mesh->dim;
632:   return(0);
633: }

637: /*@
638:   DMPlexSetDimension - Set the topological mesh dimension

640:   Collective on mesh

642:   Input Parameters:
643: + mesh - The DMPlex
644: - dim - The topological mesh dimension

646:   Level: beginner

648: .seealso: DMPlexCreate()
649: @*/
650: PetscErrorCode DMPlexSetDimension(DM dm, PetscInt dim)
651: {
652:   DM_Plex *mesh = (DM_Plex*) dm->data;

657:   mesh->dim = dim;
658:   return(0);
659: }

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

666:   Not collective

668:   Input Parameter:
669: . mesh - The DMPlex

671:   Output Parameters:
672: + pStart - The first mesh point
673: - pEnd   - The upper bound for mesh points

675:   Level: beginner

677: .seealso: DMPlexCreate(), DMPlexSetChart()
678: @*/
679: PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
680: {
681:   DM_Plex       *mesh = (DM_Plex*) dm->data;

686:   PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
687:   return(0);
688: }

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

695:   Not collective

697:   Input Parameters:
698: + mesh - The DMPlex
699: . pStart - The first mesh point
700: - pEnd   - The upper bound for mesh points

702:   Output Parameters:

704:   Level: beginner

706: .seealso: DMPlexCreate(), DMPlexGetChart()
707: @*/
708: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
709: {
710:   DM_Plex       *mesh = (DM_Plex*) dm->data;

715:   PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
716:   PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
717:   return(0);
718: }

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

725:   Not collective

727:   Input Parameters:
728: + mesh - The DMPlex
729: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

731:   Output Parameter:
732: . size - The cone size for point p

734:   Level: beginner

736: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
737: @*/
738: PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size)
739: {
740:   DM_Plex       *mesh = (DM_Plex*) dm->data;

746:   PetscSectionGetDof(mesh->coneSection, p, size);
747:   return(0);
748: }

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

755:   Not collective

757:   Input Parameters:
758: + mesh - The DMPlex
759: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
760: - size - The cone size for point p

762:   Output Parameter:

764:   Note:
765:   This should be called after DMPlexSetChart().

767:   Level: beginner

769: .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart()
770: @*/
771: PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size)
772: {
773:   DM_Plex       *mesh = (DM_Plex*) dm->data;

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

780:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
781:   return(0);
782: }

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

789:   Not collective

791:   Input Parameters:
792: + mesh - The DMPlex
793: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
794: - size - The additional cone size for point p

796:   Output Parameter:

798:   Note:
799:   This should be called after DMPlexSetChart().

801:   Level: beginner

803: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart()
804: @*/
805: PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size)
806: {
807:   DM_Plex       *mesh = (DM_Plex*) dm->data;
808:   PetscInt       csize;

813:   PetscSectionAddDof(mesh->coneSection, p, size);
814:   PetscSectionGetDof(mesh->coneSection, p, &csize);

816:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
817:   return(0);
818: }

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

825:   Not collective

827:   Input Parameters:
828: + mesh - The DMPlex
829: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

834:   Level: beginner

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

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

842: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart()
843: @*/
844: PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
845: {
846:   DM_Plex       *mesh = (DM_Plex*) dm->data;
847:   PetscInt       off;

853:   PetscSectionGetOffset(mesh->coneSection, p, &off);
854:   *cone = &mesh->cones[off];
855:   return(0);
856: }

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

863:   Not collective

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

870:   Output Parameter:

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

875:   Level: beginner

877: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
878: @*/
879: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
880: {
881:   DM_Plex       *mesh = (DM_Plex*) dm->data;
882:   PetscInt       pStart, pEnd;
883:   PetscInt       dof, off, c;

888:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
889:   PetscSectionGetDof(mesh->coneSection, p, &dof);
891:   PetscSectionGetOffset(mesh->coneSection, p, &off);
892:   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);
893:   for (c = 0; c < dof; ++c) {
894:     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);
895:     mesh->cones[off+c] = cone[c];
896:   }
897:   return(0);
898: }

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

905:   Not collective

907:   Input Parameters:
908: + mesh - The DMPlex
909: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

917:   Level: beginner

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

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

925: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
926: @*/
927: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
928: {
929:   DM_Plex       *mesh = (DM_Plex*) dm->data;
930:   PetscInt       off;

935: #if defined(PETSC_USE_DEBUG)
936:   {
937:     PetscInt dof;
938:     PetscSectionGetDof(mesh->coneSection, p, &dof);
940:   }
941: #endif
942:   PetscSectionGetOffset(mesh->coneSection, p, &off);

944:   *coneOrientation = &mesh->coneOrientations[off];
945:   return(0);
946: }

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

953:   Not collective

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

963:   Output Parameter:

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

968:   Level: beginner

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

981:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
982:   PetscSectionGetDof(mesh->coneSection, p, &dof);
984:   PetscSectionGetOffset(mesh->coneSection, p, &off);
985:   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);
986:   for (c = 0; c < dof; ++c) {
987:     PetscInt cdof, o = coneOrientation[c];

989:     PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
990:     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);
991:     mesh->coneOrientations[off+c] = o;
992:   }
993:   return(0);
994: }

998: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
999: {
1000:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1001:   PetscInt       pStart, pEnd;
1002:   PetscInt       dof, off;

1007:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1008:   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);
1009:   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);
1010:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1011:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1012:   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);
1013:   mesh->cones[off+conePos] = conePoint;
1014:   return(0);
1015: }

1019: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1020: {
1021:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1022:   PetscInt       pStart, pEnd;
1023:   PetscInt       dof, off;

1028:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1029:   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);
1030:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1031:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1032:   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);
1033:   mesh->coneOrientations[off+conePos] = coneOrientation;
1034:   return(0);
1035: }

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

1042:   Not collective

1044:   Input Parameters:
1045: + mesh - The DMPlex
1046: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1048:   Output Parameter:
1049: . size - The support size for point p

1051:   Level: beginner

1053: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1054: @*/
1055: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1056: {
1057:   DM_Plex       *mesh = (DM_Plex*) dm->data;

1063:   PetscSectionGetDof(mesh->supportSection, p, size);
1064:   return(0);
1065: }

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

1072:   Not collective

1074:   Input Parameters:
1075: + mesh - The DMPlex
1076: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1077: - size - The support size for point p

1079:   Output Parameter:

1081:   Note:
1082:   This should be called after DMPlexSetChart().

1084:   Level: beginner

1086: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1087: @*/
1088: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1089: {
1090:   DM_Plex       *mesh = (DM_Plex*) dm->data;

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

1097:   mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1098:   return(0);
1099: }

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

1106:   Not collective

1108:   Input Parameters:
1109: + mesh - The DMPlex
1110: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1115:   Level: beginner

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

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

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

1134:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1135:   *support = &mesh->supports[off];
1136:   return(0);
1137: }

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

1144:   Not collective

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

1151:   Output Parameter:

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

1156:   Level: beginner

1158: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1159: @*/
1160: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1161: {
1162:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1163:   PetscInt       pStart, pEnd;
1164:   PetscInt       dof, off, c;

1169:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1170:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1172:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1173:   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);
1174:   for (c = 0; c < dof; ++c) {
1175:     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);
1176:     mesh->supports[off+c] = support[c];
1177:   }
1178:   return(0);
1179: }

1183: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1184: {
1185:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1186:   PetscInt       pStart, pEnd;
1187:   PetscInt       dof, off;

1192:   PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1193:   PetscSectionGetDof(mesh->supportSection, p, &dof);
1194:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1195:   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);
1196:   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);
1197:   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);
1198:   mesh->supports[off+supportPos] = supportPoint;
1199:   return(0);
1200: }

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

1207:   Not collective

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

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

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

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

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

1228:   Level: beginner

1230: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1231: @*/
1232: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1233: {
1234:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1235:   PetscInt       *closure, *fifo;
1236:   const PetscInt *tmp = NULL, *tmpO = NULL;
1237:   PetscInt        tmpSize, t;
1238:   PetscInt        depth       = 0, maxSize;
1239:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1240:   PetscErrorCode  ierr;

1244:   DMPlexGetDepth(dm, &depth);
1245:   /* This is only 1-level */
1246:   if (useCone) {
1247:     DMPlexGetConeSize(dm, p, &tmpSize);
1248:     DMPlexGetCone(dm, p, &tmp);
1249:     DMPlexGetConeOrientation(dm, p, &tmpO);
1250:   } else {
1251:     DMPlexGetSupportSize(dm, p, &tmpSize);
1252:     DMPlexGetSupport(dm, p, &tmp);
1253:   }
1254:   if (depth == 1) {
1255:     if (*points) {
1256:       closure = *points;
1257:     } else {
1258:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1259:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1260:     }
1261:     closure[0] = p; closure[1] = 0;
1262:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1263:       closure[closureSize]   = tmp[t];
1264:       closure[closureSize+1] = tmpO ? tmpO[t] : 0;
1265:     }
1266:     if (numPoints) *numPoints = closureSize/2;
1267:     if (points)    *points    = closure;
1268:     return(0);
1269:   }
1270:   maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth+1);
1271:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1272:   if (*points) {
1273:     closure = *points;
1274:   } else {
1275:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1276:   }
1277:   closure[0] = p; closure[1] = 0;
1278:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1279:     const PetscInt cp = tmp[t];
1280:     const PetscInt co = tmpO ? tmpO[t] : 0;

1282:     closure[closureSize]   = cp;
1283:     closure[closureSize+1] = co;
1284:     fifo[fifoSize]         = cp;
1285:     fifo[fifoSize+1]       = co;
1286:   }
1287:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1288:   while (fifoSize - fifoStart) {
1289:     const PetscInt q   = fifo[fifoStart];
1290:     const PetscInt o   = fifo[fifoStart+1];
1291:     const PetscInt rev = o >= 0 ? 0 : 1;
1292:     const PetscInt off = rev ? -(o+1) : o;

1294:     if (useCone) {
1295:       DMPlexGetConeSize(dm, q, &tmpSize);
1296:       DMPlexGetCone(dm, q, &tmp);
1297:       DMPlexGetConeOrientation(dm, q, &tmpO);
1298:     } else {
1299:       DMPlexGetSupportSize(dm, q, &tmpSize);
1300:       DMPlexGetSupport(dm, q, &tmp);
1301:       tmpO = NULL;
1302:     }
1303:     for (t = 0; t < tmpSize; ++t) {
1304:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1305:       const PetscInt cp = tmp[i];
1306:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1307:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1308:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1309:       PetscInt       co = tmpO ? tmpO[i] : 0;
1310:       PetscInt       c;

1312:       if (rev) {
1313:         PetscInt childSize, coff;
1314:         DMPlexGetConeSize(dm, cp, &childSize);
1315:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1316:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1317:       }
1318:       /* Check for duplicate */
1319:       for (c = 0; c < closureSize; c += 2) {
1320:         if (closure[c] == cp) break;
1321:       }
1322:       if (c == closureSize) {
1323:         closure[closureSize]   = cp;
1324:         closure[closureSize+1] = co;
1325:         fifo[fifoSize]         = cp;
1326:         fifo[fifoSize+1]       = co;
1327:         closureSize           += 2;
1328:         fifoSize              += 2;
1329:       }
1330:     }
1331:     fifoStart += 2;
1332:   }
1333:   if (numPoints) *numPoints = closureSize/2;
1334:   if (points)    *points    = closure;
1335:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1336:   return(0);
1337: }

1341: /*@C
1342:   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

1344:   Not collective

1346:   Input Parameters:
1347: + mesh - The DMPlex
1348: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1349: . orientation - The orientation of the point
1350: . useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
1351: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used

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

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

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

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

1366:   Level: beginner

1368: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1369: @*/
1370: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1371: {
1372:   DM_Plex        *mesh = (DM_Plex*) dm->data;
1373:   PetscInt       *closure, *fifo;
1374:   const PetscInt *tmp = NULL, *tmpO = NULL;
1375:   PetscInt        tmpSize, t;
1376:   PetscInt        depth       = 0, maxSize;
1377:   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
1378:   PetscErrorCode  ierr;

1382:   DMPlexGetDepth(dm, &depth);
1383:   /* This is only 1-level */
1384:   if (useCone) {
1385:     DMPlexGetConeSize(dm, p, &tmpSize);
1386:     DMPlexGetCone(dm, p, &tmp);
1387:     DMPlexGetConeOrientation(dm, p, &tmpO);
1388:   } else {
1389:     DMPlexGetSupportSize(dm, p, &tmpSize);
1390:     DMPlexGetSupport(dm, p, &tmp);
1391:   }
1392:   if (depth == 1) {
1393:     if (*points) {
1394:       closure = *points;
1395:     } else {
1396:       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1397:       DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1398:     }
1399:     closure[0] = p; closure[1] = ornt;
1400:     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1401:       const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1402:       closure[closureSize]   = tmp[i];
1403:       closure[closureSize+1] = tmpO ? tmpO[i] : 0;
1404:     }
1405:     if (numPoints) *numPoints = closureSize/2;
1406:     if (points)    *points    = closure;
1407:     return(0);
1408:   }
1409:   maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth+1);
1410:   DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1411:   if (*points) {
1412:     closure = *points;
1413:   } else {
1414:     DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1415:   }
1416:   closure[0] = p; closure[1] = ornt;
1417:   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1418:     const PetscInt i  = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1419:     const PetscInt cp = tmp[i];
1420:     PetscInt       co = tmpO ? tmpO[i] : 0;

1422:     if (ornt < 0) {
1423:       PetscInt childSize, coff;
1424:       DMPlexGetConeSize(dm, cp, &childSize);
1425:       coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1426:       co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1427:     }
1428:     closure[closureSize]   = cp;
1429:     closure[closureSize+1] = co;
1430:     fifo[fifoSize]         = cp;
1431:     fifo[fifoSize+1]       = co;
1432:   }
1433:   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1434:   while (fifoSize - fifoStart) {
1435:     const PetscInt q   = fifo[fifoStart];
1436:     const PetscInt o   = fifo[fifoStart+1];
1437:     const PetscInt rev = o >= 0 ? 0 : 1;
1438:     const PetscInt off = rev ? -(o+1) : o;

1440:     if (useCone) {
1441:       DMPlexGetConeSize(dm, q, &tmpSize);
1442:       DMPlexGetCone(dm, q, &tmp);
1443:       DMPlexGetConeOrientation(dm, q, &tmpO);
1444:     } else {
1445:       DMPlexGetSupportSize(dm, q, &tmpSize);
1446:       DMPlexGetSupport(dm, q, &tmp);
1447:       tmpO = NULL;
1448:     }
1449:     for (t = 0; t < tmpSize; ++t) {
1450:       const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
1451:       const PetscInt cp = tmp[i];
1452:       /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1453:       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1454:        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1455:       PetscInt       co = tmpO ? tmpO[i] : 0;
1456:       PetscInt       c;

1458:       if (rev) {
1459:         PetscInt childSize, coff;
1460:         DMPlexGetConeSize(dm, cp, &childSize);
1461:         coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1462:         co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1463:       }
1464:       /* Check for duplicate */
1465:       for (c = 0; c < closureSize; c += 2) {
1466:         if (closure[c] == cp) break;
1467:       }
1468:       if (c == closureSize) {
1469:         closure[closureSize]   = cp;
1470:         closure[closureSize+1] = co;
1471:         fifo[fifoSize]         = cp;
1472:         fifo[fifoSize+1]       = co;
1473:         closureSize           += 2;
1474:         fifoSize              += 2;
1475:       }
1476:     }
1477:     fifoStart += 2;
1478:   }
1479:   if (numPoints) *numPoints = closureSize/2;
1480:   if (points)    *points    = closure;
1481:   DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1482:   return(0);
1483: }

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

1490:   Not collective

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

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

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

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

1508:   Level: beginner

1510: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1511: @*/
1512: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1513: {

1520:   DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1521:   if (numPoints) *numPoints = 0;
1522:   return(0);
1523: }

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

1530:   Not collective

1532:   Input Parameter:
1533: . mesh - The DMPlex

1535:   Output Parameters:
1536: + maxConeSize - The maximum number of in-edges
1537: - maxSupportSize - The maximum number of out-edges

1539:   Level: beginner

1541: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1542: @*/
1543: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1544: {
1545:   DM_Plex *mesh = (DM_Plex*) dm->data;

1549:   if (maxConeSize)    *maxConeSize    = mesh->maxConeSize;
1550:   if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1551:   return(0);
1552: }

1556: PetscErrorCode DMSetUp_Plex(DM dm)
1557: {
1558:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1559:   PetscInt       size;

1564:   PetscSectionSetUp(mesh->coneSection);
1565:   PetscSectionGetStorageSize(mesh->coneSection, &size);
1566:   PetscMalloc1(size, &mesh->cones);
1567:   PetscCalloc1(size, &mesh->coneOrientations);
1568:   if (mesh->maxSupportSize) {
1569:     PetscSectionSetUp(mesh->supportSection);
1570:     PetscSectionGetStorageSize(mesh->supportSection, &size);
1571:     PetscMalloc1(size, &mesh->supports);
1572:   }
1573:   return(0);
1574: }

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

1583:   if (subdm) {DMClone(dm, subdm);}
1584:   DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1585:   return(0);
1586: }

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

1593:   Not collective

1595:   Input Parameter:
1596: . mesh - The DMPlex

1598:   Output Parameter:

1600:   Note:
1601:   This should be called after all calls to DMPlexSetCone()

1603:   Level: beginner

1605: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1606: @*/
1607: PetscErrorCode DMPlexSymmetrize(DM dm)
1608: {
1609:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1610:   PetscInt      *offsets;
1611:   PetscInt       supportSize;
1612:   PetscInt       pStart, pEnd, p;

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

1623:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1624:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1625:     for (c = off; c < off+dof; ++c) {
1626:       PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1627:     }
1628:   }
1629:   for (p = pStart; p < pEnd; ++p) {
1630:     PetscInt dof;

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

1634:     mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1635:   }
1636:   PetscSectionSetUp(mesh->supportSection);
1637:   /* Calculate supports */
1638:   PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1639:   PetscMalloc1(supportSize, &mesh->supports);
1640:   PetscCalloc1(pEnd - pStart, &offsets);
1641:   for (p = pStart; p < pEnd; ++p) {
1642:     PetscInt dof, off, c;

1644:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1645:     PetscSectionGetOffset(mesh->coneSection, p, &off);
1646:     for (c = off; c < off+dof; ++c) {
1647:       const PetscInt q = mesh->cones[c];
1648:       PetscInt       offS;

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

1652:       mesh->supports[offS+offsets[q]] = p;
1653:       ++offsets[q];
1654:     }
1655:   }
1656:   PetscFree(offsets);
1657:   return(0);
1658: }

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

1668:   Not collective

1670:   Input Parameter:
1671: . mesh - The DMPlex

1673:   Output Parameter:

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

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

1683:   Level: beginner

1685: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1686: @*/
1687: PetscErrorCode DMPlexStratify(DM dm)
1688: {
1689:   DMLabel        label;
1690:   PetscInt       pStart, pEnd, p;
1691:   PetscInt       numRoots = 0, numLeaves = 0;

1696:   PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1697:   /* Calculate depth */
1698:   DMPlexGetChart(dm, &pStart, &pEnd);
1699:   DMPlexCreateLabel(dm, "depth");
1700:   DMPlexGetDepthLabel(dm, &label);
1701:   /* Initialize roots and count leaves */
1702:   for (p = pStart; p < pEnd; ++p) {
1703:     PetscInt coneSize, supportSize;

1705:     DMPlexGetConeSize(dm, p, &coneSize);
1706:     DMPlexGetSupportSize(dm, p, &supportSize);
1707:     if (!coneSize && supportSize) {
1708:       ++numRoots;
1709:       DMLabelSetValue(label, p, 0);
1710:     } else if (!supportSize && coneSize) {
1711:       ++numLeaves;
1712:     } else if (!supportSize && !coneSize) {
1713:       /* Isolated points */
1714:       DMLabelSetValue(label, p, 0);
1715:     }
1716:   }
1717:   if (numRoots + numLeaves == (pEnd - pStart)) {
1718:     for (p = pStart; p < pEnd; ++p) {
1719:       PetscInt coneSize, supportSize;

1721:       DMPlexGetConeSize(dm, p, &coneSize);
1722:       DMPlexGetSupportSize(dm, p, &supportSize);
1723:       if (!supportSize && coneSize) {
1724:         DMLabelSetValue(label, p, 1);
1725:       }
1726:     }
1727:   } else {
1728:     IS       pointIS;
1729:     PetscInt numPoints = 0, level = 0;

1731:     DMLabelGetStratumIS(label, level, &pointIS);
1732:     if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1733:     while (numPoints) {
1734:       const PetscInt *points;
1735:       const PetscInt  newLevel = level+1;

1737:       ISGetIndices(pointIS, &points);
1738:       for (p = 0; p < numPoints; ++p) {
1739:         const PetscInt  point = points[p];
1740:         const PetscInt *support;
1741:         PetscInt        supportSize, s;

1743:         DMPlexGetSupportSize(dm, point, &supportSize);
1744:         DMPlexGetSupport(dm, point, &support);
1745:         for (s = 0; s < supportSize; ++s) {
1746:           DMLabelSetValue(label, support[s], newLevel);
1747:         }
1748:       }
1749:       ISRestoreIndices(pointIS, &points);
1750:       ++level;
1751:       ISDestroy(&pointIS);
1752:       DMLabelGetStratumIS(label, level, &pointIS);
1753:       if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1754:       else         {numPoints = 0;}
1755:     }
1756:     ISDestroy(&pointIS);
1757:   }
1758:   PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1759:   return(0);
1760: }

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

1767:   Not Collective

1769:   Input Parameters:
1770: + dm - The DMPlex object
1771: . numPoints - The number of input points for the join
1772: - points - The input points

1774:   Output Parameters:
1775: + numCoveredPoints - The number of points in the join
1776: - coveredPoints - The points in the join

1778:   Level: intermediate

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

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

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

1788: .keywords: mesh
1789: .seealso: DMPlexRestoreJoin(), DMPlexGetMeet()
1790: @*/
1791: PetscErrorCode DMPlexGetJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1792: {
1793:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1794:   PetscInt      *join[2];
1795:   PetscInt       joinSize, i = 0;
1796:   PetscInt       dof, off, p, c, m;

1804:   DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[0]);
1805:   DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1]);
1806:   /* Copy in support of first point */
1807:   PetscSectionGetDof(mesh->supportSection, points[0], &dof);
1808:   PetscSectionGetOffset(mesh->supportSection, points[0], &off);
1809:   for (joinSize = 0; joinSize < dof; ++joinSize) {
1810:     join[i][joinSize] = mesh->supports[off+joinSize];
1811:   }
1812:   /* Check each successive support */
1813:   for (p = 1; p < numPoints; ++p) {
1814:     PetscInt newJoinSize = 0;

1816:     PetscSectionGetDof(mesh->supportSection, points[p], &dof);
1817:     PetscSectionGetOffset(mesh->supportSection, points[p], &off);
1818:     for (c = 0; c < dof; ++c) {
1819:       const PetscInt point = mesh->supports[off+c];

1821:       for (m = 0; m < joinSize; ++m) {
1822:         if (point == join[i][m]) {
1823:           join[1-i][newJoinSize++] = point;
1824:           break;
1825:         }
1826:       }
1827:     }
1828:     joinSize = newJoinSize;
1829:     i        = 1-i;
1830:   }
1831:   *numCoveredPoints = joinSize;
1832:   *coveredPoints    = join[i];
1833:   DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
1834:   return(0);
1835: }

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

1842:   Not Collective

1844:   Input Parameters:
1845: + dm - The DMPlex object
1846: . numPoints - The number of input points for the join
1847: - points - The input points

1849:   Output Parameters:
1850: + numCoveredPoints - The number of points in the join
1851: - coveredPoints - The points in the join

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

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

1859:   Level: intermediate

1861: .keywords: mesh
1862: .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet()
1863: @*/
1864: PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1865: {

1873:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
1874:   if (numCoveredPoints) *numCoveredPoints = 0;
1875:   return(0);
1876: }

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

1883:   Not Collective

1885:   Input Parameters:
1886: + dm - The DMPlex object
1887: . numPoints - The number of input points for the join
1888: - points - The input points

1890:   Output Parameters:
1891: + numCoveredPoints - The number of points in the join
1892: - coveredPoints - The points in the join

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

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

1900:   Level: intermediate

1902: .keywords: mesh
1903: .seealso: DMPlexGetJoin(), DMPlexRestoreJoin(), DMPlexGetMeet()
1904: @*/
1905: PetscErrorCode DMPlexGetFullJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1906: {
1907:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1908:   PetscInt      *offsets, **closures;
1909:   PetscInt      *join[2];
1910:   PetscInt       depth = 0, maxSize, joinSize = 0, i = 0;
1911:   PetscInt       p, d, c, m;


1920:   DMPlexGetDepth(dm, &depth);
1921:   PetscCalloc1(numPoints, &closures);
1922:   DMGetWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
1923:   maxSize = PetscPowInt(mesh->maxSupportSize,depth+1);
1924:   DMGetWorkArray(dm, maxSize, PETSC_INT, &join[0]);
1925:   DMGetWorkArray(dm, maxSize, PETSC_INT, &join[1]);

1927:   for (p = 0; p < numPoints; ++p) {
1928:     PetscInt closureSize;

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

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

1936:       DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
1937:       for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) {
1938:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
1939:           offsets[p*(depth+2)+d+1] = i;
1940:           break;
1941:         }
1942:       }
1943:       if (i == closureSize) offsets[p*(depth+2)+d+1] = i;
1944:     }
1945:     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);
1946:   }
1947:   for (d = 0; d < depth+1; ++d) {
1948:     PetscInt dof;

1950:     /* Copy in support of first point */
1951:     dof = offsets[d+1] - offsets[d];
1952:     for (joinSize = 0; joinSize < dof; ++joinSize) {
1953:       join[i][joinSize] = closures[0][(offsets[d]+joinSize)*2];
1954:     }
1955:     /* Check each successive cone */
1956:     for (p = 1; p < numPoints && joinSize; ++p) {
1957:       PetscInt newJoinSize = 0;

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

1963:         for (m = 0; m < joinSize; ++m) {
1964:           if (point == join[i][m]) {
1965:             join[1-i][newJoinSize++] = point;
1966:             break;
1967:           }
1968:         }
1969:       }
1970:       joinSize = newJoinSize;
1971:       i        = 1-i;
1972:     }
1973:     if (joinSize) break;
1974:   }
1975:   *numCoveredPoints = joinSize;
1976:   *coveredPoints    = join[i];
1977:   for (p = 0; p < numPoints; ++p) {
1978:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, NULL, &closures[p]);
1979:   }
1980:   PetscFree(closures);
1981:   DMRestoreWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
1982:   DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
1983:   return(0);
1984: }

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

1991:   Not Collective

1993:   Input Parameters:
1994: + dm - The DMPlex object
1995: . numPoints - The number of input points for the meet
1996: - points - The input points

1998:   Output Parameters:
1999: + numCoveredPoints - The number of points in the meet
2000: - coveredPoints - The points in the meet

2002:   Level: intermediate

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

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

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

2012: .keywords: mesh
2013: .seealso: DMPlexRestoreMeet(), DMPlexGetJoin()
2014: @*/
2015: PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2016: {
2017:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2018:   PetscInt      *meet[2];
2019:   PetscInt       meetSize, i = 0;
2020:   PetscInt       dof, off, p, c, m;

2028:   DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[0]);
2029:   DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1]);
2030:   /* Copy in cone of first point */
2031:   PetscSectionGetDof(mesh->coneSection, points[0], &dof);
2032:   PetscSectionGetOffset(mesh->coneSection, points[0], &off);
2033:   for (meetSize = 0; meetSize < dof; ++meetSize) {
2034:     meet[i][meetSize] = mesh->cones[off+meetSize];
2035:   }
2036:   /* Check each successive cone */
2037:   for (p = 1; p < numPoints; ++p) {
2038:     PetscInt newMeetSize = 0;

2040:     PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2041:     PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2042:     for (c = 0; c < dof; ++c) {
2043:       const PetscInt point = mesh->cones[off+c];

2045:       for (m = 0; m < meetSize; ++m) {
2046:         if (point == meet[i][m]) {
2047:           meet[1-i][newMeetSize++] = point;
2048:           break;
2049:         }
2050:       }
2051:     }
2052:     meetSize = newMeetSize;
2053:     i        = 1-i;
2054:   }
2055:   *numCoveringPoints = meetSize;
2056:   *coveringPoints    = meet[i];
2057:   DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2058:   return(0);
2059: }

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

2066:   Not Collective

2068:   Input Parameters:
2069: + dm - The DMPlex object
2070: . numPoints - The number of input points for the meet
2071: - points - The input points

2073:   Output Parameters:
2074: + numCoveredPoints - The number of points in the meet
2075: - coveredPoints - The points in the meet

2077:   Level: intermediate

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

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

2085: .keywords: mesh
2086: .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin()
2087: @*/
2088: PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2089: {

2097:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2098:   if (numCoveredPoints) *numCoveredPoints = 0;
2099:   return(0);
2100: }

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

2107:   Not Collective

2109:   Input Parameters:
2110: + dm - The DMPlex object
2111: . numPoints - The number of input points for the meet
2112: - points - The input points

2114:   Output Parameters:
2115: + numCoveredPoints - The number of points in the meet
2116: - coveredPoints - The points in the meet

2118:   Level: intermediate

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

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

2126: .keywords: mesh
2127: .seealso: DMPlexGetMeet(), DMPlexRestoreMeet(), DMPlexGetJoin()
2128: @*/
2129: PetscErrorCode DMPlexGetFullMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2130: {
2131:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2132:   PetscInt      *offsets, **closures;
2133:   PetscInt      *meet[2];
2134:   PetscInt       height = 0, maxSize, meetSize = 0, i = 0;
2135:   PetscInt       p, h, c, m;


2144:   DMPlexGetDepth(dm, &height);
2145:   PetscMalloc1(numPoints, &closures);
2146:   DMGetWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2147:   maxSize = PetscPowInt(mesh->maxConeSize,height+1);
2148:   DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[0]);
2149:   DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[1]);

2151:   for (p = 0; p < numPoints; ++p) {
2152:     PetscInt closureSize;

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

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

2160:       DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
2161:       for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) {
2162:         if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2163:           offsets[p*(height+2)+h+1] = i;
2164:           break;
2165:         }
2166:       }
2167:       if (i == closureSize) offsets[p*(height+2)+h+1] = i;
2168:     }
2169:     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);
2170:   }
2171:   for (h = 0; h < height+1; ++h) {
2172:     PetscInt dof;

2174:     /* Copy in cone of first point */
2175:     dof = offsets[h+1] - offsets[h];
2176:     for (meetSize = 0; meetSize < dof; ++meetSize) {
2177:       meet[i][meetSize] = closures[0][(offsets[h]+meetSize)*2];
2178:     }
2179:     /* Check each successive cone */
2180:     for (p = 1; p < numPoints && meetSize; ++p) {
2181:       PetscInt newMeetSize = 0;

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

2187:         for (m = 0; m < meetSize; ++m) {
2188:           if (point == meet[i][m]) {
2189:             meet[1-i][newMeetSize++] = point;
2190:             break;
2191:           }
2192:         }
2193:       }
2194:       meetSize = newMeetSize;
2195:       i        = 1-i;
2196:     }
2197:     if (meetSize) break;
2198:   }
2199:   *numCoveredPoints = meetSize;
2200:   *coveredPoints    = meet[i];
2201:   for (p = 0; p < numPoints; ++p) {
2202:     DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, NULL, &closures[p]);
2203:   }
2204:   PetscFree(closures);
2205:   DMRestoreWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2206:   DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2207:   return(0);
2208: }

2212: /*@C
2213:   DMPlexEqual - Determine if two DMs have the same topology

2215:   Not Collective

2217:   Input Parameters:
2218: + dmA - A DMPlex object
2219: - dmB - A DMPlex object

2221:   Output Parameters:
2222: . equal - PETSC_TRUE if the topologies are identical

2224:   Level: intermediate

2226:   Notes:
2227:   We are not solving graph isomorphism, so we do not permutation.

2229: .keywords: mesh
2230: .seealso: DMPlexGetCone()
2231: @*/
2232: PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal)
2233: {
2234:   PetscInt       depth, depthB, pStart, pEnd, pStartB, pEndB, p;


2242:   *equal = PETSC_FALSE;
2243:   DMPlexGetDepth(dmA, &depth);
2244:   DMPlexGetDepth(dmB, &depthB);
2245:   if (depth != depthB) return(0);
2246:   DMPlexGetChart(dmA, &pStart,  &pEnd);
2247:   DMPlexGetChart(dmB, &pStartB, &pEndB);
2248:   if ((pStart != pStartB) || (pEnd != pEndB)) return(0);
2249:   for (p = pStart; p < pEnd; ++p) {
2250:     const PetscInt *cone, *coneB, *ornt, *orntB, *support, *supportB;
2251:     PetscInt        coneSize, coneSizeB, c, supportSize, supportSizeB, s;

2253:     DMPlexGetConeSize(dmA, p, &coneSize);
2254:     DMPlexGetCone(dmA, p, &cone);
2255:     DMPlexGetConeOrientation(dmA, p, &ornt);
2256:     DMPlexGetConeSize(dmB, p, &coneSizeB);
2257:     DMPlexGetCone(dmB, p, &coneB);
2258:     DMPlexGetConeOrientation(dmB, p, &orntB);
2259:     if (coneSize != coneSizeB) return(0);
2260:     for (c = 0; c < coneSize; ++c) {
2261:       if (cone[c] != coneB[c]) return(0);
2262:       if (ornt[c] != orntB[c]) return(0);
2263:     }
2264:     DMPlexGetSupportSize(dmA, p, &supportSize);
2265:     DMPlexGetSupport(dmA, p, &support);
2266:     DMPlexGetSupportSize(dmB, p, &supportSizeB);
2267:     DMPlexGetSupport(dmB, p, &supportB);
2268:     if (supportSize != supportSizeB) return(0);
2269:     for (s = 0; s < supportSize; ++s) {
2270:       if (support[s] != supportB[s]) return(0);
2271:     }
2272:   }
2273:   *equal = PETSC_TRUE;
2274:   return(0);
2275: }

2279: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2280: {
2281:   MPI_Comm       comm;

2285:   PetscObjectGetComm((PetscObject)dm,&comm);
2287:   switch (cellDim) {
2288:   case 0:
2289:     *numFaceVertices = 0;
2290:     break;
2291:   case 1:
2292:     *numFaceVertices = 1;
2293:     break;
2294:   case 2:
2295:     switch (numCorners) {
2296:     case 3: /* triangle */
2297:       *numFaceVertices = 2; /* Edge has 2 vertices */
2298:       break;
2299:     case 4: /* quadrilateral */
2300:       *numFaceVertices = 2; /* Edge has 2 vertices */
2301:       break;
2302:     case 6: /* quadratic triangle, tri and quad cohesive Lagrange cells */
2303:       *numFaceVertices = 3; /* Edge has 3 vertices */
2304:       break;
2305:     case 9: /* quadratic quadrilateral, quadratic quad cohesive Lagrange cells */
2306:       *numFaceVertices = 3; /* Edge has 3 vertices */
2307:       break;
2308:     default:
2309:       SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim);
2310:     }
2311:     break;
2312:   case 3:
2313:     switch (numCorners) {
2314:     case 4: /* tetradehdron */
2315:       *numFaceVertices = 3; /* Face has 3 vertices */
2316:       break;
2317:     case 6: /* tet cohesive cells */
2318:       *numFaceVertices = 4; /* Face has 4 vertices */
2319:       break;
2320:     case 8: /* hexahedron */
2321:       *numFaceVertices = 4; /* Face has 4 vertices */
2322:       break;
2323:     case 9: /* tet cohesive Lagrange cells */
2324:       *numFaceVertices = 6; /* Face has 6 vertices */
2325:       break;
2326:     case 10: /* quadratic tetrahedron */
2327:       *numFaceVertices = 6; /* Face has 6 vertices */
2328:       break;
2329:     case 12: /* hex cohesive Lagrange cells */
2330:       *numFaceVertices = 6; /* Face has 6 vertices */
2331:       break;
2332:     case 18: /* quadratic tet cohesive Lagrange cells */
2333:       *numFaceVertices = 6; /* Face has 6 vertices */
2334:       break;
2335:     case 27: /* quadratic hexahedron, quadratic hex cohesive Lagrange cells */
2336:       *numFaceVertices = 9; /* Face has 9 vertices */
2337:       break;
2338:     default:
2339:       SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim);
2340:     }
2341:     break;
2342:   default:
2343:     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid cell dimension %d", cellDim);
2344:   }
2345:   return(0);
2346: }

2350: /* Trys to give the mesh a consistent orientation */
2351: PetscErrorCode DMPlexOrient(DM dm)
2352: {
2353:   PetscBT        seenCells, flippedCells, seenFaces;
2354:   PetscInt      *faceFIFO, fTop, fBottom;
2355:   PetscInt       dim, h, cStart, cEnd, c, fStart, fEnd, face, maxConeSize, *revcone, *revconeO;

2359:   /* Truth Table
2360:      mismatch    flips   do action   mismatch   flipA ^ flipB   action
2361:          F       0 flips     no         F             F           F
2362:          F       1 flip      yes        F             T           T
2363:          F       2 flips     no         T             F           T
2364:          T       0 flips     yes        T             T           F
2365:          T       1 flip      no
2366:          T       2 flips     yes
2367:   */
2368:   DMPlexGetDimension(dm, &dim);
2369:   DMPlexGetVTKCellHeight(dm, &h);
2370:   DMPlexGetHeightStratum(dm, h,   &cStart, &cEnd);
2371:   DMPlexGetHeightStratum(dm, h+1, &fStart, &fEnd);
2372:   PetscBTCreate(cEnd - cStart, &seenCells);
2373:   PetscBTMemzero(cEnd - cStart, seenCells);
2374:   PetscBTCreate(cEnd - cStart, &flippedCells);
2375:   PetscBTMemzero(cEnd - cStart, flippedCells);
2376:   PetscBTCreate(fEnd - fStart, &seenFaces);
2377:   PetscBTMemzero(fEnd - fStart, seenFaces);
2378:   PetscMalloc1((fEnd - fStart), &faceFIFO);
2379:   fTop = fBottom = 0;
2380:   /* Initialize FIFO with first cell */
2381:   if (cEnd > cStart) {
2382:     const PetscInt *cone;
2383:     PetscInt        coneSize;

2385:     DMPlexGetConeSize(dm, cStart, &coneSize);
2386:     DMPlexGetCone(dm, cStart, &cone);
2387:     for (c = 0; c < coneSize; ++c) {
2388:       faceFIFO[fBottom++] = cone[c];
2389:       PetscBTSet(seenFaces, cone[c]-fStart);
2390:     }
2391:   }
2392:   /* Consider each face in FIFO */
2393:   while (fTop < fBottom) {
2394:     const PetscInt *support, *coneA, *coneB, *coneOA, *coneOB;
2395:     PetscInt        supportSize, coneSizeA, coneSizeB, posA = -1, posB = -1;
2396:     PetscInt        seenA, flippedA, seenB, flippedB, mismatch;

2398:     face = faceFIFO[fTop++];
2399:     DMPlexGetSupportSize(dm, face, &supportSize);
2400:     DMPlexGetSupport(dm, face, &support);
2401:     if (supportSize < 2) continue;
2402:     if (supportSize != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %d", supportSize);
2403:     seenA    = PetscBTLookup(seenCells,    support[0]-cStart);
2404:     flippedA = PetscBTLookup(flippedCells, support[0]-cStart) ? 1 : 0;
2405:     seenB    = PetscBTLookup(seenCells,    support[1]-cStart);
2406:     flippedB = PetscBTLookup(flippedCells, support[1]-cStart) ? 1 : 0;

2408:     DMPlexGetConeSize(dm, support[0], &coneSizeA);
2409:     DMPlexGetConeSize(dm, support[1], &coneSizeB);
2410:     DMPlexGetCone(dm, support[0], &coneA);
2411:     DMPlexGetCone(dm, support[1], &coneB);
2412:     DMPlexGetConeOrientation(dm, support[0], &coneOA);
2413:     DMPlexGetConeOrientation(dm, support[1], &coneOB);
2414:     for (c = 0; c < coneSizeA; ++c) {
2415:       if (!PetscBTLookup(seenFaces, coneA[c]-fStart)) {
2416:         faceFIFO[fBottom++] = coneA[c];
2417:         PetscBTSet(seenFaces, coneA[c]-fStart);
2418:       }
2419:       if (coneA[c] == face) posA = c;
2420:       if (fBottom > fEnd-fStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %d was pushed exceeding capacity %d > %d", coneA[c], fBottom, fEnd-fStart);
2421:     }
2422:     if (posA < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d could not be located in cell %d", face, support[0]);
2423:     for (c = 0; c < coneSizeB; ++c) {
2424:       if (!PetscBTLookup(seenFaces, coneB[c]-fStart)) {
2425:         faceFIFO[fBottom++] = coneB[c];
2426:         PetscBTSet(seenFaces, coneB[c]-fStart);
2427:       }
2428:       if (coneB[c] == face) posB = c;
2429:       if (fBottom > fEnd-fStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %d was pushed exceeding capacity %d > %d", coneA[c], fBottom, fEnd-fStart);
2430:     }
2431:     if (posB < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d could not be located in cell %d", face, support[1]);

2433:     if (dim == 1) {
2434:       mismatch = posA == posB;
2435:     } else {
2436:       mismatch = coneOA[posA] == coneOB[posB];
2437:     }

2439:     if (mismatch ^ (flippedA ^ flippedB)) {
2440:       if (seenA && seenB) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen cells %d and %d do not match: Fault mesh is non-orientable", support[0], support[1]);
2441:       if (!seenA && !flippedA) {
2442:         PetscBTSet(flippedCells, support[0]-cStart);
2443:       } else if (!seenB && !flippedB) {
2444:         PetscBTSet(flippedCells, support[1]-cStart);
2445:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
2446:     } else if (mismatch && flippedA && flippedB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
2447:     PetscBTSet(seenCells, support[0]-cStart);
2448:     PetscBTSet(seenCells, support[1]-cStart);
2449:   }
2450:   /* Now all subdomains are oriented, but we need a consistent parallel orientation */
2451:   {
2452:     /* Find a representative face (edge) separating pairs of procs */
2453:     PetscSF            sf;
2454:     const PetscInt    *lpoints;
2455:     const PetscSFNode *rpoints;
2456:     PetscInt          *neighbors, *nranks;
2457:     PetscInt           numLeaves, numRoots, numNeighbors = 0, l, n;

2459:     DMGetPointSF(dm, &sf);
2460:     PetscSFGetGraph(sf, &numRoots, &numLeaves, &lpoints, &rpoints);
2461:     if (numLeaves >= 0) {
2462:       const PetscInt *cone, *ornt, *support;
2463:       PetscInt        coneSize, supportSize;
2464:       int            *rornt, *lornt; /* PetscSF cannot handle smaller than int */
2465:       PetscBool      *match, flipped = PETSC_FALSE;

2467:       PetscMalloc1(numLeaves,&neighbors);
2468:       /* I know this is p^2 time in general, but for bounded degree its alright */
2469:       for (l = 0; l < numLeaves; ++l) {
2470:         const PetscInt face = lpoints[l];
2471:         if ((face >= fStart) && (face < fEnd)) {
2472:           const PetscInt rank = rpoints[l].rank;
2473:           for (n = 0; n < numNeighbors; ++n) if (rank == rpoints[neighbors[n]].rank) break;
2474:           if (n >= numNeighbors) {
2475:             PetscInt supportSize;
2476:             DMPlexGetSupportSize(dm, face, &supportSize);
2477:             if (supportSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Boundary faces should see one cell, not %d", supportSize);
2478:             neighbors[numNeighbors++] = l;
2479:           }
2480:         }
2481:       }
2482:       PetscCalloc4(numNeighbors,&match,numNeighbors,&nranks,numRoots,&rornt,numRoots,&lornt);
2483:       for (face = fStart; face < fEnd; ++face) {
2484:         DMPlexGetSupportSize(dm, face, &supportSize);
2485:         if (supportSize != 1) continue;
2486:         DMPlexGetSupport(dm, face, &support);

2488:         DMPlexGetCone(dm, support[0], &cone);
2489:         DMPlexGetConeSize(dm, support[0], &coneSize);
2490:         DMPlexGetConeOrientation(dm, support[0], &ornt);
2491:         for (c = 0; c < coneSize; ++c) if (cone[c] == face) break;
2492:         if (dim == 1) {
2493:           /* Use cone position instead, shifted to -1 or 1 */
2494:           rornt[face] = c*2-1;
2495:         } else {
2496:           if (PetscBTLookup(flippedCells, support[0]-cStart)) rornt[face] = ornt[c] < 0 ? -1 :  1;
2497:           else                                                rornt[face] = ornt[c] < 0 ?  1 : -1;
2498:         }
2499:       }
2500:       /* Mark each edge with match or nomatch */
2501:       PetscSFBcastBegin(sf, MPI_INT, rornt, lornt);
2502:       PetscSFBcastEnd(sf, MPI_INT, rornt, lornt);
2503:       for (n = 0; n < numNeighbors; ++n) {
2504:         const PetscInt face = lpoints[neighbors[n]];

2506:         if (rornt[face]*lornt[face] < 0) match[n] = PETSC_TRUE;
2507:         else                             match[n] = PETSC_FALSE;
2508:         nranks[n] = rpoints[neighbors[n]].rank;
2509:       }
2510:       /* Collect the graph on 0 */
2511:       {
2512:         MPI_Comm     comm = PetscObjectComm((PetscObject) sf);
2513:         Mat          G;
2514:         PetscBT      seenProcs, flippedProcs;
2515:         PetscInt    *procFIFO, pTop, pBottom;
2516:         PetscInt    *adj = NULL;
2517:         PetscBool   *val = NULL;
2518:         PetscMPIInt *recvcounts = NULL, *displs = NULL, p;
2519:         PetscMPIInt  N = numNeighbors, numProcs = 0, rank;
2520:         PetscInt     debug = 0;

2522:         MPI_Comm_rank(comm, &rank);
2523:         if (!rank) {MPI_Comm_size(comm, &numProcs);}
2524:         PetscCalloc2(numProcs,&recvcounts,numProcs+1,&displs);
2525:         MPI_Gather(&N, 1, MPI_INT, recvcounts, 1, MPI_INT, 0, comm);
2526:         for (p = 0; p < numProcs; ++p) {
2527:           displs[p+1] = displs[p] + recvcounts[p];
2528:         }
2529:         if (!rank) {PetscMalloc2(displs[numProcs],&adj,displs[numProcs],&val);}
2530:         MPI_Gatherv(nranks, numNeighbors, MPIU_INT, adj, recvcounts, displs, MPIU_INT, 0, comm);
2531:         MPI_Gatherv(match, numNeighbors, MPIU_BOOL, val, recvcounts, displs, MPIU_BOOL, 0, comm);
2532:         if (debug) {
2533:           for (p = 0; p < numProcs; ++p) {
2534:             PetscPrintf(comm, "Proc %d:\n", p);
2535:             for (n = 0; n < recvcounts[p]; ++n) {
2536:               PetscPrintf(comm, "  edge %d (%d):\n", adj[displs[p]+n], val[displs[p]+n]);
2537:             }
2538:           }
2539:         }
2540:         /* Symmetrize the graph */
2541:         MatCreate(PETSC_COMM_SELF, &G);
2542:         MatSetSizes(G, numProcs, numProcs, numProcs, numProcs);
2543:         MatSetUp(G);
2544:         for (p = 0; p < numProcs; ++p) {
2545:           for (n = 0; n < recvcounts[p]; ++n) {
2546:             const PetscInt    r = p;
2547:             const PetscInt    q = adj[displs[p]+n];
2548:             const PetscScalar o = val[displs[p]+n] ? 1.0 : 0.0;

2550:             MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES);
2551:             MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES);
2552:           }
2553:         }
2554:         MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY);
2555:         MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY);

2557:         PetscBTCreate(numProcs, &seenProcs);
2558:         PetscBTMemzero(numProcs, seenProcs);
2559:         PetscBTCreate(numProcs, &flippedProcs);
2560:         PetscBTMemzero(numProcs, flippedProcs);
2561:         PetscMalloc1(numProcs,&procFIFO);
2562:         pTop = pBottom = 0;
2563:         for (p = 0; p < numProcs; ++p) {
2564:           if (PetscBTLookup(seenProcs, p)) continue;
2565:           /* Initialize FIFO with next proc */
2566:           procFIFO[pBottom++] = p;
2567:           PetscBTSet(seenProcs, p);
2568:           /* Consider each proc in FIFO */
2569:           while (pTop < pBottom) {
2570:             const PetscScalar *ornt;
2571:             const PetscInt    *neighbors;
2572:             PetscInt           proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors;

2574:             proc     = procFIFO[pTop++];
2575:             flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0;
2576:             MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt);
2577:             /* Loop over neighboring procs */
2578:             for (n = 0; n < numNeighbors; ++n) {
2579:               nproc    = neighbors[n];
2580:               mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1;
2581:               seen     = PetscBTLookup(seenProcs, nproc);
2582:               flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0;

2584:               if (mismatch ^ (flippedA ^ flippedB)) {
2585:                 if (seen) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen procs %d and %d do not match: Fault mesh is non-orientable", proc, nproc);
2586:                 if (!flippedB) {
2587:                   PetscBTSet(flippedProcs, nproc);
2588:               } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
2589:               } else if (mismatch && flippedA && flippedB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
2590:               if (!seen) {
2591:                 procFIFO[pBottom++] = nproc;
2592:                 PetscBTSet(seenProcs, nproc);
2593:               }
2594:             }
2595:           }
2596:         }
2597:         PetscFree(procFIFO);
2598:         MatDestroy(&G);

2600:         PetscFree2(recvcounts,displs);
2601:         PetscFree2(adj,val);
2602:         {
2603:           PetscBool *flips;

2605:           PetscMalloc1(numProcs,&flips);
2606:           for (p = 0; p < numProcs; ++p) {
2607:             flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE;
2608:             if (debug && flips[p]) {PetscPrintf(comm, "Flipping Proc %d:\n", p);}
2609:           }
2610:           MPI_Scatter(flips, 1, MPIU_BOOL, &flipped, 1, MPIU_BOOL, 0, comm);
2611:           PetscFree(flips);
2612:         }
2613:         PetscBTDestroy(&seenProcs);
2614:         PetscBTDestroy(&flippedProcs);
2615:       }
2616:       PetscFree4(match,nranks,rornt,lornt);
2617:       PetscFree(neighbors);
2618:       if (flipped) {for (c = cStart; c < cEnd; ++c) {PetscBTNegate(flippedCells, c-cStart);}}
2619:     }
2620:   }
2621:   /* Reverse flipped cells in the mesh */
2622:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2623:   DMGetWorkArray(dm, maxConeSize, PETSC_INT, &revcone);
2624:   DMGetWorkArray(dm, maxConeSize, PETSC_INT, &revconeO);
2625:   for (c = cStart; c < cEnd; ++c) {
2626:     const PetscInt *cone, *coneO, *support;
2627:     PetscInt        coneSize, supportSize, faceSize, cp, sp;

2629:     if (!PetscBTLookup(flippedCells, c-cStart)) continue;
2630:     DMPlexGetConeSize(dm, c, &coneSize);
2631:     DMPlexGetCone(dm, c, &cone);
2632:     DMPlexGetConeOrientation(dm, c, &coneO);
2633:     for (cp = 0; cp < coneSize; ++cp) {
2634:       const PetscInt rcp = coneSize-cp-1;

2636:       DMPlexGetConeSize(dm, cone[rcp], &faceSize);
2637:       revcone[cp]  = cone[rcp];
2638:       revconeO[cp] = coneO[rcp] >= 0 ? -(faceSize-coneO[rcp]) : faceSize+coneO[rcp];
2639:     }
2640:     DMPlexSetCone(dm, c, revcone);
2641:     DMPlexSetConeOrientation(dm, c, revconeO);
2642:     /* Reverse orientations of support */
2643:     faceSize = coneSize;
2644:     DMPlexGetSupportSize(dm, c, &supportSize);
2645:     DMPlexGetSupport(dm, c, &support);
2646:     for (sp = 0; sp < supportSize; ++sp) {
2647:       DMPlexGetConeSize(dm, support[sp], &coneSize);
2648:       DMPlexGetCone(dm, support[sp], &cone);
2649:       DMPlexGetConeOrientation(dm, support[sp], &coneO);
2650:       for (cp = 0; cp < coneSize; ++cp) {
2651:         if (cone[cp] != c) continue;
2652:         DMPlexInsertConeOrientation(dm, support[sp], cp, coneO[cp] >= 0 ? -(faceSize-coneO[cp]) : faceSize+coneO[cp]);
2653:       }
2654:     }
2655:   }
2656:   DMRestoreWorkArray(dm, maxConeSize, PETSC_INT, &revcone);
2657:   DMRestoreWorkArray(dm, maxConeSize, PETSC_INT, &revconeO);
2658:   PetscBTDestroy(&seenCells);
2659:   PetscBTDestroy(&flippedCells);
2660:   PetscBTDestroy(&seenFaces);
2661:   PetscFree(faceFIFO);
2662:   return(0);
2663: }

2667: PetscErrorCode DMPlexInvertCell_Internal(PetscInt dim, PetscInt numCorners, PetscInt cone[])
2668: {
2669:   int tmpc;

2672:   if (dim != 3) return(0);
2673:   switch (numCorners) {
2674:   case 4:
2675:     tmpc    = cone[0];
2676:     cone[0] = cone[1];
2677:     cone[1] = tmpc;
2678:     break;
2679:   case 8:
2680:     tmpc    = cone[1];
2681:     cone[1] = cone[3];
2682:     cone[3] = tmpc;
2683:     break;
2684:   default: break;
2685:   }
2686:   return(0);
2687: }

2691: /*@C
2692:   DMPlexInvertCell - This flips tetrahedron and hexahedron orientation since Plex stores them internally with outward normals. Other cells are left untouched.

2694:   Input Parameters:
2695: + numCorners - The number of vertices in a cell
2696: - cone - The incoming cone

2698:   Output Parameter:
2699: . cone - The inverted cone (in-place)

2701:   Level: developer

2703: .seealso: DMPlexGenerate()
2704: @*/
2705: PetscErrorCode DMPlexInvertCell(PetscInt dim, PetscInt numCorners, int cone[])
2706: {
2707:   int tmpc;

2710:   if (dim != 3) return(0);
2711:   switch (numCorners) {
2712:   case 4:
2713:     tmpc    = cone[0];
2714:     cone[0] = cone[1];
2715:     cone[1] = tmpc;
2716:     break;
2717:   case 8:
2718:     tmpc    = cone[1];
2719:     cone[1] = cone[3];
2720:     cone[3] = tmpc;
2721:     break;
2722:   default: break;
2723:   }
2724:   return(0);
2725: }

2729: /* This is to fix the tetrahedron orientation from TetGen */
2730: PETSC_UNUSED static PetscErrorCode DMPlexInvertCells_Internal(PetscInt dim, PetscInt numCells, PetscInt numCorners, int cells[])
2731: {
2732:   PetscInt       bound = numCells*numCorners, coff;

2736:   for (coff = 0; coff < bound; coff += numCorners) {
2737:     DMPlexInvertCell(dim, numCorners, &cells[coff]);
2738:   }
2739:   return(0);
2740: }

2742: #if defined(PETSC_HAVE_TRIANGLE)
2743: #include <triangle.h>

2747: PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx)
2748: {
2750:   inputCtx->numberofpoints             = 0;
2751:   inputCtx->numberofpointattributes    = 0;
2752:   inputCtx->pointlist                  = NULL;
2753:   inputCtx->pointattributelist         = NULL;
2754:   inputCtx->pointmarkerlist            = NULL;
2755:   inputCtx->numberofsegments           = 0;
2756:   inputCtx->segmentlist                = NULL;
2757:   inputCtx->segmentmarkerlist          = NULL;
2758:   inputCtx->numberoftriangleattributes = 0;
2759:   inputCtx->trianglelist               = NULL;
2760:   inputCtx->numberofholes              = 0;
2761:   inputCtx->holelist                   = NULL;
2762:   inputCtx->numberofregions            = 0;
2763:   inputCtx->regionlist                 = NULL;
2764:   return(0);
2765: }

2769: PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
2770: {
2772:   outputCtx->numberofpoints        = 0;
2773:   outputCtx->pointlist             = NULL;
2774:   outputCtx->pointattributelist    = NULL;
2775:   outputCtx->pointmarkerlist       = NULL;
2776:   outputCtx->numberoftriangles     = 0;
2777:   outputCtx->trianglelist          = NULL;
2778:   outputCtx->triangleattributelist = NULL;
2779:   outputCtx->neighborlist          = NULL;
2780:   outputCtx->segmentlist           = NULL;
2781:   outputCtx->segmentmarkerlist     = NULL;
2782:   outputCtx->numberofedges         = 0;
2783:   outputCtx->edgelist              = NULL;
2784:   outputCtx->edgemarkerlist        = NULL;
2785:   return(0);
2786: }

2790: PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
2791: {
2793:   free(outputCtx->pointlist);
2794:   free(outputCtx->pointmarkerlist);
2795:   free(outputCtx->segmentlist);
2796:   free(outputCtx->segmentmarkerlist);
2797:   free(outputCtx->edgelist);
2798:   free(outputCtx->edgemarkerlist);
2799:   free(outputCtx->trianglelist);
2800:   free(outputCtx->neighborlist);
2801:   return(0);
2802: }

2806: PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
2807: {
2808:   MPI_Comm             comm;
2809:   PetscInt             dim              = 2;
2810:   const PetscBool      createConvexHull = PETSC_FALSE;
2811:   const PetscBool      constrained      = PETSC_FALSE;
2812:   struct triangulateio in;
2813:   struct triangulateio out;
2814:   PetscInt             vStart, vEnd, v, eStart, eEnd, e;
2815:   PetscMPIInt          rank;
2816:   PetscErrorCode       ierr;

2819:   PetscObjectGetComm((PetscObject)boundary,&comm);
2820:   MPI_Comm_rank(comm, &rank);
2821:   InitInput_Triangle(&in);
2822:   InitOutput_Triangle(&out);
2823:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);

2825:   in.numberofpoints = vEnd - vStart;
2826:   if (in.numberofpoints > 0) {
2827:     PetscSection coordSection;
2828:     Vec          coordinates;
2829:     PetscScalar *array;

2831:     PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
2832:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
2833:     DMGetCoordinatesLocal(boundary, &coordinates);
2834:     DMGetCoordinateSection(boundary, &coordSection);
2835:     VecGetArray(coordinates, &array);
2836:     for (v = vStart; v < vEnd; ++v) {
2837:       const PetscInt idx = v - vStart;
2838:       PetscInt       off, d;

2840:       PetscSectionGetOffset(coordSection, v, &off);
2841:       for (d = 0; d < dim; ++d) {
2842:         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
2843:       }
2844:       DMPlexGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);
2845:     }
2846:     VecRestoreArray(coordinates, &array);
2847:   }
2848:   DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
2849:   in.numberofsegments = eEnd - eStart;
2850:   if (in.numberofsegments > 0) {
2851:     PetscMalloc1(in.numberofsegments*2, &in.segmentlist);
2852:     PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);
2853:     for (e = eStart; e < eEnd; ++e) {
2854:       const PetscInt  idx = e - eStart;
2855:       const PetscInt *cone;

2857:       DMPlexGetCone(boundary, e, &cone);

2859:       in.segmentlist[idx*2+0] = cone[0] - vStart;
2860:       in.segmentlist[idx*2+1] = cone[1] - vStart;

2862:       DMPlexGetLabelValue(boundary, "marker", e, &in.segmentmarkerlist[idx]);
2863:     }
2864:   }
2865: #if 0 /* Do not currently support holes */
2866:   PetscReal *holeCoords;
2867:   PetscInt   h, d;

2869:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
2870:   if (in.numberofholes > 0) {
2871:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
2872:     for (h = 0; h < in.numberofholes; ++h) {
2873:       for (d = 0; d < dim; ++d) {
2874:         in.holelist[h*dim+d] = holeCoords[h*dim+d];
2875:       }
2876:     }
2877:   }
2878: #endif
2879:   if (!rank) {
2880:     char args[32];

2882:     /* Take away 'Q' for verbose output */
2883:     PetscStrcpy(args, "pqezQ");
2884:     if (createConvexHull) {
2885:       PetscStrcat(args, "c");
2886:     }
2887:     if (constrained) {
2888:       PetscStrcpy(args, "zepDQ");
2889:     }
2890:     triangulate(args, &in, &out, NULL);
2891:   }
2892:   PetscFree(in.pointlist);
2893:   PetscFree(in.pointmarkerlist);
2894:   PetscFree(in.segmentlist);
2895:   PetscFree(in.segmentmarkerlist);
2896:   PetscFree(in.holelist);

2898:   {
2899:     const PetscInt numCorners  = 3;
2900:     const PetscInt numCells    = out.numberoftriangles;
2901:     const PetscInt numVertices = out.numberofpoints;
2902:     const int     *cells      = out.trianglelist;
2903:     const double  *meshCoords = out.pointlist;

2905:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
2906:     /* Set labels */
2907:     for (v = 0; v < numVertices; ++v) {
2908:       if (out.pointmarkerlist[v]) {
2909:         DMPlexSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);
2910:       }
2911:     }
2912:     if (interpolate) {
2913:       for (e = 0; e < out.numberofedges; e++) {
2914:         if (out.edgemarkerlist[e]) {
2915:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
2916:           const PetscInt *edges;
2917:           PetscInt        numEdges;

2919:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
2920:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
2921:           DMPlexSetLabelValue(*dm, "marker", edges[0], out.edgemarkerlist[e]);
2922:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
2923:         }
2924:       }
2925:     }
2926:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
2927:   }
2928: #if 0 /* Do not currently support holes */
2929:   DMPlexCopyHoles(*dm, boundary);
2930: #endif
2931:   FiniOutput_Triangle(&out);
2932:   return(0);
2933: }

2937: PetscErrorCode DMPlexRefine_Triangle(DM dm, double *maxVolumes, DM *dmRefined)
2938: {
2939:   MPI_Comm             comm;
2940:   PetscInt             dim  = 2;
2941:   struct triangulateio in;
2942:   struct triangulateio out;
2943:   PetscInt             vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
2944:   PetscMPIInt          rank;
2945:   PetscErrorCode       ierr;

2948:   PetscObjectGetComm((PetscObject)dm,&comm);
2949:   MPI_Comm_rank(comm, &rank);
2950:   InitInput_Triangle(&in);
2951:   InitOutput_Triangle(&out);
2952:   DMPlexGetDepth(dm, &depth);
2953:   MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
2954:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);

2956:   in.numberofpoints = vEnd - vStart;
2957:   if (in.numberofpoints > 0) {
2958:     PetscSection coordSection;
2959:     Vec          coordinates;
2960:     PetscScalar *array;

2962:     PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
2963:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
2964:     DMGetCoordinatesLocal(dm, &coordinates);
2965:     DMGetCoordinateSection(dm, &coordSection);
2966:     VecGetArray(coordinates, &array);
2967:     for (v = vStart; v < vEnd; ++v) {
2968:       const PetscInt idx = v - vStart;
2969:       PetscInt       off, d;

2971:       PetscSectionGetOffset(coordSection, v, &off);
2972:       for (d = 0; d < dim; ++d) {
2973:         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
2974:       }
2975:       DMPlexGetLabelValue(dm, "marker", v, &in.pointmarkerlist[idx]);
2976:     }
2977:     VecRestoreArray(coordinates, &array);
2978:   }
2979:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

2981:   in.numberofcorners   = 3;
2982:   in.numberoftriangles = cEnd - cStart;

2984:   in.trianglearealist  = (double*) maxVolumes;
2985:   if (in.numberoftriangles > 0) {
2986:     PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
2987:     for (c = cStart; c < cEnd; ++c) {
2988:       const PetscInt idx      = c - cStart;
2989:       PetscInt      *closure = NULL;
2990:       PetscInt       closureSize;

2992:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
2993:       if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
2994:       for (v = 0; v < 3; ++v) {
2995:         in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
2996:       }
2997:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
2998:     }
2999:   }
3000:   /* TODO: Segment markers are missing on input */
3001: #if 0 /* Do not currently support holes */
3002:   PetscReal *holeCoords;
3003:   PetscInt   h, d;

3005:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
3006:   if (in.numberofholes > 0) {
3007:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
3008:     for (h = 0; h < in.numberofholes; ++h) {
3009:       for (d = 0; d < dim; ++d) {
3010:         in.holelist[h*dim+d] = holeCoords[h*dim+d];
3011:       }
3012:     }
3013:   }
3014: #endif
3015:   if (!rank) {
3016:     char args[32];

3018:     /* Take away 'Q' for verbose output */
3019:     PetscStrcpy(args, "pqezQra");
3020:     triangulate(args, &in, &out, NULL);
3021:   }
3022:   PetscFree(in.pointlist);
3023:   PetscFree(in.pointmarkerlist);
3024:   PetscFree(in.segmentlist);
3025:   PetscFree(in.segmentmarkerlist);
3026:   PetscFree(in.trianglelist);

3028:   {
3029:     const PetscInt numCorners  = 3;
3030:     const PetscInt numCells    = out.numberoftriangles;
3031:     const PetscInt numVertices = out.numberofpoints;
3032:     const int     *cells      = out.trianglelist;
3033:     const double  *meshCoords = out.pointlist;
3034:     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;

3036:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
3037:     /* Set labels */
3038:     for (v = 0; v < numVertices; ++v) {
3039:       if (out.pointmarkerlist[v]) {
3040:         DMPlexSetLabelValue(*dmRefined, "marker", v+numCells, out.pointmarkerlist[v]);
3041:       }
3042:     }
3043:     if (interpolate) {
3044:       PetscInt e;

3046:       for (e = 0; e < out.numberofedges; e++) {
3047:         if (out.edgemarkerlist[e]) {
3048:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3049:           const PetscInt *edges;
3050:           PetscInt        numEdges;

3052:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3053:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3054:           DMPlexSetLabelValue(*dmRefined, "marker", edges[0], out.edgemarkerlist[e]);
3055:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3056:         }
3057:       }
3058:     }
3059:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
3060:   }
3061: #if 0 /* Do not currently support holes */
3062:   DMPlexCopyHoles(*dm, boundary);
3063: #endif
3064:   FiniOutput_Triangle(&out);
3065:   return(0);
3066: }
3067: #endif

3069: #if defined(PETSC_HAVE_TETGEN)
3070: #include <tetgen.h>
3073: PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
3074: {
3075:   MPI_Comm       comm;
3076:   const PetscInt dim  = 3;
3077:   ::tetgenio     in;
3078:   ::tetgenio     out;
3079:   PetscInt       vStart, vEnd, v, fStart, fEnd, f;
3080:   PetscMPIInt    rank;

3084:   PetscObjectGetComm((PetscObject)boundary,&comm);
3085:   MPI_Comm_rank(comm, &rank);
3086:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
3087:   in.numberofpoints = vEnd - vStart;
3088:   if (in.numberofpoints > 0) {
3089:     PetscSection coordSection;
3090:     Vec          coordinates;
3091:     PetscScalar *array;

3093:     in.pointlist       = new double[in.numberofpoints*dim];
3094:     in.pointmarkerlist = new int[in.numberofpoints];

3096:     DMGetCoordinatesLocal(boundary, &coordinates);
3097:     DMGetCoordinateSection(boundary, &coordSection);
3098:     VecGetArray(coordinates, &array);
3099:     for (v = vStart; v < vEnd; ++v) {
3100:       const PetscInt idx = v - vStart;
3101:       PetscInt       off, d;

3103:       PetscSectionGetOffset(coordSection, v, &off);
3104:       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3105:       DMPlexGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);
3106:     }
3107:     VecRestoreArray(coordinates, &array);
3108:   }
3109:   DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);

3111:   in.numberoffacets = fEnd - fStart;
3112:   if (in.numberoffacets > 0) {
3113:     in.facetlist       = new tetgenio::facet[in.numberoffacets];
3114:     in.facetmarkerlist = new int[in.numberoffacets];
3115:     for (f = fStart; f < fEnd; ++f) {
3116:       const PetscInt idx     = f - fStart;
3117:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v;

3119:       in.facetlist[idx].numberofpolygons = 1;
3120:       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
3121:       in.facetlist[idx].numberofholes    = 0;
3122:       in.facetlist[idx].holelist         = NULL;

3124:       DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3125:       for (p = 0; p < numPoints*2; p += 2) {
3126:         const PetscInt point = points[p];
3127:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
3128:       }

3130:       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
3131:       poly->numberofvertices = numVertices;
3132:       poly->vertexlist       = new int[poly->numberofvertices];
3133:       for (v = 0; v < numVertices; ++v) {
3134:         const PetscInt vIdx = points[v] - vStart;
3135:         poly->vertexlist[v] = vIdx;
3136:       }
3137:       DMPlexGetLabelValue(boundary, "marker", f, &in.facetmarkerlist[idx]);
3138:       DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3139:     }
3140:   }
3141:   if (!rank) {
3142:     char args[32];

3144:     /* Take away 'Q' for verbose output */
3145:     PetscStrcpy(args, "pqezQ");
3146:     ::tetrahedralize(args, &in, &out);
3147:   }
3148:   {
3149:     const PetscInt numCorners  = 4;
3150:     const PetscInt numCells    = out.numberoftetrahedra;
3151:     const PetscInt numVertices = out.numberofpoints;
3152:     const double   *meshCoords = out.pointlist;
3153:     int            *cells      = out.tetrahedronlist;

3155:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
3156:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
3157:     /* Set labels */
3158:     for (v = 0; v < numVertices; ++v) {
3159:       if (out.pointmarkerlist[v]) {
3160:         DMPlexSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);
3161:       }
3162:     }
3163:     if (interpolate) {
3164:       PetscInt e;

3166:       for (e = 0; e < out.numberofedges; e++) {
3167:         if (out.edgemarkerlist[e]) {
3168:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3169:           const PetscInt *edges;
3170:           PetscInt        numEdges;

3172:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
3173:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3174:           DMPlexSetLabelValue(*dm, "marker", edges[0], out.edgemarkerlist[e]);
3175:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
3176:         }
3177:       }
3178:       for (f = 0; f < out.numberoftrifaces; f++) {
3179:         if (out.trifacemarkerlist[f]) {
3180:           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
3181:           const PetscInt *faces;
3182:           PetscInt        numFaces;

3184:           DMPlexGetJoin(*dm, 3, vertices, &numFaces, &faces);
3185:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3186:           DMPlexSetLabelValue(*dm, "marker", faces[0], out.trifacemarkerlist[f]);
3187:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
3188:         }
3189:       }
3190:     }
3191:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
3192:   }
3193:   return(0);
3194: }

3198: PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
3199: {
3200:   MPI_Comm       comm;
3201:   const PetscInt dim  = 3;
3202:   ::tetgenio     in;
3203:   ::tetgenio     out;
3204:   PetscInt       vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
3205:   PetscMPIInt    rank;

3209:   PetscObjectGetComm((PetscObject)dm,&comm);
3210:   MPI_Comm_rank(comm, &rank);
3211:   DMPlexGetDepth(dm, &depth);
3212:   MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
3213:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);

3215:   in.numberofpoints = vEnd - vStart;
3216:   if (in.numberofpoints > 0) {
3217:     PetscSection coordSection;
3218:     Vec          coordinates;
3219:     PetscScalar *array;

3221:     in.pointlist       = new double[in.numberofpoints*dim];
3222:     in.pointmarkerlist = new int[in.numberofpoints];

3224:     DMGetCoordinatesLocal(dm, &coordinates);
3225:     DMGetCoordinateSection(dm, &coordSection);
3226:     VecGetArray(coordinates, &array);
3227:     for (v = vStart; v < vEnd; ++v) {
3228:       const PetscInt idx = v - vStart;
3229:       PetscInt       off, d;

3231:       PetscSectionGetOffset(coordSection, v, &off);
3232:       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3233:       DMPlexGetLabelValue(dm, "marker", v, &in.pointmarkerlist[idx]);
3234:     }
3235:     VecRestoreArray(coordinates, &array);
3236:   }
3237:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

3239:   in.numberofcorners       = 4;
3240:   in.numberoftetrahedra    = cEnd - cStart;
3241:   in.tetrahedronvolumelist = (double*) maxVolumes;
3242:   if (in.numberoftetrahedra > 0) {
3243:     in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
3244:     for (c = cStart; c < cEnd; ++c) {
3245:       const PetscInt idx      = c - cStart;
3246:       PetscInt      *closure = NULL;
3247:       PetscInt       closureSize;

3249:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3250:       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
3251:       for (v = 0; v < 4; ++v) {
3252:         in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
3253:       }
3254:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3255:     }
3256:   }
3257:   /* TODO: Put in boundary faces with markers */
3258:   if (!rank) {
3259:     char args[32];

3261:     /* Take away 'Q' for verbose output */
3262:     /*PetscStrcpy(args, "qezQra"); */
3263:     PetscStrcpy(args, "qezraVVVV");
3264:     ::tetrahedralize(args, &in, &out);
3265:   }
3266:   in.tetrahedronvolumelist = NULL;

3268:   {
3269:     const PetscInt numCorners  = 4;
3270:     const PetscInt numCells    = out.numberoftetrahedra;
3271:     const PetscInt numVertices = out.numberofpoints;
3272:     const double   *meshCoords = out.pointlist;
3273:     int            *cells      = out.tetrahedronlist;

3275:     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;

3277:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
3278:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
3279:     /* Set labels */
3280:     for (v = 0; v < numVertices; ++v) {
3281:       if (out.pointmarkerlist[v]) {
3282:         DMPlexSetLabelValue(*dmRefined, "marker", v+numCells, out.pointmarkerlist[v]);
3283:       }
3284:     }
3285:     if (interpolate) {
3286:       PetscInt e, f;

3288:       for (e = 0; e < out.numberofedges; e++) {
3289:         if (out.edgemarkerlist[e]) {
3290:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3291:           const PetscInt *edges;
3292:           PetscInt        numEdges;

3294:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3295:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3296:           DMPlexSetLabelValue(*dmRefined, "marker", edges[0], out.edgemarkerlist[e]);
3297:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3298:         }
3299:       }
3300:       for (f = 0; f < out.numberoftrifaces; f++) {
3301:         if (out.trifacemarkerlist[f]) {
3302:           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
3303:           const PetscInt *faces;
3304:           PetscInt        numFaces;

3306:           DMPlexGetJoin(*dmRefined, 3, vertices, &numFaces, &faces);
3307:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3308:           DMPlexSetLabelValue(*dmRefined, "marker", faces[0], out.trifacemarkerlist[f]);
3309:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
3310:         }
3311:       }
3312:     }
3313:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
3314:   }
3315:   return(0);
3316: }
3317: #endif

3319: #if defined(PETSC_HAVE_CTETGEN)
3320: #include <ctetgen.h>

3324: PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
3325: {
3326:   MPI_Comm       comm;
3327:   const PetscInt dim  = 3;
3328:   PLC           *in, *out;
3329:   PetscInt       verbose = 0, vStart, vEnd, v, fStart, fEnd, f;
3330:   PetscMPIInt    rank;

3334:   PetscObjectGetComm((PetscObject)boundary,&comm);
3335:   PetscOptionsGetInt(((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
3336:   MPI_Comm_rank(comm, &rank);
3337:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
3338:   PLCCreate(&in);
3339:   PLCCreate(&out);

3341:   in->numberofpoints = vEnd - vStart;
3342:   if (in->numberofpoints > 0) {
3343:     PetscSection coordSection;
3344:     Vec          coordinates;
3345:     PetscScalar *array;

3347:     PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
3348:     PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);
3349:     DMGetCoordinatesLocal(boundary, &coordinates);
3350:     DMGetCoordinateSection(boundary, &coordSection);
3351:     VecGetArray(coordinates, &array);
3352:     for (v = vStart; v < vEnd; ++v) {
3353:       const PetscInt idx = v - vStart;
3354:       PetscInt       off, d, m;

3356:       PetscSectionGetOffset(coordSection, v, &off);
3357:       for (d = 0; d < dim; ++d) {
3358:         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3359:       }
3360:       DMPlexGetLabelValue(boundary, "marker", v, &m);

3362:       in->pointmarkerlist[idx] = (int) m;
3363:     }
3364:     VecRestoreArray(coordinates, &array);
3365:   }
3366:   DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);

3368:   in->numberoffacets = fEnd - fStart;
3369:   if (in->numberoffacets > 0) {
3370:     PetscMalloc1(in->numberoffacets, &in->facetlist);
3371:     PetscMalloc1(in->numberoffacets,   &in->facetmarkerlist);
3372:     for (f = fStart; f < fEnd; ++f) {
3373:       const PetscInt idx     = f - fStart;
3374:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, m;
3375:       polygon       *poly;

3377:       in->facetlist[idx].numberofpolygons = 1;

3379:       PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);

3381:       in->facetlist[idx].numberofholes    = 0;
3382:       in->facetlist[idx].holelist         = NULL;

3384:       DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3385:       for (p = 0; p < numPoints*2; p += 2) {
3386:         const PetscInt point = points[p];
3387:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
3388:       }

3390:       poly                   = in->facetlist[idx].polygonlist;
3391:       poly->numberofvertices = numVertices;
3392:       PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
3393:       for (v = 0; v < numVertices; ++v) {
3394:         const PetscInt vIdx = points[v] - vStart;
3395:         poly->vertexlist[v] = vIdx;
3396:       }
3397:       DMPlexGetLabelValue(boundary, "marker", f, &m);
3398:       in->facetmarkerlist[idx] = (int) m;
3399:       DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3400:     }
3401:   }
3402:   if (!rank) {
3403:     TetGenOpts t;

3405:     TetGenOptsInitialize(&t);
3406:     t.in        = boundary; /* Should go away */
3407:     t.plc       = 1;
3408:     t.quality   = 1;
3409:     t.edgesout  = 1;
3410:     t.zeroindex = 1;
3411:     t.quiet     = 1;
3412:     t.verbose   = verbose;
3413:     TetGenCheckOpts(&t);
3414:     TetGenTetrahedralize(&t, in, out);
3415:   }
3416:   {
3417:     const PetscInt numCorners  = 4;
3418:     const PetscInt numCells    = out->numberoftetrahedra;
3419:     const PetscInt numVertices = out->numberofpoints;
3420:     const double   *meshCoords = out->pointlist;
3421:     int            *cells      = out->tetrahedronlist;

3423:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
3424:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
3425:     /* Set labels */
3426:     for (v = 0; v < numVertices; ++v) {
3427:       if (out->pointmarkerlist[v]) {
3428:         DMPlexSetLabelValue(*dm, "marker", v+numCells, out->pointmarkerlist[v]);
3429:       }
3430:     }
3431:     if (interpolate) {
3432:       PetscInt e;

3434:       for (e = 0; e < out->numberofedges; e++) {
3435:         if (out->edgemarkerlist[e]) {
3436:           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
3437:           const PetscInt *edges;
3438:           PetscInt        numEdges;

3440:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
3441:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3442:           DMPlexSetLabelValue(*dm, "marker", edges[0], out->edgemarkerlist[e]);
3443:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
3444:         }
3445:       }
3446:       for (f = 0; f < out->numberoftrifaces; f++) {
3447:         if (out->trifacemarkerlist[f]) {
3448:           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
3449:           const PetscInt *faces;
3450:           PetscInt        numFaces;

3452:           DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
3453:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3454:           DMPlexSetLabelValue(*dm, "marker", faces[0], out->trifacemarkerlist[f]);
3455:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
3456:         }
3457:       }
3458:     }
3459:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
3460:   }

3462:   PLCDestroy(&in);
3463:   PLCDestroy(&out);
3464:   return(0);
3465: }

3469: PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
3470: {
3471:   MPI_Comm       comm;
3472:   const PetscInt dim  = 3;
3473:   PLC           *in, *out;
3474:   PetscInt       verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
3475:   PetscMPIInt    rank;

3479:   PetscObjectGetComm((PetscObject)dm,&comm);
3480:   PetscOptionsGetInt(((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
3481:   MPI_Comm_rank(comm, &rank);
3482:   DMPlexGetDepth(dm, &depth);
3483:   MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
3484:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3485:   PLCCreate(&in);
3486:   PLCCreate(&out);

3488:   in->numberofpoints = vEnd - vStart;
3489:   if (in->numberofpoints > 0) {
3490:     PetscSection coordSection;
3491:     Vec          coordinates;
3492:     PetscScalar *array;

3494:     PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
3495:     PetscMalloc1(in->numberofpoints,       &in->pointmarkerlist);
3496:     DMGetCoordinatesLocal(dm, &coordinates);
3497:     DMGetCoordinateSection(dm, &coordSection);
3498:     VecGetArray(coordinates, &array);
3499:     for (v = vStart; v < vEnd; ++v) {
3500:       const PetscInt idx = v - vStart;
3501:       PetscInt       off, d, m;

3503:       PetscSectionGetOffset(coordSection, v, &off);
3504:       for (d = 0; d < dim; ++d) {
3505:         in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3506:       }
3507:       DMPlexGetLabelValue(dm, "marker", v, &m);

3509:       in->pointmarkerlist[idx] = (int) m;
3510:     }
3511:     VecRestoreArray(coordinates, &array);
3512:   }
3513:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

3515:   in->numberofcorners       = 4;
3516:   in->numberoftetrahedra    = cEnd - cStart;
3517:   in->tetrahedronvolumelist = maxVolumes;
3518:   if (in->numberoftetrahedra > 0) {
3519:     PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
3520:     for (c = cStart; c < cEnd; ++c) {
3521:       const PetscInt idx      = c - cStart;
3522:       PetscInt      *closure = NULL;
3523:       PetscInt       closureSize;

3525:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3526:       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
3527:       for (v = 0; v < 4; ++v) {
3528:         in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
3529:       }
3530:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3531:     }
3532:   }
3533:   if (!rank) {
3534:     TetGenOpts t;

3536:     TetGenOptsInitialize(&t);

3538:     t.in        = dm; /* Should go away */
3539:     t.refine    = 1;
3540:     t.varvolume = 1;
3541:     t.quality   = 1;
3542:     t.edgesout  = 1;
3543:     t.zeroindex = 1;
3544:     t.quiet     = 1;
3545:     t.verbose   = verbose; /* Change this */

3547:     TetGenCheckOpts(&t);
3548:     TetGenTetrahedralize(&t, in, out);
3549:   }
3550:   {
3551:     const PetscInt numCorners  = 4;
3552:     const PetscInt numCells    = out->numberoftetrahedra;
3553:     const PetscInt numVertices = out->numberofpoints;
3554:     const double   *meshCoords = out->pointlist;
3555:     int            *cells      = out->tetrahedronlist;
3556:     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;

3558:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
3559:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
3560:     /* Set labels */
3561:     for (v = 0; v < numVertices; ++v) {
3562:       if (out->pointmarkerlist[v]) {
3563:         DMPlexSetLabelValue(*dmRefined, "marker", v+numCells, out->pointmarkerlist[v]);
3564:       }
3565:     }
3566:     if (interpolate) {
3567:       PetscInt e, f;

3569:       for (e = 0; e < out->numberofedges; e++) {
3570:         if (out->edgemarkerlist[e]) {
3571:           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
3572:           const PetscInt *edges;
3573:           PetscInt        numEdges;

3575:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3576:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3577:           DMPlexSetLabelValue(*dmRefined, "marker", edges[0], out->edgemarkerlist[e]);
3578:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3579:         }
3580:       }
3581:       for (f = 0; f < out->numberoftrifaces; f++) {
3582:         if (out->trifacemarkerlist[f]) {
3583:           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
3584:           const PetscInt *faces;
3585:           PetscInt        numFaces;

3587:           DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
3588:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3589:           DMPlexSetLabelValue(*dmRefined, "marker", faces[0], out->trifacemarkerlist[f]);
3590:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
3591:         }
3592:       }
3593:     }
3594:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
3595:   }
3596:   PLCDestroy(&in);
3597:   PLCDestroy(&out);
3598:   return(0);
3599: }
3600: #endif

3604: /*@C
3605:   DMPlexGenerate - Generates a mesh.

3607:   Not Collective

3609:   Input Parameters:
3610: + boundary - The DMPlex boundary object
3611: . name - The mesh generation package name
3612: - interpolate - Flag to create intermediate mesh elements

3614:   Output Parameter:
3615: . mesh - The DMPlex object

3617:   Level: intermediate

3619: .keywords: mesh, elements
3620: .seealso: DMPlexCreate(), DMRefine()
3621: @*/
3622: PetscErrorCode DMPlexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
3623: {
3624:   PetscInt       dim;
3625:   char           genname[1024];
3626:   PetscBool      isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;

3632:   DMPlexGetDimension(boundary, &dim);
3633:   PetscOptionsGetString(((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);
3634:   if (flg) name = genname;
3635:   if (name) {
3636:     PetscStrcmp(name, "triangle", &isTriangle);
3637:     PetscStrcmp(name, "tetgen",   &isTetgen);
3638:     PetscStrcmp(name, "ctetgen",  &isCTetgen);
3639:   }
3640:   switch (dim) {
3641:   case 1:
3642:     if (!name || isTriangle) {
3643: #if defined(PETSC_HAVE_TRIANGLE)
3644:       DMPlexGenerate_Triangle(boundary, interpolate, mesh);
3645: #else
3646:       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
3647: #endif
3648:     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
3649:     break;
3650:   case 2:
3651:     if (!name || isCTetgen) {
3652: #if defined(PETSC_HAVE_CTETGEN)
3653:       DMPlexGenerate_CTetgen(boundary, interpolate, mesh);
3654: #else
3655:       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
3656: #endif
3657:     } else if (isTetgen) {
3658: #if defined(PETSC_HAVE_TETGEN)
3659:       DMPlexGenerate_Tetgen(boundary, interpolate, mesh);
3660: #else
3661:       SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
3662: #endif
3663:     } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
3664:     break;
3665:   default:
3666:     SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
3667:   }
3668:   return(0);
3669: }

3673: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
3674: {
3675:   PetscReal      refinementLimit;
3676:   PetscInt       dim, cStart, cEnd;
3677:   char           genname[1024], *name = NULL;
3678:   PetscBool      isUniform, isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;

3682:   DMPlexGetRefinementUniform(dm, &isUniform);
3683:   if (isUniform) {
3684:     CellRefiner cellRefiner;

3686:     DMPlexGetCellRefiner_Internal(dm, &cellRefiner);
3687:     DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);
3688:     DMPlexCopyBoundary(dm, *dmRefined);
3689:     return(0);
3690:   }
3691:   DMPlexGetRefinementLimit(dm, &refinementLimit);
3692:   if (refinementLimit == 0.0) return(0);
3693:   DMPlexGetDimension(dm, &dim);
3694:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
3695:   PetscOptionsGetString(((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);
3696:   if (flg) name = genname;
3697:   if (name) {
3698:     PetscStrcmp(name, "triangle", &isTriangle);
3699:     PetscStrcmp(name, "tetgen",   &isTetgen);
3700:     PetscStrcmp(name, "ctetgen",  &isCTetgen);
3701:   }
3702:   switch (dim) {
3703:   case 2:
3704:     if (!name || isTriangle) {
3705: #if defined(PETSC_HAVE_TRIANGLE)
3706:       double  *maxVolumes;
3707:       PetscInt c;

3709:       PetscMalloc1((cEnd - cStart), &maxVolumes);
3710:       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
3711:       DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);
3712:       PetscFree(maxVolumes);
3713: #else
3714:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
3715: #endif
3716:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
3717:     break;
3718:   case 3:
3719:     if (!name || isCTetgen) {
3720: #if defined(PETSC_HAVE_CTETGEN)
3721:       PetscReal *maxVolumes;
3722:       PetscInt   c;

3724:       PetscMalloc1((cEnd - cStart), &maxVolumes);
3725:       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
3726:       DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);
3727: #else
3728:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
3729: #endif
3730:     } else if (isTetgen) {
3731: #if defined(PETSC_HAVE_TETGEN)
3732:       double  *maxVolumes;
3733:       PetscInt c;

3735:       PetscMalloc1((cEnd - cStart), &maxVolumes);
3736:       for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
3737:       DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);
3738: #else
3739:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
3740: #endif
3741:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
3742:     break;
3743:   default:
3744:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
3745:   }
3746:   DMPlexCopyBoundary(dm, *dmRefined);
3747:   return(0);
3748: }

3752: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
3753: {
3754:   DM             cdm = dm;
3755:   PetscInt       r;
3756:   PetscBool      isUniform;

3760:   DMPlexGetRefinementUniform(dm, &isUniform);
3761:   if (!isUniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Non-uniform refinement is incompatible with the hierarchy");
3762:   for (r = 0; r < nlevels; ++r) {
3763:     CellRefiner cellRefiner;

3765:     DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);
3766:     DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);
3767:     DMPlexSetCoarseDM(dmRefined[r], cdm);
3768:     cdm  = dmRefined[r];
3769:   }
3770:   return(0);
3771: }

3775: PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened)
3776: {
3777:   DM_Plex       *mesh = (DM_Plex*) dm->data;

3781:   PetscObjectReference((PetscObject) mesh->coarseMesh);
3782:   *dmCoarsened = mesh->coarseMesh;
3783:   return(0);
3784: }

3788: PetscErrorCode DMPlexLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
3789: {
3790:   PetscInt d;

3793:   if (!dm->maxCell) {
3794:     for (d = 0; d < dim; ++d) out[d] = in[d];
3795:   } else {
3796:     for (d = 0; d < dim; ++d) {
3797:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
3798:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
3799:       } else {
3800:         out[d] = in[d];
3801:       }
3802:     }
3803:   }
3804:   return(0);
3805: }

3809: PetscErrorCode DMPlexLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
3810: {
3811:   PetscInt d;

3814:   if (!dm->maxCell) {
3815:     for (d = 0; d < dim; ++d) out[d] += in[d];
3816:   } else {
3817:     for (d = 0; d < dim; ++d) {
3818:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
3819:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
3820:       } else {
3821:         out[d] += in[d];
3822:       }
3823:     }
3824:   }
3825:   return(0);
3826: }

3830: PetscErrorCode DMPlexLocalizeCoordinates(DM dm)
3831: {
3832:   PetscSection   coordSection, cSection;
3833:   Vec            coordinates,  cVec;
3834:   PetscScalar   *coords, *coords2, *anchor;
3835:   PetscInt       Nc, cStart, cEnd, c, vStart, vEnd, v, dof, d, off, off2, bs, coordSize;

3840:   if (!dm->maxCell) return(0);
3841:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
3842:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3843:   DMGetCoordinatesLocal(dm, &coordinates);
3844:   DMGetCoordinateSection(dm, &coordSection);
3845:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
3846:   PetscSectionSetNumFields(cSection, 1);
3847:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
3848:   PetscSectionSetFieldComponents(cSection, 0, Nc);
3849:   PetscSectionSetChart(cSection, cStart, vEnd);
3850:   for (v = vStart; v < vEnd; ++v) {
3851:     PetscSectionGetDof(coordSection, v, &dof);
3852:     PetscSectionSetDof(cSection,     v,  dof);
3853:     PetscSectionSetFieldDof(cSection, v, 0, dof);
3854:   }
3855:   for (c = cStart; c < cEnd; ++c) {
3856:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, NULL);
3857:     PetscSectionSetDof(cSection, c, dof);
3858:     PetscSectionSetFieldDof(cSection, c, 0, dof);
3859:   }
3860:   PetscSectionSetUp(cSection);
3861:   PetscSectionGetStorageSize(cSection, &coordSize);
3862:   VecCreate(PetscObjectComm((PetscObject) dm), &cVec);
3863:   VecGetBlockSize(coordinates, &bs);
3864:   VecSetBlockSize(cVec,         bs);
3865:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
3866:   VecSetType(cVec,VECSTANDARD);
3867:   VecGetArray(coordinates, &coords);
3868:   VecGetArray(cVec,        &coords2);
3869:   for (v = vStart; v < vEnd; ++v) {
3870:     PetscSectionGetDof(coordSection, v, &dof);
3871:     PetscSectionGetOffset(coordSection, v, &off);
3872:     PetscSectionGetOffset(cSection,     v, &off2);
3873:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
3874:   }
3875:   DMGetWorkArray(dm, 3, PETSC_SCALAR, &anchor);
3876:   for (c = cStart; c < cEnd; ++c) {
3877:     PetscScalar *cellCoords = NULL;
3878:     PetscInt     b;

3880:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
3881:     PetscSectionGetOffset(cSection, c, &off2);
3882:     for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
3883:     for (d = 0; d < dof/bs; ++d) {DMPlexLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
3884:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
3885:   }
3886:   DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);
3887:   VecRestoreArray(coordinates, &coords);
3888:   VecRestoreArray(cVec,        &coords2);
3889:   DMSetCoordinateSection(dm, cSection);
3890:   DMSetCoordinatesLocal(dm, cVec);
3891:   VecDestroy(&cVec);
3892:   PetscSectionDestroy(&cSection);
3893:   return(0);
3894: }

3898: /*@
3899:   DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point

3901:   Not Collective

3903:   Input Parameter:
3904: . dm    - The DMPlex object

3906:   Output Parameter:
3907: . depthLabel - The DMLabel recording point depth

3909:   Level: developer

3911: .keywords: mesh, points
3912: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
3913: @*/
3914: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
3915: {
3916:   DM_Plex       *mesh = (DM_Plex*) dm->data;

3922:   if (!mesh->depthLabel) {DMPlexGetLabel(dm, "depth", &mesh->depthLabel);}
3923:   *depthLabel = mesh->depthLabel;
3924:   return(0);
3925: }

3929: /*@
3930:   DMPlexGetDepth - Get the depth of the DAG representing this mesh

3932:   Not Collective

3934:   Input Parameter:
3935: . dm    - The DMPlex object

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

3940:   Level: developer

3942: .keywords: mesh, points
3943: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
3944: @*/
3945: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
3946: {
3947:   DMLabel        label;
3948:   PetscInt       d = 0;

3954:   DMPlexGetDepthLabel(dm, &label);
3955:   if (label) {DMLabelGetNumValues(label, &d);}
3956:   *depth = d-1;
3957:   return(0);
3958: }

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

3965:   Not Collective

3967:   Input Parameters:
3968: + dm           - The DMPlex object
3969: - stratumValue - The requested depth

3971:   Output Parameters:
3972: + start - The first point at this depth
3973: - end   - One beyond the last point at this depth

3975:   Level: developer

3977: .keywords: mesh, points
3978: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
3979: @*/
3980: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
3981: {
3982:   DMLabel        label;
3983:   PetscInt       pStart, pEnd;

3990:   DMPlexGetChart(dm, &pStart, &pEnd);
3991:   if (pStart == pEnd) return(0);
3992:   if (stratumValue < 0) {
3993:     if (start) *start = pStart;
3994:     if (end)   *end   = pEnd;
3995:     return(0);
3996:   }
3997:   DMPlexGetDepthLabel(dm, &label);
3998:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
3999:   DMLabelGetStratumBounds(label, stratumValue, start, end);
4000:   return(0);
4001: }

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

4008:   Not Collective

4010:   Input Parameters:
4011: + dm           - The DMPlex object
4012: - stratumValue - The requested height

4014:   Output Parameters:
4015: + start - The first point at this height
4016: - end   - One beyond the last point at this height

4018:   Level: developer

4020: .keywords: mesh, points
4021: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
4022: @*/
4023: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
4024: {
4025:   DMLabel        label;
4026:   PetscInt       depth, pStart, pEnd;

4033:   DMPlexGetChart(dm, &pStart, &pEnd);
4034:   if (pStart == pEnd) return(0);
4035:   if (stratumValue < 0) {
4036:     if (start) *start = pStart;
4037:     if (end)   *end   = pEnd;
4038:     return(0);
4039:   }
4040:   DMPlexGetDepthLabel(dm, &label);
4041:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
4042:   DMLabelGetNumValues(label, &depth);
4043:   DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
4044:   return(0);
4045: }

4049: /* Set the number of dof on each point and separate by fields */
4050: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
4051: {
4052:   PetscInt      *numDofTot;
4053:   PetscInt       depth, pStart = 0, pEnd = 0;
4054:   PetscInt       p, d, dep, f;

4058:   PetscMalloc1((dim+1), &numDofTot);
4059:   for (d = 0; d <= dim; ++d) {
4060:     numDofTot[d] = 0;
4061:     for (f = 0; f < numFields; ++f) numDofTot[d] += numDof[f*(dim+1)+d];
4062:   }
4063:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
4064:   if (numFields > 0) {
4065:     PetscSectionSetNumFields(*section, numFields);
4066:     if (numComp) {
4067:       for (f = 0; f < numFields; ++f) {
4068:         PetscSectionSetFieldComponents(*section, f, numComp[f]);
4069:       }
4070:     }
4071:   }
4072:   DMPlexGetChart(dm, &pStart, &pEnd);
4073:   PetscSectionSetChart(*section, pStart, pEnd);
4074:   DMPlexGetDepth(dm, &depth);
4075:   for (dep = 0; dep <= depth; ++dep) {
4076:     d    = dim == depth ? dep : (!dep ? 0 : dim);
4077:     DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
4078:     for (p = pStart; p < pEnd; ++p) {
4079:       for (f = 0; f < numFields; ++f) {
4080:         PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
4081:       }
4082:       PetscSectionSetDof(*section, p, numDofTot[d]);
4083:     }
4084:   }
4085:   PetscFree(numDofTot);
4086:   return(0);
4087: }

4091: /* Set the number of dof on each point and separate by fields
4092:    If constDof is PETSC_DETERMINE, constrain every dof on the point
4093: */
4094: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC,const PetscInt bcField[],const IS bcPoints[], PetscInt constDof, PetscSection section)
4095: {
4096:   PetscInt       numFields;
4097:   PetscInt       bc;

4101:   PetscSectionGetNumFields(section, &numFields);
4102:   for (bc = 0; bc < numBC; ++bc) {
4103:     PetscInt        field = 0;
4104:     const PetscInt *idx;
4105:     PetscInt        n, i;

4107:     if (numFields) field = bcField[bc];
4108:     ISGetLocalSize(bcPoints[bc], &n);
4109:     ISGetIndices(bcPoints[bc], &idx);
4110:     for (i = 0; i < n; ++i) {
4111:       const PetscInt p        = idx[i];
4112:       PetscInt       numConst = constDof;

4114:       /* Constrain every dof on the point */
4115:       if (numConst < 0) {
4116:         if (numFields) {
4117:           PetscSectionGetFieldDof(section, p, field, &numConst);
4118:         } else {
4119:           PetscSectionGetDof(section, p, &numConst);
4120:         }
4121:       }
4122:       if (numFields) {
4123:         PetscSectionAddFieldConstraintDof(section, p, field, numConst);
4124:       }
4125:       PetscSectionAddConstraintDof(section, p, numConst);
4126:     }
4127:     ISRestoreIndices(bcPoints[bc], &idx);
4128:   }
4129:   return(0);
4130: }

4134: /* Set the constrained indices on each point and separate by fields */
4135: PetscErrorCode DMPlexCreateSectionBCIndicesAll(DM dm, PetscSection section)
4136: {
4137:   PetscInt      *maxConstraints;
4138:   PetscInt       numFields, f, pStart = 0, pEnd = 0, p;

4142:   PetscSectionGetNumFields(section, &numFields);
4143:   PetscSectionGetChart(section, &pStart, &pEnd);
4144:   PetscMalloc1((numFields+1), &maxConstraints);
4145:   for (f = 0; f <= numFields; ++f) maxConstraints[f] = 0;
4146:   for (p = pStart; p < pEnd; ++p) {
4147:     PetscInt cdof;

4149:     if (numFields) {
4150:       for (f = 0; f < numFields; ++f) {
4151:         PetscSectionGetFieldConstraintDof(section, p, f, &cdof);
4152:         maxConstraints[f] = PetscMax(maxConstraints[f], cdof);
4153:       }
4154:     } else {
4155:       PetscSectionGetConstraintDof(section, p, &cdof);
4156:       maxConstraints[0] = PetscMax(maxConstraints[0], cdof);
4157:     }
4158:   }
4159:   for (f = 0; f < numFields; ++f) {
4160:     maxConstraints[numFields] += maxConstraints[f];
4161:   }
4162:   if (maxConstraints[numFields]) {
4163:     PetscInt *indices;

4165:     PetscMalloc1(maxConstraints[numFields], &indices);
4166:     for (p = pStart; p < pEnd; ++p) {
4167:       PetscInt cdof, d;

4169:       PetscSectionGetConstraintDof(section, p, &cdof);
4170:       if (cdof) {
4171:         if (cdof > maxConstraints[numFields]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, "Likely memory corruption, point %D cDof %D > maxConstraints %D", p, cdof, maxConstraints[numFields]);
4172:         if (numFields) {
4173:           PetscInt numConst = 0, foff = 0;

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

4178:             PetscSectionGetFieldDof(section, p, f, &fdof);
4179:             PetscSectionGetFieldConstraintDof(section, p, f, &cfdof);
4180:             /* Change constraint numbering from absolute local dof number to field relative local dof number */
4181:             for (d = 0; d < cfdof; ++d) indices[numConst+d] = d;
4182:             PetscSectionSetFieldConstraintIndices(section, p, f, &indices[numConst]);
4183:             for (d = 0; d < cfdof; ++d) indices[numConst+d] += foff;
4184:             numConst += cfdof;
4185:             foff     += fdof;
4186:           }
4187:           if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
4188:         } else {
4189:           for (d = 0; d < cdof; ++d) indices[d] = d;
4190:         }
4191:         PetscSectionSetConstraintIndices(section, p, indices);
4192:       }
4193:     }
4194:     PetscFree(indices);
4195:   }
4196:   PetscFree(maxConstraints);
4197:   return(0);
4198: }

4202: /* Set the constrained field indices on each point */
4203: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt field, IS bcPoints, IS constraintIndices, PetscSection section)
4204: {
4205:   const PetscInt *points, *indices;
4206:   PetscInt        numFields, maxDof, numPoints, p, numConstraints;
4207:   PetscErrorCode  ierr;

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

4213:   ISGetLocalSize(bcPoints, &numPoints);
4214:   ISGetIndices(bcPoints, &points);
4215:   if (!constraintIndices) {
4216:     PetscInt *idx, i;

4218:     PetscSectionGetMaxDof(section, &maxDof);
4219:     PetscMalloc1(maxDof, &idx);
4220:     for (i = 0; i < maxDof; ++i) idx[i] = i;
4221:     for (p = 0; p < numPoints; ++p) {
4222:       PetscSectionSetFieldConstraintIndices(section, points[p], field, idx);
4223:     }
4224:     PetscFree(idx);
4225:   } else {
4226:     ISGetLocalSize(constraintIndices, &numConstraints);
4227:     ISGetIndices(constraintIndices, &indices);
4228:     for (p = 0; p < numPoints; ++p) {
4229:       PetscInt fcdof;

4231:       PetscSectionGetFieldConstraintDof(section, points[p], field, &fcdof);
4232:       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);
4233:       PetscSectionSetFieldConstraintIndices(section, points[p], field, indices);
4234:     }
4235:     ISRestoreIndices(constraintIndices, &indices);
4236:   }
4237:   ISRestoreIndices(bcPoints, &points);
4238:   return(0);
4239: }

4243: /* Set the constrained indices on each point and separate by fields */
4244: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
4245: {
4246:   PetscInt      *indices;
4247:   PetscInt       numFields, maxDof, f, pStart = 0, pEnd = 0, p;

4251:   PetscSectionGetMaxDof(section, &maxDof);
4252:   PetscMalloc1(maxDof, &indices);
4253:   PetscSectionGetNumFields(section, &numFields);
4254:   if (!numFields) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This function only works after users have set field constraint indices.");
4255:   PetscSectionGetChart(section, &pStart, &pEnd);
4256:   for (p = pStart; p < pEnd; ++p) {
4257:     PetscInt cdof, d;

4259:     PetscSectionGetConstraintDof(section, p, &cdof);
4260:     if (cdof) {
4261:       PetscInt numConst = 0, foff = 0;

4263:       for (f = 0; f < numFields; ++f) {
4264:         const PetscInt *fcind;
4265:         PetscInt        fdof, fcdof;

4267:         PetscSectionGetFieldDof(section, p, f, &fdof);
4268:         PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
4269:         if (fcdof) {PetscSectionGetFieldConstraintIndices(section, p, f, &fcind);}
4270:         /* Change constraint numbering from field relative local dof number to absolute local dof number */
4271:         for (d = 0; d < fcdof; ++d) indices[numConst+d] = fcind[d]+foff;
4272:         foff     += fdof;
4273:         numConst += fcdof;
4274:       }
4275:       if (cdof != numConst) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
4276:       PetscSectionSetConstraintIndices(section, p, indices);
4277:     }
4278:   }
4279:   PetscFree(indices);
4280:   return(0);
4281: }

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

4288:   Not Collective

4290:   Input Parameters:
4291: + dm        - The DMPlex object
4292: . dim       - The spatial dimension of the problem
4293: . numFields - The number of fields in the problem
4294: . numComp   - An array of size numFields that holds the number of components for each field
4295: . numDof    - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
4296: . numBC     - The number of boundary conditions
4297: . bcField   - An array of size numBC giving the field number for each boundry condition
4298: . bcPoints  - An array of size numBC giving an IS holding the sieve points to which each boundary condition applies
4299: - perm      - Optional permutation of the chart, or NULL

4301:   Output Parameter:
4302: . section - The PetscSection object

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

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

4309:   Level: developer

4311:   Fortran Notes:
4312:   A Fortran 90 version is available as DMPlexCreateSectionF90()

4314: .keywords: mesh, elements
4315: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
4316: @*/
4317: 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)
4318: {

4322:   DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
4323:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcPoints, PETSC_DETERMINE, *section);
4324:   if (perm) {PetscSectionSetPermutation(*section, perm);}
4325:   PetscSectionSetUp(*section);
4326:   if (numBC) {DMPlexCreateSectionBCIndicesAll(dm, *section);}
4327:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
4328:   return(0);
4329: }

4333: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
4334: {
4335:   PetscSection   section;

4339:   DMClone(dm, cdm);
4340:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
4341:   DMSetDefaultSection(*cdm, section);
4342:   PetscSectionDestroy(&section);
4343:   return(0);
4344: }

4348: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
4349: {
4350:   DM_Plex *mesh = (DM_Plex*) dm->data;

4354:   if (section) *section = mesh->coneSection;
4355:   return(0);
4356: }

4360: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
4361: {
4362:   DM_Plex *mesh = (DM_Plex*) dm->data;

4366:   if (section) *section = mesh->supportSection;
4367:   return(0);
4368: }

4372: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
4373: {
4374:   DM_Plex *mesh = (DM_Plex*) dm->data;

4378:   if (cones) *cones = mesh->cones;
4379:   return(0);
4380: }

4384: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
4385: {
4386:   DM_Plex *mesh = (DM_Plex*) dm->data;

4390:   if (coneOrientations) *coneOrientations = mesh->coneOrientations;
4391:   return(0);
4392: }

4394: /******************************** FEM Support **********************************/

4398: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
4399: {
4400:   PetscScalar    *array, *vArray;
4401:   const PetscInt *cone, *coneO;
4402:   PetscInt        pStart, pEnd, p, numPoints, size = 0, offset = 0;
4403:   PetscErrorCode  ierr;

4406:   PetscSectionGetChart(section, &pStart, &pEnd);
4407:   DMPlexGetConeSize(dm, point, &numPoints);
4408:   DMPlexGetCone(dm, point, &cone);
4409:   DMPlexGetConeOrientation(dm, point, &coneO);
4410:   if (!values || !*values) {
4411:     if ((point >= pStart) && (point < pEnd)) {
4412:       PetscInt dof;

4414:       PetscSectionGetDof(section, point, &dof);
4415:       size += dof;
4416:     }
4417:     for (p = 0; p < numPoints; ++p) {
4418:       const PetscInt cp = cone[p];
4419:       PetscInt       dof;

4421:       if ((cp < pStart) || (cp >= pEnd)) continue;
4422:       PetscSectionGetDof(section, cp, &dof);
4423:       size += dof;
4424:     }
4425:     if (!values) {
4426:       if (csize) *csize = size;
4427:       return(0);
4428:     }
4429:     DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
4430:   } else {
4431:     array = *values;
4432:   }
4433:   size = 0;
4434:   VecGetArray(v, &vArray);
4435:   if ((point >= pStart) && (point < pEnd)) {
4436:     PetscInt     dof, off, d;
4437:     PetscScalar *varr;

4439:     PetscSectionGetDof(section, point, &dof);
4440:     PetscSectionGetOffset(section, point, &off);
4441:     varr = &vArray[off];
4442:     for (d = 0; d < dof; ++d, ++offset) {
4443:       array[offset] = varr[d];
4444:     }
4445:     size += dof;
4446:   }
4447:   for (p = 0; p < numPoints; ++p) {
4448:     const PetscInt cp = cone[p];
4449:     PetscInt       o  = coneO[p];
4450:     PetscInt       dof, off, d;
4451:     PetscScalar   *varr;

4453:     if ((cp < pStart) || (cp >= pEnd)) continue;
4454:     PetscSectionGetDof(section, cp, &dof);
4455:     PetscSectionGetOffset(section, cp, &off);
4456:     varr = &vArray[off];
4457:     if (o >= 0) {
4458:       for (d = 0; d < dof; ++d, ++offset) {
4459:         array[offset] = varr[d];
4460:       }
4461:     } else {
4462:       for (d = dof-1; d >= 0; --d, ++offset) {
4463:         array[offset] = varr[d];
4464:       }
4465:     }
4466:     size += dof;
4467:   }
4468:   VecRestoreArray(v, &vArray);
4469:   if (!*values) {
4470:     if (csize) *csize = size;
4471:     *values = array;
4472:   } else {
4473:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
4474:     *csize = size;
4475:   }
4476:   return(0);
4477: }

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

4487:   *size = 0;
4488:   for (p = 0; p < numPoints*2; p += 2) {
4489:     const PetscInt point = points[p];
4490:     const PetscInt o     = points[p+1];
4491:     PetscInt       dof, off, d;
4492:     const PetscScalar *varr;

4494:     PetscSectionGetDof(section, point, &dof);
4495:     PetscSectionGetOffset(section, point, &off);
4496:     varr = &vArray[off];
4497:     if (o >= 0) {
4498:       for (d = 0; d < dof; ++d, ++offset)    array[offset] = varr[d];
4499:     } else {
4500:       for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
4501:     }
4502:   }
4503:   *size = offset;
4504:   return(0);
4505: }

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

4515:   *size = 0;
4516:   for (f = 0; f < numFields; ++f) {
4517:     PetscInt fcomp, p;

4519:     PetscSectionGetFieldComponents(section, f, &fcomp);
4520:     for (p = 0; p < numPoints*2; p += 2) {
4521:       const PetscInt point = points[p];
4522:       const PetscInt o     = points[p+1];
4523:       PetscInt       fdof, foff, d, c;
4524:       const PetscScalar *varr;

4526:       PetscSectionGetFieldDof(section, point, f, &fdof);
4527:       PetscSectionGetFieldOffset(section, point, f, &foff);
4528:       varr = &vArray[foff];
4529:       if (o >= 0) {
4530:         for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
4531:       } else {
4532:         for (d = fdof/fcomp-1; d >= 0; --d) {
4533:           for (c = 0; c < fcomp; ++c, ++offset) {
4534:             array[offset] = varr[d*fcomp+c];
4535:           }
4536:         }
4537:       }
4538:     }
4539:   }
4540:   *size = offset;
4541:   return(0);
4542: }

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

4549:   Not collective

4551:   Input Parameters:
4552: + dm - The DM
4553: . section - The section describing the layout in v, or NULL to use the default section
4554: . v - The local vector
4555: - point - The sieve point in the DM

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

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

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

4567:   Level: intermediate

4569: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
4570: @*/
4571: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
4572: {
4573:   PetscSection    clSection;
4574:   IS              clPoints;
4575:   PetscScalar    *array, *vArray;
4576:   PetscInt       *points = NULL;
4577:   const PetscInt *clp;
4578:   PetscInt        depth, numFields, numPoints, size;
4579:   PetscErrorCode  ierr;

4583:   if (!section) {DMGetDefaultSection(dm, &section);}
4586:   DMPlexGetDepth(dm, &depth);
4587:   PetscSectionGetNumFields(section, &numFields);
4588:   if (depth == 1 && numFields < 2) {
4589:     DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
4590:     return(0);
4591:   }
4592:   /* Get points */
4593:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4594:   if (!clPoints) {
4595:     PetscInt pStart, pEnd, p, q;

4597:     PetscSectionGetChart(section, &pStart, &pEnd);
4598:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4599:     /* Compress out points not in the section */
4600:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4601:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4602:         points[q*2]   = points[p];
4603:         points[q*2+1] = points[p+1];
4604:         ++q;
4605:       }
4606:     }
4607:     numPoints = q;
4608:   } else {
4609:     PetscInt dof, off;

4611:     PetscSectionGetDof(clSection, point, &dof);
4612:     PetscSectionGetOffset(clSection, point, &off);
4613:     ISGetIndices(clPoints, &clp);
4614:     numPoints = dof/2;
4615:     points    = (PetscInt *) &clp[off];
4616:   }
4617:   /* Get array */
4618:   if (!values || !*values) {
4619:     PetscInt asize = 0, dof, p;

4621:     for (p = 0; p < numPoints*2; p += 2) {
4622:       PetscSectionGetDof(section, points[p], &dof);
4623:       asize += dof;
4624:     }
4625:     if (!values) {
4626:       if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4627:       else           {ISRestoreIndices(clPoints, &clp);}
4628:       if (csize) *csize = asize;
4629:       return(0);
4630:     }
4631:     DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
4632:   } else {
4633:     array = *values;
4634:   }
4635:   VecGetArray(v, &vArray);
4636:   /* Get values */
4637:   if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);}
4638:   else               {DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);}
4639:   /* Cleanup points */
4640:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4641:   else           {ISRestoreIndices(clPoints, &clp);}
4642:   /* Cleanup array */
4643:   VecRestoreArray(v, &vArray);
4644:   if (!*values) {
4645:     if (csize) *csize = size;
4646:     *values = array;
4647:   } else {
4648:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
4649:     *csize = size;
4650:   }
4651:   return(0);
4652: }

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

4659:   Not collective

4661:   Input Parameters:
4662: + dm - The DM
4663: . section - The section describing the layout in v, or NULL to use the default section
4664: . v - The local vector
4665: . point - The sieve point in the DM
4666: . csize - The number of values in the closure, or NULL
4667: - values - The array of values, which is a borrowed array and should not be freed

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

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

4675:   Level: intermediate

4677: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
4678: @*/
4679: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
4680: {
4681:   PetscInt       size = 0;

4685:   /* Should work without recalculating size */
4686:   DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
4687:   return(0);
4688: }

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

4695: 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[])
4696: {
4697:   PetscInt        cdof;   /* The number of constraints on this point */
4698:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
4699:   PetscScalar    *a;
4700:   PetscInt        off, cind = 0, k;
4701:   PetscErrorCode  ierr;

4704:   PetscSectionGetConstraintDof(section, point, &cdof);
4705:   PetscSectionGetOffset(section, point, &off);
4706:   a    = &array[off];
4707:   if (!cdof || setBC) {
4708:     if (orientation >= 0) {
4709:       for (k = 0; k < dof; ++k) {
4710:         fuse(&a[k], values[k]);
4711:       }
4712:     } else {
4713:       for (k = 0; k < dof; ++k) {
4714:         fuse(&a[k], values[dof-k-1]);
4715:       }
4716:     }
4717:   } else {
4718:     PetscSectionGetConstraintIndices(section, point, &cdofs);
4719:     if (orientation >= 0) {
4720:       for (k = 0; k < dof; ++k) {
4721:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4722:         fuse(&a[k], values[k]);
4723:       }
4724:     } else {
4725:       for (k = 0; k < dof; ++k) {
4726:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4727:         fuse(&a[k], values[dof-k-1]);
4728:       }
4729:     }
4730:   }
4731:   return(0);
4732: }

4736: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
4737: {
4738:   PetscInt        cdof;   /* The number of constraints on this point */
4739:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
4740:   PetscScalar    *a;
4741:   PetscInt        off, cind = 0, k;
4742:   PetscErrorCode  ierr;

4745:   PetscSectionGetConstraintDof(section, point, &cdof);
4746:   PetscSectionGetOffset(section, point, &off);
4747:   a    = &array[off];
4748:   if (cdof) {
4749:     PetscSectionGetConstraintIndices(section, point, &cdofs);
4750:     if (orientation >= 0) {
4751:       for (k = 0; k < dof; ++k) {
4752:         if ((cind < cdof) && (k == cdofs[cind])) {
4753:           fuse(&a[k], values[k]);
4754:           ++cind;
4755:         }
4756:       }
4757:     } else {
4758:       for (k = 0; k < dof; ++k) {
4759:         if ((cind < cdof) && (k == cdofs[cind])) {
4760:           fuse(&a[k], values[dof-k-1]);
4761:           ++cind;
4762:         }
4763:       }
4764:     }
4765:   }
4766:   return(0);
4767: }

4771: 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[])
4772: {
4773:   PetscScalar    *a;
4774:   PetscInt        fdof, foff, fcdof, foffset = *offset;
4775:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4776:   PetscInt        cind = 0, k, c;
4777:   PetscErrorCode  ierr;

4780:   PetscSectionGetFieldDof(section, point, f, &fdof);
4781:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
4782:   PetscSectionGetFieldOffset(section, point, f, &foff);
4783:   a    = &array[foff];
4784:   if (!fcdof || setBC) {
4785:     if (o >= 0) {
4786:       for (k = 0; k < fdof; ++k) fuse(&a[k], values[foffset+k]);
4787:     } else {
4788:       for (k = fdof/fcomp-1; k >= 0; --k) {
4789:         for (c = 0; c < fcomp; ++c) {
4790:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
4791:         }
4792:       }
4793:     }
4794:   } else {
4795:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
4796:     if (o >= 0) {
4797:       for (k = 0; k < fdof; ++k) {
4798:         if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
4799:         fuse(&a[k], values[foffset+k]);
4800:       }
4801:     } else {
4802:       for (k = fdof/fcomp-1; k >= 0; --k) {
4803:         for (c = 0; c < fcomp; ++c) {
4804:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
4805:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
4806:         }
4807:       }
4808:     }
4809:   }
4810:   *offset += fdof;
4811:   return(0);
4812: }

4816: 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[])
4817: {
4818:   PetscScalar    *a;
4819:   PetscInt        fdof, foff, fcdof, foffset = *offset;
4820:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4821:   PetscInt        cind = 0, k, c;
4822:   PetscErrorCode  ierr;

4825:   PetscSectionGetFieldDof(section, point, f, &fdof);
4826:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
4827:   PetscSectionGetFieldOffset(section, point, f, &foff);
4828:   a    = &array[foff];
4829:   if (fcdof) {
4830:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
4831:     if (o >= 0) {
4832:       for (k = 0; k < fdof; ++k) {
4833:         if ((cind < fcdof) && (k == fcdofs[cind])) {
4834:           fuse(&a[k], values[foffset+k]);
4835:           ++cind;
4836:         }
4837:       }
4838:     } else {
4839:       for (k = fdof/fcomp-1; k >= 0; --k) {
4840:         for (c = 0; c < fcomp; ++c) {
4841:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {
4842:             fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
4843:             ++cind;
4844:           }
4845:         }
4846:       }
4847:     }
4848:   }
4849:   *offset += fdof;
4850:   return(0);
4851: }

4855: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
4856: {
4857:   PetscScalar    *array;
4858:   const PetscInt *cone, *coneO;
4859:   PetscInt        pStart, pEnd, p, numPoints, off, dof;
4860:   PetscErrorCode  ierr;

4863:   PetscSectionGetChart(section, &pStart, &pEnd);
4864:   DMPlexGetConeSize(dm, point, &numPoints);
4865:   DMPlexGetCone(dm, point, &cone);
4866:   DMPlexGetConeOrientation(dm, point, &coneO);
4867:   VecGetArray(v, &array);
4868:   for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
4869:     const PetscInt cp = !p ? point : cone[p-1];
4870:     const PetscInt o  = !p ? 0     : coneO[p-1];

4872:     if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
4873:     PetscSectionGetDof(section, cp, &dof);
4874:     /* ADD_VALUES */
4875:     {
4876:       const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
4877:       PetscScalar    *a;
4878:       PetscInt        cdof, coff, cind = 0, k;

4880:       PetscSectionGetConstraintDof(section, cp, &cdof);
4881:       PetscSectionGetOffset(section, cp, &coff);
4882:       a    = &array[coff];
4883:       if (!cdof) {
4884:         if (o >= 0) {
4885:           for (k = 0; k < dof; ++k) {
4886:             a[k] += values[off+k];
4887:           }
4888:         } else {
4889:           for (k = 0; k < dof; ++k) {
4890:             a[k] += values[off+dof-k-1];
4891:           }
4892:         }
4893:       } else {
4894:         PetscSectionGetConstraintIndices(section, cp, &cdofs);
4895:         if (o >= 0) {
4896:           for (k = 0; k < dof; ++k) {
4897:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4898:             a[k] += values[off+k];
4899:           }
4900:         } else {
4901:           for (k = 0; k < dof; ++k) {
4902:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4903:             a[k] += values[off+dof-k-1];
4904:           }
4905:         }
4906:       }
4907:     }
4908:   }
4909:   VecRestoreArray(v, &array);
4910:   return(0);
4911: }

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

4918:   Not collective

4920:   Input Parameters:
4921: + dm - The DM
4922: . section - The section describing the layout in v, or NULL to use the default section
4923: . v - The local vector
4924: . point - The sieve point in the DM
4925: . values - The array of values
4926: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

4931:   Level: intermediate

4933: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
4934: @*/
4935: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
4936: {
4937:   PetscSection    clSection;
4938:   IS              clPoints;
4939:   PetscScalar    *array;
4940:   PetscInt       *points = NULL;
4941:   const PetscInt *clp;
4942:   PetscInt        depth, numFields, numPoints, p;
4943:   PetscErrorCode  ierr;

4947:   if (!section) {DMGetDefaultSection(dm, &section);}
4950:   DMPlexGetDepth(dm, &depth);
4951:   PetscSectionGetNumFields(section, &numFields);
4952:   if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
4953:     DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
4954:     return(0);
4955:   }
4956:   /* Get points */
4957:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4958:   if (!clPoints) {
4959:     PetscInt pStart, pEnd, q;

4961:     PetscSectionGetChart(section, &pStart, &pEnd);
4962:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4963:     /* Compress out points not in the section */
4964:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4965:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4966:         points[q*2]   = points[p];
4967:         points[q*2+1] = points[p+1];
4968:         ++q;
4969:       }
4970:     }
4971:     numPoints = q;
4972:   } else {
4973:     PetscInt dof, off;

4975:     PetscSectionGetDof(clSection, point, &dof);
4976:     PetscSectionGetOffset(clSection, point, &off);
4977:     ISGetIndices(clPoints, &clp);
4978:     numPoints = dof/2;
4979:     points    = (PetscInt *) &clp[off];
4980:   }
4981:   /* Get array */
4982:   VecGetArray(v, &array);
4983:   /* Get values */
4984:   if (numFields > 0) {
4985:     PetscInt offset = 0, fcomp, f;
4986:     for (f = 0; f < numFields; ++f) {
4987:       PetscSectionGetFieldComponents(section, f, &fcomp);
4988:       switch (mode) {
4989:       case INSERT_VALUES:
4990:         for (p = 0; p < numPoints*2; p += 2) {
4991:           const PetscInt point = points[p];
4992:           const PetscInt o     = points[p+1];
4993:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
4994:         } break;
4995:       case INSERT_ALL_VALUES:
4996:         for (p = 0; p < numPoints*2; p += 2) {
4997:           const PetscInt point = points[p];
4998:           const PetscInt o     = points[p+1];
4999:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
5000:         } break;
5001:       case INSERT_BC_VALUES:
5002:         for (p = 0; p < numPoints*2; p += 2) {
5003:           const PetscInt point = points[p];
5004:           const PetscInt o     = points[p+1];
5005:           updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
5006:         } break;
5007:       case ADD_VALUES:
5008:         for (p = 0; p < numPoints*2; p += 2) {
5009:           const PetscInt point = points[p];
5010:           const PetscInt o     = points[p+1];
5011:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
5012:         } break;
5013:       case ADD_ALL_VALUES:
5014:         for (p = 0; p < numPoints*2; p += 2) {
5015:           const PetscInt point = points[p];
5016:           const PetscInt o     = points[p+1];
5017:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
5018:         } break;
5019:       default:
5020:         SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
5021:       }
5022:     }
5023:   } else {
5024:     PetscInt dof, off;

5026:     switch (mode) {
5027:     case INSERT_VALUES:
5028:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5029:         PetscInt o = points[p+1];
5030:         PetscSectionGetDof(section, points[p], &dof);
5031:         updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
5032:       } break;
5033:     case INSERT_ALL_VALUES:
5034:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5035:         PetscInt o = points[p+1];
5036:         PetscSectionGetDof(section, points[p], &dof);
5037:         updatePoint_private(section, points[p], dof, insert, PETSC_TRUE,  o, &values[off], array);
5038:       } break;
5039:     case INSERT_BC_VALUES:
5040:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5041:         PetscInt o = points[p+1];
5042:         PetscSectionGetDof(section, points[p], &dof);
5043:         updatePointBC_private(section, points[p], dof, insert,  o, &values[off], array);
5044:       } break;
5045:     case ADD_VALUES:
5046:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5047:         PetscInt o = points[p+1];
5048:         PetscSectionGetDof(section, points[p], &dof);
5049:         updatePoint_private(section, points[p], dof, add,    PETSC_FALSE, o, &values[off], array);
5050:       } break;
5051:     case ADD_ALL_VALUES:
5052:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5053:         PetscInt o = points[p+1];
5054:         PetscSectionGetDof(section, points[p], &dof);
5055:         updatePoint_private(section, points[p], dof, add,    PETSC_TRUE,  o, &values[off], array);
5056:       } break;
5057:     default:
5058:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
5059:     }
5060:   }
5061:   /* Cleanup points */
5062:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
5063:   else           {ISRestoreIndices(clPoints, &clp);}
5064:   /* Cleanup array */
5065:   VecRestoreArray(v, &array);
5066:   return(0);
5067: }

5071: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
5072: {
5073:   PetscSection    clSection;
5074:   IS              clPoints;
5075:   PetscScalar    *array;
5076:   PetscInt       *points = NULL;
5077:   const PetscInt *clp;
5078:   PetscInt        numFields, numPoints, p;
5079:   PetscInt        offset = 0, fcomp, f;
5080:   PetscErrorCode  ierr;

5084:   if (!section) {DMGetDefaultSection(dm, &section);}
5087:   PetscSectionGetNumFields(section, &numFields);
5088:   /* Get points */
5089:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
5090:   if (!clPoints) {
5091:     PetscInt pStart, pEnd, q;

5093:     PetscSectionGetChart(section, &pStart, &pEnd);
5094:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
5095:     /* Compress out points not in the section */
5096:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
5097:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
5098:         points[q*2]   = points[p];
5099:         points[q*2+1] = points[p+1];
5100:         ++q;
5101:       }
5102:     }
5103:     numPoints = q;
5104:   } else {
5105:     PetscInt dof, off;

5107:     PetscSectionGetDof(clSection, point, &dof);
5108:     PetscSectionGetOffset(clSection, point, &off);
5109:     ISGetIndices(clPoints, &clp);
5110:     numPoints = dof/2;
5111:     points    = (PetscInt *) &clp[off];
5112:   }
5113:   /* Get array */
5114:   VecGetArray(v, &array);
5115:   /* Get values */
5116:   for (f = 0; f < numFields; ++f) {
5117:     PetscSectionGetFieldComponents(section, f, &fcomp);
5118:     if (!fieldActive[f]) {
5119:       for (p = 0; p < numPoints*2; p += 2) {
5120:         PetscInt fdof;
5121:         PetscSectionGetFieldDof(section, points[p], f, &fdof);
5122:         offset += fdof;
5123:       }
5124:       continue;
5125:     }
5126:     switch (mode) {
5127:     case INSERT_VALUES:
5128:       for (p = 0; p < numPoints*2; p += 2) {
5129:         const PetscInt point = points[p];
5130:         const PetscInt o     = points[p+1];
5131:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
5132:       } break;
5133:     case INSERT_ALL_VALUES:
5134:       for (p = 0; p < numPoints*2; p += 2) {
5135:         const PetscInt point = points[p];
5136:         const PetscInt o     = points[p+1];
5137:         updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
5138:         } break;
5139:     case INSERT_BC_VALUES:
5140:       for (p = 0; p < numPoints*2; p += 2) {
5141:         const PetscInt point = points[p];
5142:         const PetscInt o     = points[p+1];
5143:         updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
5144:       } break;
5145:     case ADD_VALUES:
5146:       for (p = 0; p < numPoints*2; p += 2) {
5147:         const PetscInt point = points[p];
5148:         const PetscInt o     = points[p+1];
5149:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
5150:       } break;
5151:     case ADD_ALL_VALUES:
5152:       for (p = 0; p < numPoints*2; p += 2) {
5153:         const PetscInt point = points[p];
5154:         const PetscInt o     = points[p+1];
5155:         updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
5156:       } break;
5157:     default:
5158:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
5159:     }
5160:   }
5161:   /* Cleanup points */
5162:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
5163:   else           {ISRestoreIndices(clPoints, &clp);}
5164:   /* Cleanup array */
5165:   VecRestoreArray(v, &array);
5166:   return(0);
5167: }

5171: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
5172: {
5173:   PetscMPIInt    rank;
5174:   PetscInt       i, j;

5178:   MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
5179:   PetscViewerASCIIPrintf(viewer, "[%D]mat for sieve point %D\n", rank, point);
5180:   for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%D]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
5181:   for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%D]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
5182:   numCIndices = numCIndices ? numCIndices : numRIndices;
5183:   for (i = 0; i < numRIndices; i++) {
5184:     PetscViewerASCIIPrintf(viewer, "[%D]", rank);
5185:     for (j = 0; j < numCIndices; j++) {
5186: #if defined(PETSC_USE_COMPLEX)
5187:       PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
5188: #else
5189:       PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
5190: #endif
5191:     }
5192:     PetscViewerASCIIPrintf(viewer, "\n");
5193:   }
5194:   return(0);
5195: }

5199: /* . off - The global offset of this point */
5200: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
5201: {
5202:   PetscInt        dof;    /* The number of unknowns on this point */
5203:   PetscInt        cdof;   /* The number of constraints on this point */
5204:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
5205:   PetscInt        cind = 0, k;
5206:   PetscErrorCode  ierr;

5209:   PetscSectionGetDof(section, point, &dof);
5210:   PetscSectionGetConstraintDof(section, point, &cdof);
5211:   if (!cdof || setBC) {
5212:     if (orientation >= 0) {
5213:       for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
5214:     } else {
5215:       for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
5216:     }
5217:   } else {
5218:     PetscSectionGetConstraintIndices(section, point, &cdofs);
5219:     if (orientation >= 0) {
5220:       for (k = 0; k < dof; ++k) {
5221:         if ((cind < cdof) && (k == cdofs[cind])) {
5222:           /* Insert check for returning constrained indices */
5223:           indices[*loff+k] = -(off+k+1);
5224:           ++cind;
5225:         } else {
5226:           indices[*loff+k] = off+k-cind;
5227:         }
5228:       }
5229:     } else {
5230:       for (k = 0; k < dof; ++k) {
5231:         if ((cind < cdof) && (k == cdofs[cind])) {
5232:           /* Insert check for returning constrained indices */
5233:           indices[*loff+dof-k-1] = -(off+k+1);
5234:           ++cind;
5235:         } else {
5236:           indices[*loff+dof-k-1] = off+k-cind;
5237:         }
5238:       }
5239:     }
5240:   }
5241:   *loff += dof;
5242:   return(0);
5243: }

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

5254:   PetscSectionGetNumFields(section, &numFields);
5255:   for (f = 0, foff = 0; f < numFields; ++f) {
5256:     PetscInt        fdof, fcomp, cfdof;
5257:     const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
5258:     PetscInt        cind = 0, k, c;

5260:     PetscSectionGetFieldComponents(section, f, &fcomp);
5261:     PetscSectionGetFieldDof(section, point, f, &fdof);
5262:     PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
5263:     if (!cfdof || setBC) {
5264:       if (orientation >= 0) {
5265:         for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
5266:       } else {
5267:         for (k = fdof/fcomp-1; k >= 0; --k) {
5268:           for (c = 0; c < fcomp; ++c) {
5269:             indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
5270:           }
5271:         }
5272:       }
5273:     } else {
5274:       PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
5275:       if (orientation >= 0) {
5276:         for (k = 0; k < fdof; ++k) {
5277:           if ((cind < cfdof) && (k == fcdofs[cind])) {
5278:             indices[foffs[f]+k] = -(off+foff+k+1);
5279:             ++cind;
5280:           } else {
5281:             indices[foffs[f]+k] = off+foff+k-cind;
5282:           }
5283:         }
5284:       } else {
5285:         for (k = fdof/fcomp-1; k >= 0; --k) {
5286:           for (c = 0; c < fcomp; ++c) {
5287:             if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
5288:               indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
5289:               ++cind;
5290:             } else {
5291:               indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
5292:             }
5293:           }
5294:         }
5295:       }
5296:     }
5297:     foff     += fdof - cfdof;
5298:     foffs[f] += fdof;
5299:   }
5300:   return(0);
5301: }

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

5308:   Not collective

5310:   Input Parameters:
5311: + dm - The DM
5312: . section - The section describing the layout in v, or NULL to use the default section
5313: . globalSection - The section describing the layout in v, or NULL to use the default global section
5314: . A - The matrix
5315: . point - The sieve point in the DM
5316: . values - The array of values
5317: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

5322:   Level: intermediate

5324: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
5325: @*/
5326: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
5327: {
5328:   DM_Plex        *mesh   = (DM_Plex*) dm->data;
5329:   PetscSection    clSection;
5330:   IS              clPoints;
5331:   PetscInt       *points = NULL;
5332:   const PetscInt *clp;
5333:   PetscInt       *indices;
5334:   PetscInt        offsets[32];
5335:   PetscInt        numFields, numPoints, numIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
5336:   PetscErrorCode  ierr;

5340:   if (!section) {DMGetDefaultSection(dm, &section);}
5342:   if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
5345:   PetscSectionGetNumFields(section, &numFields);
5346:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
5347:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
5348:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
5349:   if (!clPoints) {
5350:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
5351:     /* Compress out points not in the section */
5352:     PetscSectionGetChart(section, &pStart, &pEnd);
5353:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
5354:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
5355:         points[q*2]   = points[p];
5356:         points[q*2+1] = points[p+1];
5357:         ++q;
5358:       }
5359:     }
5360:     numPoints = q;
5361:   } else {
5362:     PetscInt dof, off;

5364:     PetscSectionGetDof(clSection, point, &dof);
5365:     numPoints = dof/2;
5366:     PetscSectionGetOffset(clSection, point, &off);
5367:     ISGetIndices(clPoints, &clp);
5368:     points = (PetscInt *) &clp[off];
5369:   }
5370:   for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
5371:     PetscInt fdof;

5373:     PetscSectionGetDof(section, points[p], &dof);
5374:     for (f = 0; f < numFields; ++f) {
5375:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
5376:       offsets[f+1] += fdof;
5377:     }
5378:     numIndices += dof;
5379:   }
5380:   for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];

5382:   if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
5383:   DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
5384:   if (numFields) {
5385:     for (p = 0; p < numPoints*2; p += 2) {
5386:       PetscInt o = points[p+1];
5387:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
5388:       indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
5389:     }
5390:   } else {
5391:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
5392:       PetscInt o = points[p+1];
5393:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
5394:       indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
5395:     }
5396:   }
5397:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
5398:   MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
5399:   if (ierr) {
5400:     PetscMPIInt    rank;
5401:     PetscErrorCode ierr2;

5403:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
5404:     ierr2 = (*PetscErrorPrintf)("[%D]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
5405:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
5406:     ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
5407: 
5408:   }
5409:   if (!clPoints) {
5410:     DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
5411:   } else {
5412:     ISRestoreIndices(clPoints, &clp);
5413:   }
5414:   DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
5415:   return(0);
5416: }

5420: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
5421: {
5422:   DM_Plex        *mesh   = (DM_Plex*) dmf->data;
5423:   PetscInt       *fpoints = NULL, *ftotpoints = NULL;
5424:   PetscInt       *cpoints = NULL;
5425:   PetscInt       *findices, *cindices;
5426:   PetscInt        foffsets[32], coffsets[32];
5427:   CellRefiner     cellRefiner;
5428:   PetscInt        numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
5429:   PetscErrorCode  ierr;

5434:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
5436:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
5438:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
5440:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
5443:   PetscSectionGetNumFields(fsection, &numFields);
5444:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
5445:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
5446:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
5447:   /* Column indices */
5448:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5449:   maxFPoints = numCPoints;
5450:   /* Compress out points not in the section */
5451:   /*   TODO: Squeeze out points with 0 dof as well */
5452:   PetscSectionGetChart(csection, &pStart, &pEnd);
5453:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
5454:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
5455:       cpoints[q*2]   = cpoints[p];
5456:       cpoints[q*2+1] = cpoints[p+1];
5457:       ++q;
5458:     }
5459:   }
5460:   numCPoints = q;
5461:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
5462:     PetscInt fdof;

5464:     PetscSectionGetDof(csection, cpoints[p], &dof);
5465:     if (!dof) continue;
5466:     for (f = 0; f < numFields; ++f) {
5467:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
5468:       coffsets[f+1] += fdof;
5469:     }
5470:     numCIndices += dof;
5471:   }
5472:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
5473:   /* Row indices */
5474:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
5475:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
5476:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
5477:   for (r = 0, q = 0; r < numSubcells; ++r) {
5478:     /* TODO Map from coarse to fine cells */
5479:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
5480:     /* Compress out points not in the section */
5481:     PetscSectionGetChart(fsection, &pStart, &pEnd);
5482:     for (p = 0; p < numFPoints*2; p += 2) {
5483:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
5484:         PetscSectionGetDof(fsection, fpoints[p], &dof);
5485:         if (!dof) continue;
5486:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
5487:         if (s < q) continue;
5488:         ftotpoints[q*2]   = fpoints[p];
5489:         ftotpoints[q*2+1] = fpoints[p+1];
5490:         ++q;
5491:       }
5492:     }
5493:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
5494:   }
5495:   numFPoints = q;
5496:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
5497:     PetscInt fdof;

5499:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
5500:     if (!dof) continue;
5501:     for (f = 0; f < numFields; ++f) {
5502:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
5503:       foffsets[f+1] += fdof;
5504:     }
5505:     numFIndices += dof;
5506:   }
5507:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

5509:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
5510:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
5511:   DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
5512:   DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
5513:   if (numFields) {
5514:     for (p = 0; p < numFPoints*2; p += 2) {
5515:       PetscInt o = ftotpoints[p+1];
5516:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5517:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
5518:     }
5519:     for (p = 0; p < numCPoints*2; p += 2) {
5520:       PetscInt o = cpoints[p+1];
5521:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5522:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5523:     }
5524:   } else {
5525:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5526:       PetscInt o = ftotpoints[p+1];
5527:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5528:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5529:     }
5530:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5531:       PetscInt o = cpoints[p+1];
5532:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5533:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5534:     }
5535:   }
5536:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
5537:   MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
5538:   if (ierr) {
5539:     PetscMPIInt    rank;
5540:     PetscErrorCode ierr2;

5542:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
5543:     ierr2 = (*PetscErrorPrintf)("[%D]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
5544:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
5545:     ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
5546:     ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
5547: 
5548:   }
5549:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5550:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5551:   DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
5552:   DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
5553:   return(0);
5554: }

5558: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
5559: {
5560:   PetscInt      *fpoints = NULL, *ftotpoints = NULL;
5561:   PetscInt      *cpoints = NULL;
5562:   PetscInt       foffsets[32], coffsets[32];
5563:   CellRefiner    cellRefiner;
5564:   PetscInt       numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;

5570:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
5572:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
5574:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
5576:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
5578:   PetscSectionGetNumFields(fsection, &numFields);
5579:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
5580:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
5581:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
5582:   /* Column indices */
5583:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5584:   maxFPoints = numCPoints;
5585:   /* Compress out points not in the section */
5586:   /*   TODO: Squeeze out points with 0 dof as well */
5587:   PetscSectionGetChart(csection, &pStart, &pEnd);
5588:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
5589:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
5590:       cpoints[q*2]   = cpoints[p];
5591:       cpoints[q*2+1] = cpoints[p+1];
5592:       ++q;
5593:     }
5594:   }
5595:   numCPoints = q;
5596:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
5597:     PetscInt fdof;

5599:     PetscSectionGetDof(csection, cpoints[p], &dof);
5600:     if (!dof) continue;
5601:     for (f = 0; f < numFields; ++f) {
5602:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
5603:       coffsets[f+1] += fdof;
5604:     }
5605:     numCIndices += dof;
5606:   }
5607:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
5608:   /* Row indices */
5609:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
5610:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
5611:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
5612:   for (r = 0, q = 0; r < numSubcells; ++r) {
5613:     /* TODO Map from coarse to fine cells */
5614:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
5615:     /* Compress out points not in the section */
5616:     PetscSectionGetChart(fsection, &pStart, &pEnd);
5617:     for (p = 0; p < numFPoints*2; p += 2) {
5618:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
5619:         PetscSectionGetDof(fsection, fpoints[p], &dof);
5620:         if (!dof) continue;
5621:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
5622:         if (s < q) continue;
5623:         ftotpoints[q*2]   = fpoints[p];
5624:         ftotpoints[q*2+1] = fpoints[p+1];
5625:         ++q;
5626:       }
5627:     }
5628:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
5629:   }
5630:   numFPoints = q;
5631:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
5632:     PetscInt fdof;

5634:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
5635:     if (!dof) continue;
5636:     for (f = 0; f < numFields; ++f) {
5637:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
5638:       foffsets[f+1] += fdof;
5639:     }
5640:     numFIndices += dof;
5641:   }
5642:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

5644:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
5645:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
5646:   if (numFields) {
5647:     for (p = 0; p < numFPoints*2; p += 2) {
5648:       PetscInt o = ftotpoints[p+1];
5649:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5650:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
5651:     }
5652:     for (p = 0; p < numCPoints*2; p += 2) {
5653:       PetscInt o = cpoints[p+1];
5654:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5655:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5656:     }
5657:   } else {
5658:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5659:       PetscInt o = ftotpoints[p+1];
5660:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5661:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5662:     }
5663:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5664:       PetscInt o = cpoints[p+1];
5665:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5666:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5667:     }
5668:   }
5669:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5670:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5671:   return(0);
5672: }

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

5679:   Input Parameter:
5680: . dm - The DMPlex object

5682:   Output Parameters:
5683: + cMax - The first hybrid cell
5684: . fMax - The first hybrid face
5685: . eMax - The first hybrid edge
5686: - vMax - The first hybrid vertex

5688:   Level: developer

5690: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
5691: @*/
5692: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
5693: {
5694:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5695:   PetscInt       dim;

5700:   DMPlexGetDimension(dm, &dim);
5701:   if (cMax) *cMax = mesh->hybridPointMax[dim];
5702:   if (fMax) *fMax = mesh->hybridPointMax[dim-1];
5703:   if (eMax) *eMax = mesh->hybridPointMax[1];
5704:   if (vMax) *vMax = mesh->hybridPointMax[0];
5705:   return(0);
5706: }

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

5713:   Input Parameters:
5714: . dm   - The DMPlex object
5715: . cMax - The first hybrid cell
5716: . fMax - The first hybrid face
5717: . eMax - The first hybrid edge
5718: - vMax - The first hybrid vertex

5720:   Level: developer

5722: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5723: @*/
5724: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5725: {
5726:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5727:   PetscInt       dim;

5732:   DMPlexGetDimension(dm, &dim);
5733:   if (cMax >= 0) mesh->hybridPointMax[dim]   = cMax;
5734:   if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5735:   if (eMax >= 0) mesh->hybridPointMax[1]     = eMax;
5736:   if (vMax >= 0) mesh->hybridPointMax[0]     = vMax;
5737:   return(0);
5738: }

5742: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5743: {
5744:   DM_Plex *mesh = (DM_Plex*) dm->data;

5749:   *cellHeight = mesh->vtkCellHeight;
5750:   return(0);
5751: }

5755: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5756: {
5757:   DM_Plex *mesh = (DM_Plex*) dm->data;

5761:   mesh->vtkCellHeight = cellHeight;
5762:   return(0);
5763: }

5767: /* We can easily have a form that takes an IS instead */
5768: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5769: {
5770:   PetscSection   section, globalSection;
5771:   PetscInt      *numbers, p;

5775:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
5776:   PetscSectionSetChart(section, pStart, pEnd);
5777:   for (p = pStart; p < pEnd; ++p) {
5778:     PetscSectionSetDof(section, p, 1);
5779:   }
5780:   PetscSectionSetUp(section);
5781:   PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, &globalSection);
5782:   PetscMalloc1((pEnd - pStart), &numbers);
5783:   for (p = pStart; p < pEnd; ++p) {
5784:     PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5785:     if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5786:     else                       numbers[p-pStart] += shift;
5787:   }
5788:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5789:   if (globalSize) {
5790:     PetscLayout layout;
5791:     PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5792:     PetscLayoutGetSize(layout, globalSize);
5793:     PetscLayoutDestroy(&layout);
5794:   }
5795:   PetscSectionDestroy(&section);
5796:   PetscSectionDestroy(&globalSection);
5797:   return(0);
5798: }

5802: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5803: {
5804:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5805:   PetscInt       cellHeight, cStart, cEnd, cMax;

5810:   if (!mesh->globalCellNumbers) {
5811:     DMPlexGetVTKCellHeight(dm, &cellHeight);
5812:     DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5813:     DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5814:     if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
5815:     DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, &mesh->globalCellNumbers);
5816:   }
5817:   *globalCellNumbers = mesh->globalCellNumbers;
5818:   return(0);
5819: }

5823: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5824: {
5825:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5826:   PetscInt       vStart, vEnd, vMax;

5831:   if (!mesh->globalVertexNumbers) {
5832:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5833:     DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5834:     if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
5835:     DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, &mesh->globalVertexNumbers);
5836:   }
5837:   *globalVertexNumbers = mesh->globalVertexNumbers;
5838:   return(0);
5839: }

5843: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5844: {
5845:   IS             nums[4];
5846:   PetscInt       depths[4];
5847:   PetscInt       depth, d, shift = 0;

5852:   DMPlexGetDepth(dm, &depth);
5853:   depths[0] = depth; depths[1] = 0;
5854:   for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5855:   for (d = 0; d <= depth; ++d) {
5856:     PetscInt pStart, pEnd, gsize;

5858:     DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5859:     DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5860:     shift += gsize;
5861:   }
5862:   ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5863:   for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5864:   return(0);
5865: }


5870: /*@C
5871:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
5872:   the local section and an SF describing the section point overlap.

5874:   Input Parameters:
5875:   + s - The PetscSection for the local field layout
5876:   . sf - The SF describing parallel layout of the section points
5877:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
5878:   . label - The label specifying the points
5879:   - labelValue - The label stratum specifying the points

5881:   Output Parameter:
5882:   . gsection - The PetscSection for the global field layout

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

5886:   Level: developer

5888: .seealso: PetscSectionCreate()
5889: @*/
5890: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
5891: {
5892:   PetscInt      *neg = NULL, *tmpOff = NULL;
5893:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

5897:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
5898:   PetscSectionGetChart(s, &pStart, &pEnd);
5899:   PetscSectionSetChart(*gsection, pStart, pEnd);
5900:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
5901:   if (nroots >= 0) {
5902:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
5903:     PetscCalloc1(nroots, &neg);
5904:     if (nroots > pEnd-pStart) {
5905:       PetscCalloc1(nroots, &tmpOff);
5906:     } else {
5907:       tmpOff = &(*gsection)->atlasDof[-pStart];
5908:     }
5909:   }
5910:   /* Mark ghost points with negative dof */
5911:   for (p = pStart; p < pEnd; ++p) {
5912:     PetscInt value;

5914:     DMLabelGetValue(label, p, &value);
5915:     if (value != labelValue) continue;
5916:     PetscSectionGetDof(s, p, &dof);
5917:     PetscSectionSetDof(*gsection, p, dof);
5918:     PetscSectionGetConstraintDof(s, p, &cdof);
5919:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
5920:     if (neg) neg[p] = -(dof+1);
5921:   }
5922:   PetscSectionSetUpBC(*gsection);
5923:   if (nroots >= 0) {
5924:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5925:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5926:     if (nroots > pEnd-pStart) {
5927:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
5928:     }
5929:   }
5930:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
5931:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5932:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
5933:     (*gsection)->atlasOff[p] = off;
5934:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
5935:   }
5936:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
5937:   globalOff -= off;
5938:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5939:     (*gsection)->atlasOff[p] += globalOff;
5940:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
5941:   }
5942:   /* Put in negative offsets for ghost points */
5943:   if (nroots >= 0) {
5944:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5945:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5946:     if (nroots > pEnd-pStart) {
5947:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
5948:     }
5949:   }
5950:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
5951:   PetscFree(neg);
5952:   return(0);
5953: }

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

5960:   Input Parameters:
5961:   + dm - The DMPlex object

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

5965:   Level: developer

5967: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5968: @*/
5969: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5970: {
5971:   PetscSection    coneSection, supportSection;
5972:   const PetscInt *cone, *support;
5973:   PetscInt        coneSize, c, supportSize, s;
5974:   PetscInt        pStart, pEnd, p, csize, ssize;
5975:   PetscErrorCode  ierr;

5979:   DMPlexGetConeSection(dm, &coneSection);
5980:   DMPlexGetSupportSection(dm, &supportSection);
5981:   /* Check that point p is found in the support of its cone points, and vice versa */
5982:   DMPlexGetChart(dm, &pStart, &pEnd);
5983:   for (p = pStart; p < pEnd; ++p) {
5984:     DMPlexGetConeSize(dm, p, &coneSize);
5985:     DMPlexGetCone(dm, p, &cone);
5986:     for (c = 0; c < coneSize; ++c) {
5987:       PetscBool dup = PETSC_FALSE;
5988:       PetscInt  d;
5989:       for (d = c-1; d >= 0; --d) {
5990:         if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5991:       }
5992:       DMPlexGetSupportSize(dm, cone[c], &supportSize);
5993:       DMPlexGetSupport(dm, cone[c], &support);
5994:       for (s = 0; s < supportSize; ++s) {
5995:         if (support[s] == p) break;
5996:       }
5997:       if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5998:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5999:         for (s = 0; s < coneSize; ++s) {
6000:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
6001:         }
6002:         PetscPrintf(PETSC_COMM_SELF, "\n");
6003:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
6004:         for (s = 0; s < supportSize; ++s) {
6005:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
6006:         }
6007:         PetscPrintf(PETSC_COMM_SELF, "\n");
6008:         if (dup) {
6009:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
6010:         } else {
6011:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
6012:         }
6013:       }
6014:     }
6015:     DMPlexGetSupportSize(dm, p, &supportSize);
6016:     DMPlexGetSupport(dm, p, &support);
6017:     for (s = 0; s < supportSize; ++s) {
6018:       DMPlexGetConeSize(dm, support[s], &coneSize);
6019:       DMPlexGetCone(dm, support[s], &cone);
6020:       for (c = 0; c < coneSize; ++c) {
6021:         if (cone[c] == p) break;
6022:       }
6023:       if (c >= coneSize) {
6024:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
6025:         for (c = 0; c < supportSize; ++c) {
6026:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
6027:         }
6028:         PetscPrintf(PETSC_COMM_SELF, "\n");
6029:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
6030:         for (c = 0; c < coneSize; ++c) {
6031:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
6032:         }
6033:         PetscPrintf(PETSC_COMM_SELF, "\n");
6034:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
6035:       }
6036:     }
6037:   }
6038:   PetscSectionGetStorageSize(coneSection, &csize);
6039:   PetscSectionGetStorageSize(supportSection, &ssize);
6040:   if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
6041:   return(0);
6042: }

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

6049:   Input Parameters:
6050: + dm - The DMPlex object
6051: . isSimplex - Are the cells simplices or tensor products
6052: - cellHeight - Normally 0

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

6056:   Level: developer

6058: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
6059: @*/
6060: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
6061: {
6062:   PetscInt       dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;

6067:   DMPlexGetDimension(dm, &dim);
6068:   switch (dim) {
6069:   case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
6070:   case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
6071:   case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
6072:   default:
6073:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
6074:   }
6075:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
6076:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
6077:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
6078:   cMax = cMax >= 0 ? cMax : cEnd;
6079:   for (c = cStart; c < cMax; ++c) {
6080:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

6082:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6083:     for (cl = 0; cl < closureSize*2; cl += 2) {
6084:       const PetscInt p = closure[cl];
6085:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
6086:     }
6087:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6088:     if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has  %d vertices != %d", c, coneSize, numCorners);
6089:   }
6090:   for (c = cMax; c < cEnd; ++c) {
6091:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

6093:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6094:     for (cl = 0; cl < closureSize*2; cl += 2) {
6095:       const PetscInt p = closure[cl];
6096:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
6097:     }
6098:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6099:     if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has  %d vertices > %d", c, coneSize, numHybridCorners);
6100:   }
6101:   return(0);
6102: }

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

6109:   Input Parameters:
6110: + dm - The DMPlex object
6111: . isSimplex - Are the cells simplices or tensor products
6112: - cellHeight - Normally 0

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

6116:   Level: developer

6118: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
6119: @*/
6120: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
6121: {
6122:   PetscInt       pMax[4];
6123:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, h;

6128:   DMPlexGetDimension(dm, &dim);
6129:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
6130:   DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
6131:   for (h = cellHeight; h < dim; ++h) {
6132:     DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
6133:     for (c = cStart; c < cEnd; ++c) {
6134:       const PetscInt *cone, *ornt, *faces;
6135:       PetscInt        numFaces, faceSize, coneSize,f;
6136:       PetscInt       *closure = NULL, closureSize, cl, numCorners = 0;

6138:       if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
6139:       DMPlexGetConeSize(dm, c, &coneSize);
6140:       DMPlexGetCone(dm, c, &cone);
6141:       DMPlexGetConeOrientation(dm, c, &ornt);
6142:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6143:       for (cl = 0; cl < closureSize*2; cl += 2) {
6144:         const PetscInt p = closure[cl];
6145:         if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
6146:       }
6147:       DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
6148:       if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
6149:       for (f = 0; f < numFaces; ++f) {
6150:         PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;

6152:         DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
6153:         for (cl = 0; cl < fclosureSize*2; cl += 2) {
6154:           const PetscInt p = fclosure[cl];
6155:           if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
6156:         }
6157:         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);
6158:         for (v = 0; v < fnumCorners; ++v) {
6159:           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]);
6160:         }
6161:         DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
6162:       }
6163:       DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
6164:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6165:     }
6166:   }
6167:   return(0);
6168: }

6172: /* Pointwise interpolation
6173:      Just code FEM for now
6174:      u^f = I u^c
6175:      sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
6176:      u^f_i = sum_j psi^f_i I phi^c_j u^c_j
6177:      I_{ij} = psi^f_i phi^c_j
6178: */
6179: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
6180: {
6181:   PetscSection   gsc, gsf;
6182:   PetscInt       m, n;
6183:   void          *ctx;

6187:   /*
6188:   Loop over coarse cells
6189:     Loop over coarse basis functions
6190:       Loop over fine cells in coarse cell
6191:         Loop over fine dual basis functions
6192:           Evaluate coarse basis on fine dual basis quad points
6193:           Sum
6194:           Update local element matrix
6195:     Accumulate to interpolation matrix

6197:    Can extend PetscFEIntegrateJacobian_Basic() to do a specialized cell loop
6198:   */
6199:   DMGetDefaultGlobalSection(dmFine, &gsf);
6200:   PetscSectionGetConstrainedStorageSize(gsf, &m);
6201:   DMGetDefaultGlobalSection(dmCoarse, &gsc);
6202:   PetscSectionGetConstrainedStorageSize(gsc, &n);
6203:   /* We need to preallocate properly */
6204:   MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
6205:   MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
6206:   MatSetType(*interpolation, dmCoarse->mattype);
6207:   DMGetApplicationContext(dmFine, &ctx);
6208:   DMPlexComputeInterpolatorFEM(dmCoarse, dmFine, *interpolation, ctx);
6209:   /* Use naive scaling */
6210:   DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
6211:   return(0);
6212: }

6216: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, VecScatter *ctx)
6217: {

6221:   DMPlexComputeInjectorFEM(dmCoarse, dmFine, ctx, NULL);
6222:   return(0);
6223: }

6227: /* Pointwise interpolation
6228:      Just code FEM for now
6229:      u^f = I u^c
6230:      sum_k u^f_k phi^f_k = I sum_l u^c_l phi^c_l
6231:      u^f_i = sum_l int psi^f_i I phi^c_l u^c_l
6232:      I_{ij} = int psi^f_i phi^c_j
6233: */
6234: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
6235: {
6236:   PetscSection   section;
6237:   IS            *bcPoints;
6238:   PetscInt      *bcFields, *numComp, *numDof;
6239:   PetscInt       depth, dim, numBd, numBC = 0, numFields, bd, bc, f;

6243:   /* Handle boundary conditions */
6244:   DMPlexGetDepth(dm, &depth);
6245:   DMPlexGetDimension(dm, &dim);
6246:   DMPlexGetNumBoundary(dm, &numBd);
6247:   for (bd = 0; bd < numBd; ++bd) {
6248:     PetscBool isEssential;
6249:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
6250:     if (isEssential) ++numBC;
6251:   }
6252:   PetscMalloc2(numBC,&bcFields,numBC,&bcPoints);
6253:   for (bd = 0, bc = 0; bd < numBd; ++bd) {
6254:     const char     *bdLabel;
6255:     DMLabel         label;
6256:     const PetscInt *values;
6257:     PetscInt        bd2, field, numValues;
6258:     PetscBool       isEssential, duplicate = PETSC_FALSE;

6260:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);
6261:     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
6262:     DMPlexGetLabel(dm, bdLabel, &label);
6263:     /* Only want to do this for FEM, and only once */
6264:     for (bd2 = 0; bd2 < bd; ++bd2) {
6265:       const char *bdname;
6266:       DMPlexGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL);
6267:       PetscStrcmp(bdname, bdLabel, &duplicate);
6268:       if (duplicate) break;
6269:     }
6270:     if (!duplicate) {
6271:       DMPlexLabelComplete(dm, label);
6272:       DMPlexLabelAddCells(dm, label);
6273:     }
6274:     /* Filter out cells, if you actually want to constraint cells you need to do things by hand right now */
6275:     if (isEssential) {
6276:       IS              tmp;
6277:       PetscInt       *newidx;
6278:       const PetscInt *idx;
6279:       PetscInt        cStart, cEnd, n, p, newn = 0;

6281:       bcFields[bc] = field;
6282:       DMPlexGetStratumIS(dm, bdLabel, values[0], &tmp);
6283:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
6284:       ISGetLocalSize(tmp, &n);
6285:       ISGetIndices(tmp, &idx);
6286:       for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
6287:       PetscMalloc1(newn,&newidx);
6288:       newn = 0;
6289:       for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
6290:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
6291:       ISRestoreIndices(tmp, &idx);
6292:       ISDestroy(&tmp);
6293:     }
6294:   }
6295:   /* Handle discretization */
6296:   DMGetNumFields(dm, &numFields);
6297:   PetscMalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
6298:   for (f = 0; f < numFields; ++f) {
6299:     PetscFE         fe;
6300:     const PetscInt *numFieldDof;
6301:     PetscInt        d;

6303:     DMGetField(dm, f, (PetscObject *) &fe);
6304:     PetscFEGetNumComponents(fe, &numComp[f]);
6305:     PetscFEGetNumDof(fe, &numFieldDof);
6306:     for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
6307:   }
6308:   for (f = 0; f < numFields; ++f) {
6309:     PetscInt d;
6310:     for (d = 1; d < dim; ++d) {
6311:       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.");
6312:     }
6313:   }
6314:   DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, NULL, &section);
6315:   for (f = 0; f < numFields; ++f) {
6316:     PetscFE     fe;
6317:     const char *name;

6319:     DMGetField(dm, f, (PetscObject *) &fe);
6320:     PetscObjectGetName((PetscObject) fe, &name);
6321:     PetscSectionSetFieldName(section, f, name);
6322:   }
6323:   DMSetDefaultSection(dm, section);
6324:   PetscSectionDestroy(&section);
6325:   for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);}
6326:   PetscFree2(bcFields,bcPoints);
6327:   PetscFree2(numComp,numDof);
6328:   return(0);
6329: }

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

6336:   Input Parameter:
6337: . dm - The DMPlex object

6339:   Output Parameter:
6340: . cdm - The coarse DM

6342:   Level: intermediate

6344: .seealso: DMPlexSetCoarseDM()
6345: @*/
6346: PetscErrorCode DMPlexGetCoarseDM(DM dm, DM *cdm)
6347: {
6351:   *cdm = ((DM_Plex *) dm->data)->coarseMesh;
6352:   return(0);
6353: }

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

6360:   Input Parameters:
6361: + dm - The DMPlex object
6362: - cdm - The coarse DM

6364:   Level: intermediate

6366: .seealso: DMPlexGetCoarseDM()
6367: @*/
6368: PetscErrorCode DMPlexSetCoarseDM(DM dm, DM cdm)
6369: {
6370:   DM_Plex       *mesh;

6376:   mesh = (DM_Plex *) dm->data;
6377:   DMDestroy(&mesh->coarseMesh);
6378:   mesh->coarseMesh = cdm;
6379:   PetscObjectReference((PetscObject) mesh->coarseMesh);
6380:   return(0);
6381: }