Actual source code: plex.c

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

  6: /* Logging support */
  7: PetscLogEvent DMPLEX_Interpolate, PETSCPARTITIONER_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeOverlap, DMPLEX_DistributeField, DMPLEX_DistributeData, DMPLEX_Migrate, DMPLEX_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM, DMPLEX_InterpolatorFEM, DMPLEX_InjectorFEM, DMPLEX_IntegralFEM, DMPLEX_CreateGmsh;

  9: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
 10: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
 11: PETSC_EXTERN PetscErrorCode VecLoad_Default(Vec, PetscViewer);

 15: PetscErrorCode DMPlexGetFieldType_Internal(DM dm, PetscSection section, PetscInt field, PetscInt *sStart, PetscInt *sEnd, PetscViewerVTKFieldType *ft)
 16: {
 17:   PetscInt       dim, pStart, pEnd, vStart, vEnd, cStart, cEnd, cEndInterior, vdof = 0, cdof = 0;

 21:   *ft  = PETSC_VTK_POINT_FIELD;
 22:   DMGetDimension(dm, &dim);
 23:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 24:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 25:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
 26:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
 27:   PetscSectionGetChart(section, &pStart, &pEnd);
 28:   if (field >= 0) {
 29:     if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetFieldDof(section, vStart, field, &vdof);}
 30:     if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetFieldDof(section, cStart, field, &cdof);}
 31:   } else {
 32:     if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetDof(section, vStart, &vdof);}
 33:     if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetDof(section, cStart, &cdof);}
 34:   }
 35:   if (vdof) {
 36:     *sStart = vStart;
 37:     *sEnd   = vEnd;
 38:     if (vdof == dim) *ft = PETSC_VTK_POINT_VECTOR_FIELD;
 39:     else             *ft = PETSC_VTK_POINT_FIELD;
 40:   } else if (cdof) {
 41:     *sStart = cStart;
 42:     *sEnd   = cEnd;
 43:     if (cdof == dim) *ft = PETSC_VTK_CELL_VECTOR_FIELD;
 44:     else             *ft = PETSC_VTK_CELL_FIELD;
 45:   } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Could not classify input Vec for VTK");
 46:   return(0);
 47: }

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

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

 67:     DMGetNumFields(dm, &numFields);
 68:     if (numFields) {
 69:       PetscObject fe;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

565: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
566: {
567:   PetscBool      iascii, ishdf5, isvtk;

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

595: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
596: {
597:   PetscBool      isbinary, ishdf5;

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

618: static PetscErrorCode BoundaryDestroy(DMBoundary *boundary)
619: {
620:   DMBoundary     b, next;

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

640: PetscErrorCode DMDestroy_Plex(DM dm)
641: {
642:   DM_Plex       *mesh = (DM_Plex*) dm->data;
643:   PlexLabel      next = mesh->labels;

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

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

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

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

715:     if (bs < 0) {
716:       if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
717:         PetscInt pStart, pEnd, p, dof, cdof;

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

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

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

763:   Not collective

765:   Input Parameter:
766: . mesh - The DMPlex

768:   Output Parameters:
769: + pStart - The first mesh point
770: - pEnd   - The upper bound for mesh points

772:   Level: beginner

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

783:   PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
784:   return(0);
785: }

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

792:   Not collective

794:   Input Parameters:
795: + mesh - The DMPlex
796: . pStart - The first mesh point
797: - pEnd   - The upper bound for mesh points

799:   Output Parameters:

801:   Level: beginner

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

812:   PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
813:   PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
814:   return(0);
815: }

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

822:   Not collective

824:   Input Parameters:
825: + mesh - The DMPlex
826: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

828:   Output Parameter:
829: . size - The cone size for point p

831:   Level: beginner

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

843:   PetscSectionGetDof(mesh->coneSection, p, size);
844:   return(0);
845: }

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

852:   Not collective

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

859:   Output Parameter:

861:   Note:
862:   This should be called after DMPlexSetChart().

864:   Level: beginner

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

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

877:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
878:   return(0);
879: }

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

886:   Not collective

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

893:   Output Parameter:

895:   Note:
896:   This should be called after DMPlexSetChart().

898:   Level: beginner

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

910:   PetscSectionAddDof(mesh->coneSection, p, size);
911:   PetscSectionGetDof(mesh->coneSection, p, &csize);

913:   mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
914:   return(0);
915: }

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

922:   Not collective

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

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

931:   Level: beginner

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

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

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

950:   PetscSectionGetOffset(mesh->coneSection, p, &off);
951:   *cone = &mesh->cones[off];
952:   return(0);
953: }

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

960:   Not collective

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

967:   Output Parameter:

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

972:   Level: beginner

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

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

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

1002:   Not collective

1004:   Input Parameters:
1005: + mesh - The DMPlex
1006: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

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

1014:   Level: beginner

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

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

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

1032: #if defined(PETSC_USE_DEBUG)
1033:   {
1034:     PetscInt dof;
1035:     PetscSectionGetDof(mesh->coneSection, p, &dof);
1037:   }
1038: #endif
1039:   PetscSectionGetOffset(mesh->coneSection, p, &off);

1041:   *coneOrientation = &mesh->coneOrientations[off];
1042:   return(0);
1043: }

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

1050:   Not collective

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

1060:   Output Parameter:

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

1065:   Level: beginner

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

1078:   PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1079:   PetscSectionGetDof(mesh->coneSection, p, &dof);
1081:   PetscSectionGetOffset(mesh->coneSection, p, &off);
1082:   if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1083:   for (c = 0; c < dof; ++c) {
1084:     PetscInt cdof, o = coneOrientation[c];

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

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

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

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

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

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

1139:   Not collective

1141:   Input Parameters:
1142: + mesh - The DMPlex
1143: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()

1145:   Output Parameter:
1146: . size - The support size for point p

1148:   Level: beginner

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

1160:   PetscSectionGetDof(mesh->supportSection, p, size);
1161:   return(0);
1162: }

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

1169:   Not collective

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

1176:   Output Parameter:

1178:   Note:
1179:   This should be called after DMPlexSetChart().

1181:   Level: beginner

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

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

1194:   mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1195:   return(0);
1196: }

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

1203:   Not collective

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

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

1212:   Level: beginner

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

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

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

1231:   PetscSectionGetOffset(mesh->supportSection, p, &off);
1232:   *support = &mesh->supports[off];
1233:   return(0);
1234: }

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

1241:   Not collective

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

1248:   Output Parameter:

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

1253:   Level: beginner

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

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

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

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

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

1304:   Not collective

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

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

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

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

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

1325:   Level: beginner

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

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

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

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

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

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

1446: /*@C
1447:   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

1449:   Not collective

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

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

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

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

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

1471:   Level: beginner

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

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

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

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

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

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

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

1603:   Not collective

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

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

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

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

1621:   Level: beginner

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

1633:   DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1634:   if (numPoints) *numPoints = 0;
1635:   return(0);
1636: }

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

1643:   Not collective

1645:   Input Parameter:
1646: . mesh - The DMPlex

1648:   Output Parameters:
1649: + maxConeSize - The maximum number of in-edges
1650: - maxSupportSize - The maximum number of out-edges

1652:   Level: beginner

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

1662:   if (maxConeSize)    *maxConeSize    = mesh->maxConeSize;
1663:   if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1664:   return(0);
1665: }

1669: PetscErrorCode DMSetUp_Plex(DM dm)
1670: {
1671:   DM_Plex       *mesh = (DM_Plex*) dm->data;
1672:   PetscInt       size;

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

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

1696:   if (subdm) {DMClone(dm, subdm);}
1697:   DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1698:   return(0);
1699: }

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

1706:   Not collective

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

1711:   Output Parameter:

1713:   Note:
1714:   This should be called after all calls to DMPlexSetCone()

1716:   Level: beginner

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

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

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

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

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

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

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

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

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

1781:   Not collective

1783:   Input Parameter:
1784: . mesh - The DMPlex

1786:   Output Parameter:

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

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

1796:   Level: beginner

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

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

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

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

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

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

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

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

1880:   Not Collective

1882:   Input Parameters:
1883: + dm - The DMPlex object
1884: . numPoints - The number of input points for the join
1885: - points - The input points

1887:   Output Parameters:
1888: + numCoveredPoints - The number of points in the join
1889: - coveredPoints - The points in the join

1891:   Level: intermediate

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

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

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

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

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

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

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

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

1955:   Not Collective

1957:   Input Parameters:
1958: + dm - The DMPlex object
1959: . numPoints - The number of input points for the join
1960: - points - The input points

1962:   Output Parameters:
1963: + numCoveredPoints - The number of points in the join
1964: - coveredPoints - The points in the join

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

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

1972:   Level: intermediate

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

1986:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
1987:   if (numCoveredPoints) *numCoveredPoints = 0;
1988:   return(0);
1989: }

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

1996:   Not Collective

1998:   Input Parameters:
1999: + dm - The DMPlex object
2000: . numPoints - The number of input points for the join
2001: - points - The input points

2003:   Output Parameters:
2004: + numCoveredPoints - The number of points in the join
2005: - coveredPoints - The points in the join

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

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

2013:   Level: intermediate

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


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

2041:   for (p = 0; p < numPoints; ++p) {
2042:     PetscInt closureSize;

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

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

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

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

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

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

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

2105:   Not Collective

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

2112:   Output Parameters:
2113: + numCoveredPoints - The number of points in the meet
2114: - coveredPoints - The points in the meet

2116:   Level: intermediate

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

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: DMPlexRestoreMeet(), DMPlexGetJoin()
2128: @*/
2129: PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2130: {
2131:   DM_Plex       *mesh = (DM_Plex*) dm->data;
2132:   PetscInt      *meet[2];
2133:   PetscInt       meetSize, i = 0;
2134:   PetscInt       dof, off, p, c, m;

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

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

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

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

2180:   Not Collective

2182:   Input Parameters:
2183: + dm - The DMPlex object
2184: . numPoints - The number of input points for the meet
2185: - points - The input points

2187:   Output Parameters:
2188: + numCoveredPoints - The number of points in the meet
2189: - coveredPoints - The points in the meet

2191:   Level: intermediate

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

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

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

2211:   DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2212:   if (numCoveredPoints) *numCoveredPoints = 0;
2213:   return(0);
2214: }

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

2221:   Not Collective

2223:   Input Parameters:
2224: + dm - The DMPlex object
2225: . numPoints - The number of input points for the meet
2226: - points - The input points

2228:   Output Parameters:
2229: + numCoveredPoints - The number of points in the meet
2230: - coveredPoints - The points in the meet

2232:   Level: intermediate

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

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

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


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

2266:   for (p = 0; p < numPoints; ++p) {
2267:     PetscInt closureSize;

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

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

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

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

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

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

2327: /*@C
2328:   DMPlexEqual - Determine if two DMs have the same topology

2330:   Not Collective

2332:   Input Parameters:
2333: + dmA - A DMPlex object
2334: - dmB - A DMPlex object

2336:   Output Parameters:
2337: . equal - PETSC_TRUE if the topologies are identical

2339:   Level: intermediate

2341:   Notes:
2342:   We are not solving graph isomorphism, so we do not permutation.

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


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

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

2394: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2395: {
2396:   MPI_Comm       comm;

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

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

2468:   Input Parameters:
2469: + dm     - The DM
2470: - in     - The input coordinate point (dim numbers)

2472:   Output Parameter:
2473: . out - The localized coordinate point

2475:   Level: developer

2477: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeAddCoordinate()
2478: @*/
2479: PetscErrorCode DMPlexLocalizeCoordinate(DM dm, const PetscScalar in[], PetscScalar out[])
2480: {
2481:   PetscInt       dim, d;

2485:   DMGetCoordinateDim(dm, &dim);
2486:   if (!dm->maxCell) {
2487:     for (d = 0; d < dim; ++d) out[d] = in[d];
2488:   } else {
2489:     for (d = 0; d < dim; ++d) {
2490:       out[d] = in[d] - dm->L[d]*floor(PetscRealPart(in[d])/dm->L[d]);
2491:     }
2492:   }
2493:   return(0);
2494: }

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

2501:   Input Parameters:
2502: + dm     - The DM
2503: . dim    - The spatial dimension
2504: . anchor - The anchor point, the input point can be no more than maxCell away from it
2505: - in     - The input coordinate point (dim numbers)

2507:   Output Parameter:
2508: . out - The localized coordinate point

2510:   Level: developer

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

2514: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeAddCoordinate()
2515: */
2516: PetscErrorCode DMPlexLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2517: {
2518:   PetscInt d;

2521:   if (!dm->maxCell) {
2522:     for (d = 0; d < dim; ++d) out[d] = in[d];
2523:   } else {
2524:     for (d = 0; d < dim; ++d) {
2525:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2526:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2527:       } else {
2528:         out[d] = in[d];
2529:       }
2530:     }
2531:   }
2532:   return(0);
2533: }
2536: PetscErrorCode DMPlexLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
2537: {
2538:   PetscInt d;

2541:   if (!dm->maxCell) {
2542:     for (d = 0; d < dim; ++d) out[d] = in[d];
2543:   } else {
2544:     for (d = 0; d < dim; ++d) {
2545:       if (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d]) {
2546:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
2547:       } else {
2548:         out[d] = in[d];
2549:       }
2550:     }
2551:   }
2552:   return(0);
2553: }

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

2560:   Input Parameters:
2561: + dm     - The DM
2562: . dim    - The spatial dimension
2563: . anchor - The anchor point, the input point can be no more than maxCell away from it
2564: . in     - The input coordinate delta (dim numbers)
2565: - out    - The input coordinate point (dim numbers)

2567:   Output Parameter:
2568: . out    - The localized coordinate in + out

2570:   Level: developer

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

2574: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeCoordinate()
2575: */
2576: PetscErrorCode DMPlexLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2577: {
2578:   PetscInt d;

2581:   if (!dm->maxCell) {
2582:     for (d = 0; d < dim; ++d) out[d] += in[d];
2583:   } else {
2584:     for (d = 0; d < dim; ++d) {
2585:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2586:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2587:       } else {
2588:         out[d] += in[d];
2589:       }
2590:     }
2591:   }
2592:   return(0);
2593: }

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

2600:   Input Parameter:
2601: . dm - The DM

2603:   Level: developer

2605: .seealso: DMPlexLocalizeCoordinate(), DMPlexLocalizeAddCoordinate()
2606: @*/
2607: PetscErrorCode DMPlexLocalizeCoordinates(DM dm)
2608: {
2609:   PetscSection   coordSection, cSection;
2610:   Vec            coordinates,  cVec;
2611:   PetscScalar   *coords, *coords2, *anchor;
2612:   PetscInt       Nc, cStart, cEnd, c, vStart, vEnd, v, dof, d, off, off2, bs, coordSize;

2617:   if (!dm->maxCell) return(0);
2618:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2619:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2620:   DMGetCoordinatesLocal(dm, &coordinates);
2621:   DMGetCoordinateSection(dm, &coordSection);
2622:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
2623:   PetscSectionSetNumFields(cSection, 1);
2624:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
2625:   PetscSectionSetFieldComponents(cSection, 0, Nc);
2626:   PetscSectionSetChart(cSection, cStart, vEnd);
2627:   for (v = vStart; v < vEnd; ++v) {
2628:     PetscSectionGetDof(coordSection, v, &dof);
2629:     PetscSectionSetDof(cSection,     v,  dof);
2630:     PetscSectionSetFieldDof(cSection, v, 0, dof);
2631:   }
2632:   for (c = cStart; c < cEnd; ++c) {
2633:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, NULL);
2634:     PetscSectionSetDof(cSection, c, dof);
2635:     PetscSectionSetFieldDof(cSection, c, 0, dof);
2636:   }
2637:   PetscSectionSetUp(cSection);
2638:   PetscSectionGetStorageSize(cSection, &coordSize);
2639:   VecCreate(PetscObjectComm((PetscObject) dm), &cVec);
2640:   VecGetBlockSize(coordinates, &bs);
2641:   VecSetBlockSize(cVec,         bs);
2642:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
2643:   VecSetType(cVec,VECSTANDARD);
2644:   VecGetArray(coordinates, &coords);
2645:   VecGetArray(cVec,        &coords2);
2646:   for (v = vStart; v < vEnd; ++v) {
2647:     PetscSectionGetDof(coordSection, v, &dof);
2648:     PetscSectionGetOffset(coordSection, v, &off);
2649:     PetscSectionGetOffset(cSection,     v, &off2);
2650:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
2651:   }
2652:   DMGetWorkArray(dm, 3, PETSC_SCALAR, &anchor);
2653:   for (c = cStart; c < cEnd; ++c) {
2654:     PetscScalar *cellCoords = NULL;
2655:     PetscInt     b;

2657:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2658:     PetscSectionGetOffset(cSection, c, &off2);
2659:     for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
2660:     for (d = 0; d < dof/bs; ++d) {DMPlexLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
2661:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2662:   }
2663:   DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);
2664:   VecRestoreArray(coordinates, &coords);
2665:   VecRestoreArray(cVec,        &coords2);
2666:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
2667:   DMSetCoordinatesLocal(dm, cVec);
2668:   VecDestroy(&cVec);
2669:   PetscSectionDestroy(&cSection);
2670:   return(0);
2671: }

2675: /*@
2676:   DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point

2678:   Not Collective

2680:   Input Parameter:
2681: . dm    - The DMPlex object

2683:   Output Parameter:
2684: . depthLabel - The DMLabel recording point depth

2686:   Level: developer

2688: .keywords: mesh, points
2689: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2690: @*/
2691: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
2692: {
2693:   DM_Plex       *mesh = (DM_Plex*) dm->data;

2699:   if (!mesh->depthLabel) {DMPlexGetLabel(dm, "depth", &mesh->depthLabel);}
2700:   *depthLabel = mesh->depthLabel;
2701:   return(0);
2702: }

2706: /*@
2707:   DMPlexGetDepth - Get the depth of the DAG representing this mesh

2709:   Not Collective

2711:   Input Parameter:
2712: . dm    - The DMPlex object

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

2717:   Level: developer

2719: .keywords: mesh, points
2720: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2721: @*/
2722: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
2723: {
2724:   DMLabel        label;
2725:   PetscInt       d = 0;

2731:   DMPlexGetDepthLabel(dm, &label);
2732:   if (label) {DMLabelGetNumValues(label, &d);}
2733:   *depth = d-1;
2734:   return(0);
2735: }

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

2742:   Not Collective

2744:   Input Parameters:
2745: + dm           - The DMPlex object
2746: - stratumValue - The requested depth

2748:   Output Parameters:
2749: + start - The first point at this depth
2750: - end   - One beyond the last point at this depth

2752:   Level: developer

2754: .keywords: mesh, points
2755: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
2756: @*/
2757: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2758: {
2759:   DMLabel        label;
2760:   PetscInt       pStart, pEnd;

2767:   DMPlexGetChart(dm, &pStart, &pEnd);
2768:   if (pStart == pEnd) return(0);
2769:   if (stratumValue < 0) {
2770:     if (start) *start = pStart;
2771:     if (end)   *end   = pEnd;
2772:     return(0);
2773:   }
2774:   DMPlexGetDepthLabel(dm, &label);
2775:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2776:   DMLabelGetStratumBounds(label, stratumValue, start, end);
2777:   return(0);
2778: }

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

2785:   Not Collective

2787:   Input Parameters:
2788: + dm           - The DMPlex object
2789: - stratumValue - The requested height

2791:   Output Parameters:
2792: + start - The first point at this height
2793: - end   - One beyond the last point at this height

2795:   Level: developer

2797: .keywords: mesh, points
2798: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
2799: @*/
2800: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2801: {
2802:   DMLabel        label;
2803:   PetscInt       depth, pStart, pEnd;

2810:   DMPlexGetChart(dm, &pStart, &pEnd);
2811:   if (pStart == pEnd) return(0);
2812:   if (stratumValue < 0) {
2813:     if (start) *start = pStart;
2814:     if (end)   *end   = pEnd;
2815:     return(0);
2816:   }
2817:   DMPlexGetDepthLabel(dm, &label);
2818:   if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2819:   DMLabelGetNumValues(label, &depth);
2820:   DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
2821:   return(0);
2822: }

2826: /* Set the number of dof on each point and separate by fields */
2827: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
2828: {
2829:   PetscInt      *pMax;
2830:   PetscInt       depth, pStart = 0, pEnd = 0;
2831:   PetscInt       Nf, p, d, dep, f;
2832:   PetscBool     *isFE;

2836:   PetscMalloc1(numFields, &isFE);
2837:   DMGetNumFields(dm, &Nf);
2838:   for (f = 0; f < numFields; ++f) {
2839:     PetscObject  obj;
2840:     PetscClassId id;

2842:     isFE[f] = PETSC_FALSE;
2843:     if (f >= Nf) continue;
2844:     DMGetField(dm, f, &obj);
2845:     PetscObjectGetClassId(obj, &id);
2846:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
2847:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
2848:   }
2849:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
2850:   if (numFields > 0) {
2851:     PetscSectionSetNumFields(*section, numFields);
2852:     if (numComp) {
2853:       for (f = 0; f < numFields; ++f) {
2854:         PetscSectionSetFieldComponents(*section, f, numComp[f]);
2855:       }
2856:     }
2857:   }
2858:   DMPlexGetChart(dm, &pStart, &pEnd);
2859:   PetscSectionSetChart(*section, pStart, pEnd);
2860:   DMPlexGetDepth(dm, &depth);
2861:   PetscMalloc1(depth+1,&pMax);
2862:   DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2863:   for (dep = 0; dep <= depth; ++dep) {
2864:     d    = dim == depth ? dep : (!dep ? 0 : dim);
2865:     DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
2866:     pMax[dep] = pMax[dep] < 0 ? pEnd : pMax[dep];
2867:     for (p = pStart; p < pEnd; ++p) {
2868:       PetscInt tot = 0;

2870:       for (f = 0; f < numFields; ++f) {
2871:         if (isFE[f] && p >= pMax[dep]) continue;
2872:         PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
2873:         tot += numDof[f*(dim+1)+d];
2874:       }
2875:       PetscSectionSetDof(*section, p, tot);
2876:     }
2877:   }
2878:   PetscFree(pMax);
2879:   PetscFree(isFE);
2880:   return(0);
2881: }

2885: /* Set the number of dof on each point and separate by fields
2886:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
2887: */
2888: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2889: {
2890:   PetscInt       numFields;
2891:   PetscInt       bc;
2892:   PetscSection   aSec;

2896:   PetscSectionGetNumFields(section, &numFields);
2897:   for (bc = 0; bc < numBC; ++bc) {
2898:     PetscInt        field = 0;
2899:     const PetscInt *comp;
2900:     const PetscInt *idx;
2901:     PetscInt        Nc = -1, n, i;

2903:     if (numFields) field = bcField[bc];
2904:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2905:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2906:     ISGetLocalSize(bcPoints[bc], &n);
2907:     ISGetIndices(bcPoints[bc], &idx);
2908:     for (i = 0; i < n; ++i) {
2909:       const PetscInt p = idx[i];
2910:       PetscInt       numConst;

2912:       if (numFields) {
2913:         PetscSectionGetFieldDof(section, p, field, &numConst);
2914:       } else {
2915:         PetscSectionGetDof(section, p, &numConst);
2916:       }
2917:       /* If Nc < 0, constrain every dof on the point */
2918:       if (Nc > 0) numConst = PetscMin(numConst, Nc);
2919:       if (numFields) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
2920:       PetscSectionAddConstraintDof(section, p, numConst);
2921:     }
2922:     ISRestoreIndices(bcPoints[bc], &idx);
2923:     if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2924:   }
2925:   DMPlexGetAnchors(dm, &aSec, NULL);
2926:   if (aSec) {
2927:     PetscInt aStart, aEnd, a;

2929:     PetscSectionGetChart(aSec, &aStart, &aEnd);
2930:     for (a = aStart; a < aEnd; a++) {
2931:       PetscInt dof, f;

2933:       PetscSectionGetDof(aSec, a, &dof);
2934:       if (dof) {
2935:         /* if there are point-to-point constraints, then all dofs are constrained */
2936:         PetscSectionGetDof(section, a, &dof);
2937:         PetscSectionSetConstraintDof(section, a, dof);
2938:         for (f = 0; f < numFields; f++) {
2939:           PetscSectionGetFieldDof(section, a, f, &dof);
2940:           PetscSectionSetFieldConstraintDof(section, a, f, dof);
2941:         }
2942:       }
2943:     }
2944:   }
2945:   return(0);
2946: }

2950: /* Set the constrained field indices on each point
2951:    If bcComps is NULL or the IS is NULL, constrain every dof on the point
2952: */
2953: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2954: {
2955:   PetscSection   aSec;
2956:   PetscInt      *indices;
2957:   PetscInt       numFields, maxDof, pStart, pEnd, p, bc, f, d;

2961:   PetscSectionGetNumFields(section, &numFields);
2962:   if (!numFields) return(0);
2963:   /* Initialize all field indices to -1 */
2964:   PetscSectionGetChart(section, &pStart, &pEnd);
2965:   PetscSectionGetMaxDof(section, &maxDof);
2966:   PetscMalloc1(maxDof, &indices);
2967:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
2968:   for (p = pStart; p < pEnd; ++p) for (f = 0; f < numFields; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
2969:   /* Handle BC constraints */
2970:   for (bc = 0; bc < numBC; ++bc) {
2971:     const PetscInt  field = bcField[bc];
2972:     const PetscInt *comp, *idx;
2973:     PetscInt        Nc = -1, n, i;

2975:     if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2976:     if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2977:     ISGetLocalSize(bcPoints[bc], &n);
2978:     ISGetIndices(bcPoints[bc], &idx);
2979:     for (i = 0; i < n; ++i) {
2980:       const PetscInt  p = idx[i];
2981:       const PetscInt *find;
2982:       PetscInt        fcdof, c;

2984:       PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
2985:       if (Nc < 0) {
2986:         for (d = 0; d < fcdof; ++d) indices[d] = d;
2987:       } else {
2988:         PetscSectionGetFieldConstraintIndices(section, p, field, &find);
2989:         for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
2990:         for (c = 0; c < Nc; ++c) indices[d+c] = comp[c];
2991:         PetscSortInt(d+Nc, indices);
2992:         for (c = d+Nc; c < fcdof; ++c) indices[c] = -1;
2993:       }
2994:       PetscSectionSetFieldConstraintIndices(section, p, field, indices);
2995:     }
2996:     if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2997:     ISRestoreIndices(bcPoints[bc], &idx);
2998:   }
2999:   /* Handle anchors */
3000:   DMPlexGetAnchors(dm, &aSec, NULL);
3001:   if (aSec) {
3002:     PetscInt aStart, aEnd, a;

3004:     for (d = 0; d < maxDof; ++d) indices[d] = d;
3005:     PetscSectionGetChart(aSec, &aStart, &aEnd);
3006:     for (a = aStart; a < aEnd; a++) {
3007:       PetscInt dof, fdof, f;

3009:       PetscSectionGetDof(aSec, a, &dof);
3010:       if (dof) {
3011:         /* if there are point-to-point constraints, then all dofs are constrained */
3012:         for (f = 0; f < numFields; f++) {
3013:           PetscSectionGetFieldDof(section, a, f, &fdof);
3014:           PetscSectionSetFieldConstraintIndices(section, a, f, indices);
3015:         }
3016:       }
3017:     }
3018:   }
3019:   PetscFree(indices);
3020:   return(0);
3021: }

3025: /* Set the constrained indices on each point */
3026: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
3027: {
3028:   PetscInt      *indices;
3029:   PetscInt       numFields, maxDof, pStart, pEnd, p, f, d;

3033:   PetscSectionGetNumFields(section, &numFields);
3034:   PetscSectionGetMaxDof(section, &maxDof);
3035:   PetscSectionGetChart(section, &pStart, &pEnd);
3036:   PetscMalloc1(maxDof, &indices);
3037:   for (d = 0; d < maxDof; ++d) indices[d] = -1;
3038:   for (p = pStart; p < pEnd; ++p) {
3039:     PetscInt cdof, d;

3041:     PetscSectionGetConstraintDof(section, p, &cdof);
3042:     if (cdof) {
3043:       if (numFields) {
3044:         PetscInt numConst = 0, foff = 0;

3046:         for (f = 0; f < numFields; ++f) {
3047:           const PetscInt *find;
3048:           PetscInt        fcdof, fdof;

3050:           PetscSectionGetFieldDof(section, p, f, &fdof);
3051:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
3052:           /* Change constraint numbering from field component to local dof number */
3053:           PetscSectionGetFieldConstraintIndices(section, p, f, &find);
3054:           for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
3055:           numConst += fcdof;
3056:           foff     += fdof;
3057:         }
3058:         if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
3059:       } else {
3060:         for (d = 0; d < cdof; ++d) indices[d] = d;
3061:       }
3062:       PetscSectionSetConstraintIndices(section, p, indices);
3063:     }
3064:   }
3065:   PetscFree(indices);
3066:   return(0);
3067: }

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

3074:   Not Collective

3076:   Input Parameters:
3077: + dm        - The DMPlex object
3078: . dim       - The spatial dimension of the problem
3079: . numFields - The number of fields in the problem
3080: . numComp   - An array of size numFields that holds the number of components for each field
3081: . numDof    - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
3082: . numBC     - The number of boundary conditions
3083: . bcField   - An array of size numBC giving the field number for each boundry condition
3084: . bcComps   - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
3085: . bcPoints  - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
3086: - perm      - Optional permutation of the chart, or NULL

3088:   Output Parameter:
3089: . section - The PetscSection object

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

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

3096:   Level: developer

3098:   Fortran Notes:
3099:   A Fortran 90 version is available as DMPlexCreateSectionF90()

3101: .keywords: mesh, elements
3102: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
3103: @*/
3104: PetscErrorCode DMPlexCreateSection(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], IS perm, PetscSection *section)
3105: {
3106:   PetscSection   aSec;

3110:   DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
3111:   DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
3112:   if (perm) {PetscSectionSetPermutation(*section, perm);}
3113:   PetscSectionSetUp(*section);
3114:   DMPlexGetAnchors(dm,&aSec,NULL);
3115:   if (numBC || aSec) {
3116:     DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
3117:     DMPlexCreateSectionBCIndices(dm, *section);
3118:   }
3119:   PetscSectionViewFromOptions(*section,NULL,"-section_view");
3120:   return(0);
3121: }

3125: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
3126: {
3127:   PetscSection   section;

3131:   DMClone(dm, cdm);
3132:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
3133:   DMSetDefaultSection(*cdm, section);
3134:   PetscSectionDestroy(&section);
3135:   return(0);
3136: }

3140: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
3141: {
3142:   DM_Plex *mesh = (DM_Plex*) dm->data;

3146:   if (section) *section = mesh->coneSection;
3147:   return(0);
3148: }

3152: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
3153: {
3154:   DM_Plex *mesh = (DM_Plex*) dm->data;

3158:   if (section) *section = mesh->supportSection;
3159:   return(0);
3160: }

3164: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
3165: {
3166:   DM_Plex *mesh = (DM_Plex*) dm->data;

3170:   if (cones) *cones = mesh->cones;
3171:   return(0);
3172: }

3176: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
3177: {
3178:   DM_Plex *mesh = (DM_Plex*) dm->data;

3182:   if (coneOrientations) *coneOrientations = mesh->coneOrientations;
3183:   return(0);
3184: }

3186: /******************************** FEM Support **********************************/

3190: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3191: {
3192:   PetscScalar    *array, *vArray;
3193:   const PetscInt *cone, *coneO;
3194:   PetscInt        pStart, pEnd, p, numPoints, size = 0, offset = 0;
3195:   PetscErrorCode  ierr;

3198:   PetscSectionGetChart(section, &pStart, &pEnd);
3199:   DMPlexGetConeSize(dm, point, &numPoints);
3200:   DMPlexGetCone(dm, point, &cone);
3201:   DMPlexGetConeOrientation(dm, point, &coneO);
3202:   if (!values || !*values) {
3203:     if ((point >= pStart) && (point < pEnd)) {
3204:       PetscInt dof;

3206:       PetscSectionGetDof(section, point, &dof);
3207:       size += dof;
3208:     }
3209:     for (p = 0; p < numPoints; ++p) {
3210:       const PetscInt cp = cone[p];
3211:       PetscInt       dof;

3213:       if ((cp < pStart) || (cp >= pEnd)) continue;
3214:       PetscSectionGetDof(section, cp, &dof);
3215:       size += dof;
3216:     }
3217:     if (!values) {
3218:       if (csize) *csize = size;
3219:       return(0);
3220:     }
3221:     DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
3222:   } else {
3223:     array = *values;
3224:   }
3225:   size = 0;
3226:   VecGetArray(v, &vArray);
3227:   if ((point >= pStart) && (point < pEnd)) {
3228:     PetscInt     dof, off, d;
3229:     PetscScalar *varr;

3231:     PetscSectionGetDof(section, point, &dof);
3232:     PetscSectionGetOffset(section, point, &off);
3233:     varr = &vArray[off];
3234:     for (d = 0; d < dof; ++d, ++offset) {
3235:       array[offset] = varr[d];
3236:     }
3237:     size += dof;
3238:   }
3239:   for (p = 0; p < numPoints; ++p) {
3240:     const PetscInt cp = cone[p];
3241:     PetscInt       o  = coneO[p];
3242:     PetscInt       dof, off, d;
3243:     PetscScalar   *varr;

3245:     if ((cp < pStart) || (cp >= pEnd)) continue;
3246:     PetscSectionGetDof(section, cp, &dof);
3247:     PetscSectionGetOffset(section, cp, &off);
3248:     varr = &vArray[off];
3249:     if (o >= 0) {
3250:       for (d = 0; d < dof; ++d, ++offset) {
3251:         array[offset] = varr[d];
3252:       }
3253:     } else {
3254:       for (d = dof-1; d >= 0; --d, ++offset) {
3255:         array[offset] = varr[d];
3256:       }
3257:     }
3258:     size += dof;
3259:   }
3260:   VecRestoreArray(v, &vArray);
3261:   if (!*values) {
3262:     if (csize) *csize = size;
3263:     *values = array;
3264:   } else {
3265:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3266:     *csize = size;
3267:   }
3268:   return(0);
3269: }

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

3279:   *size = 0;
3280:   for (p = 0; p < numPoints*2; p += 2) {
3281:     const PetscInt point = points[p];
3282:     const PetscInt o     = points[p+1];
3283:     PetscInt       dof, off, d;
3284:     const PetscScalar *varr;

3286:     PetscSectionGetDof(section, point, &dof);
3287:     PetscSectionGetOffset(section, point, &off);
3288:     varr = &vArray[off];
3289:     if (o >= 0) {
3290:       for (d = 0; d < dof; ++d, ++offset)    array[offset] = varr[d];
3291:     } else {
3292:       for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
3293:     }
3294:   }
3295:   *size = offset;
3296:   return(0);
3297: }

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

3307:   *size = 0;
3308:   for (f = 0; f < numFields; ++f) {
3309:     PetscInt fcomp, p;

3311:     PetscSectionGetFieldComponents(section, f, &fcomp);
3312:     for (p = 0; p < numPoints*2; p += 2) {
3313:       const PetscInt point = points[p];
3314:       const PetscInt o     = points[p+1];
3315:       PetscInt       fdof, foff, d, c;
3316:       const PetscScalar *varr;

3318:       PetscSectionGetFieldDof(section, point, f, &fdof);
3319:       PetscSectionGetFieldOffset(section, point, f, &foff);
3320:       varr = &vArray[foff];
3321:       if (o >= 0) {
3322:         for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
3323:       } else {
3324:         for (d = fdof/fcomp-1; d >= 0; --d) {
3325:           for (c = 0; c < fcomp; ++c, ++offset) {
3326:             array[offset] = varr[d*fcomp+c];
3327:           }
3328:         }
3329:       }
3330:     }
3331:   }
3332:   *size = offset;
3333:   return(0);
3334: }

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

3341:   Not collective

3343:   Input Parameters:
3344: + dm - The DM
3345: . section - The section describing the layout in v, or NULL to use the default section
3346: . v - The local vector
3347: - point - The sieve point in the DM

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

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

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

3359:   Level: intermediate

3361: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3362: @*/
3363: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3364: {
3365:   PetscSection    clSection;
3366:   IS              clPoints;
3367:   PetscScalar    *array, *vArray;
3368:   PetscInt       *points = NULL;
3369:   const PetscInt *clp;
3370:   PetscInt        depth, numFields, numPoints, size;
3371:   PetscErrorCode  ierr;

3375:   if (!section) {DMGetDefaultSection(dm, &section);}
3378:   DMPlexGetDepth(dm, &depth);
3379:   PetscSectionGetNumFields(section, &numFields);
3380:   if (depth == 1 && numFields < 2) {
3381:     DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
3382:     return(0);
3383:   }
3384:   /* Get points */
3385:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3386:   if (!clPoints) {
3387:     PetscInt pStart, pEnd, p, q;

3389:     PetscSectionGetChart(section, &pStart, &pEnd);
3390:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3391:     /* Compress out points not in the section */
3392:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3393:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3394:         points[q*2]   = points[p];
3395:         points[q*2+1] = points[p+1];
3396:         ++q;
3397:       }
3398:     }
3399:     numPoints = q;
3400:   } else {
3401:     PetscInt dof, off;

3403:     PetscSectionGetDof(clSection, point, &dof);
3404:     PetscSectionGetOffset(clSection, point, &off);
3405:     ISGetIndices(clPoints, &clp);
3406:     numPoints = dof/2;
3407:     points    = (PetscInt *) &clp[off];
3408:   }
3409:   /* Get array */
3410:   if (!values || !*values) {
3411:     PetscInt asize = 0, dof, p;

3413:     for (p = 0; p < numPoints*2; p += 2) {
3414:       PetscSectionGetDof(section, points[p], &dof);
3415:       asize += dof;
3416:     }
3417:     if (!values) {
3418:       if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3419:       else           {ISRestoreIndices(clPoints, &clp);}
3420:       if (csize) *csize = asize;
3421:       return(0);
3422:     }
3423:     DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
3424:   } else {
3425:     array = *values;
3426:   }
3427:   VecGetArray(v, &vArray);
3428:   /* Get values */
3429:   if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);}
3430:   else               {DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);}
3431:   /* Cleanup points */
3432:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3433:   else           {ISRestoreIndices(clPoints, &clp);}
3434:   /* Cleanup array */
3435:   VecRestoreArray(v, &vArray);
3436:   if (!*values) {
3437:     if (csize) *csize = size;
3438:     *values = array;
3439:   } else {
3440:     if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3441:     *csize = size;
3442:   }
3443:   return(0);
3444: }

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

3451:   Not collective

3453:   Input Parameters:
3454: + dm - The DM
3455: . section - The section describing the layout in v, or NULL to use the default section
3456: . v - The local vector
3457: . point - The sieve point in the DM
3458: . csize - The number of values in the closure, or NULL
3459: - values - The array of values, which is a borrowed array and should not be freed

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

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

3467:   Level: intermediate

3469: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3470: @*/
3471: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3472: {
3473:   PetscInt       size = 0;

3477:   /* Should work without recalculating size */
3478:   DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
3479:   return(0);
3480: }

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

3487: 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[])
3488: {
3489:   PetscInt        cdof;   /* The number of constraints on this point */
3490:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3491:   PetscScalar    *a;
3492:   PetscInt        off, cind = 0, k;
3493:   PetscErrorCode  ierr;

3496:   PetscSectionGetConstraintDof(section, point, &cdof);
3497:   PetscSectionGetOffset(section, point, &off);
3498:   a    = &array[off];
3499:   if (!cdof || setBC) {
3500:     if (orientation >= 0) {
3501:       for (k = 0; k < dof; ++k) {
3502:         fuse(&a[k], values[k]);
3503:       }
3504:     } else {
3505:       for (k = 0; k < dof; ++k) {
3506:         fuse(&a[k], values[dof-k-1]);
3507:       }
3508:     }
3509:   } else {
3510:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3511:     if (orientation >= 0) {
3512:       for (k = 0; k < dof; ++k) {
3513:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3514:         fuse(&a[k], values[k]);
3515:       }
3516:     } else {
3517:       for (k = 0; k < dof; ++k) {
3518:         if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3519:         fuse(&a[k], values[dof-k-1]);
3520:       }
3521:     }
3522:   }
3523:   return(0);
3524: }

3528: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3529: {
3530:   PetscInt        cdof;   /* The number of constraints on this point */
3531:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3532:   PetscScalar    *a;
3533:   PetscInt        off, cind = 0, k;
3534:   PetscErrorCode  ierr;

3537:   PetscSectionGetConstraintDof(section, point, &cdof);
3538:   PetscSectionGetOffset(section, point, &off);
3539:   a    = &array[off];
3540:   if (cdof) {
3541:     PetscSectionGetConstraintIndices(section, point, &cdofs);
3542:     if (orientation >= 0) {
3543:       for (k = 0; k < dof; ++k) {
3544:         if ((cind < cdof) && (k == cdofs[cind])) {
3545:           fuse(&a[k], values[k]);
3546:           ++cind;
3547:         }
3548:       }
3549:     } else {
3550:       for (k = 0; k < dof; ++k) {
3551:         if ((cind < cdof) && (k == cdofs[cind])) {
3552:           fuse(&a[k], values[dof-k-1]);
3553:           ++cind;
3554:         }
3555:       }
3556:     }
3557:   }
3558:   return(0);
3559: }

3563: 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[])
3564: {
3565:   PetscScalar    *a;
3566:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3567:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3568:   PetscInt        cind = 0, k, c;
3569:   PetscErrorCode  ierr;

3572:   PetscSectionGetFieldDof(section, point, f, &fdof);
3573:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3574:   PetscSectionGetFieldOffset(section, point, f, &foff);
3575:   a    = &array[foff];
3576:   if (!fcdof || setBC) {
3577:     if (o >= 0) {
3578:       for (k = 0; k < fdof; ++k) fuse(&a[k], values[foffset+k]);
3579:     } else {
3580:       for (k = fdof/fcomp-1; k >= 0; --k) {
3581:         for (c = 0; c < fcomp; ++c) {
3582:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3583:         }
3584:       }
3585:     }
3586:   } else {
3587:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3588:     if (o >= 0) {
3589:       for (k = 0; k < fdof; ++k) {
3590:         if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
3591:         fuse(&a[k], values[foffset+k]);
3592:       }
3593:     } else {
3594:       for (k = fdof/fcomp-1; k >= 0; --k) {
3595:         for (c = 0; c < fcomp; ++c) {
3596:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
3597:           fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3598:         }
3599:       }
3600:     }
3601:   }
3602:   *offset += fdof;
3603:   return(0);
3604: }

3608: 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[])
3609: {
3610:   PetscScalar    *a;
3611:   PetscInt        fdof, foff, fcdof, foffset = *offset;
3612:   const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3613:   PetscInt        cind = 0, k, c;
3614:   PetscErrorCode  ierr;

3617:   PetscSectionGetFieldDof(section, point, f, &fdof);
3618:   PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3619:   PetscSectionGetFieldOffset(section, point, f, &foff);
3620:   a    = &array[foff];
3621:   if (fcdof) {
3622:     PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3623:     if (o >= 0) {
3624:       for (k = 0; k < fdof; ++k) {
3625:         if ((cind < fcdof) && (k == fcdofs[cind])) {
3626:           fuse(&a[k], values[foffset+k]);
3627:           ++cind;
3628:         }
3629:       }
3630:     } else {
3631:       for (k = fdof/fcomp-1; k >= 0; --k) {
3632:         for (c = 0; c < fcomp; ++c) {
3633:           if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {
3634:             fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3635:             ++cind;
3636:           }
3637:         }
3638:       }
3639:     }
3640:   }
3641:   *offset += fdof;
3642:   return(0);
3643: }

3647: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3648: {
3649:   PetscScalar    *array;
3650:   const PetscInt *cone, *coneO;
3651:   PetscInt        pStart, pEnd, p, numPoints, off, dof;
3652:   PetscErrorCode  ierr;

3655:   PetscSectionGetChart(section, &pStart, &pEnd);
3656:   DMPlexGetConeSize(dm, point, &numPoints);
3657:   DMPlexGetCone(dm, point, &cone);
3658:   DMPlexGetConeOrientation(dm, point, &coneO);
3659:   VecGetArray(v, &array);
3660:   for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
3661:     const PetscInt cp = !p ? point : cone[p-1];
3662:     const PetscInt o  = !p ? 0     : coneO[p-1];

3664:     if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
3665:     PetscSectionGetDof(section, cp, &dof);
3666:     /* ADD_VALUES */
3667:     {
3668:       const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3669:       PetscScalar    *a;
3670:       PetscInt        cdof, coff, cind = 0, k;

3672:       PetscSectionGetConstraintDof(section, cp, &cdof);
3673:       PetscSectionGetOffset(section, cp, &coff);
3674:       a    = &array[coff];
3675:       if (!cdof) {
3676:         if (o >= 0) {
3677:           for (k = 0; k < dof; ++k) {
3678:             a[k] += values[off+k];
3679:           }
3680:         } else {
3681:           for (k = 0; k < dof; ++k) {
3682:             a[k] += values[off+dof-k-1];
3683:           }
3684:         }
3685:       } else {
3686:         PetscSectionGetConstraintIndices(section, cp, &cdofs);
3687:         if (o >= 0) {
3688:           for (k = 0; k < dof; ++k) {
3689:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3690:             a[k] += values[off+k];
3691:           }
3692:         } else {
3693:           for (k = 0; k < dof; ++k) {
3694:             if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3695:             a[k] += values[off+dof-k-1];
3696:           }
3697:         }
3698:       }
3699:     }
3700:   }
3701:   VecRestoreArray(v, &array);
3702:   return(0);
3703: }

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

3710:   Not collective

3712:   Input Parameters:
3713: + dm - The DM
3714: . section - The section describing the layout in v, or NULL to use the default section
3715: . v - The local vector
3716: . point - The sieve point in the DM
3717: . values - The array of values
3718: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

3723:   Level: intermediate

3725: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
3726: @*/
3727: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3728: {
3729:   PetscSection    clSection;
3730:   IS              clPoints;
3731:   PetscScalar    *array;
3732:   PetscInt       *points = NULL;
3733:   const PetscInt *clp;
3734:   PetscInt        depth, numFields, numPoints, p;
3735:   PetscErrorCode  ierr;

3739:   if (!section) {DMGetDefaultSection(dm, &section);}
3742:   DMPlexGetDepth(dm, &depth);
3743:   PetscSectionGetNumFields(section, &numFields);
3744:   if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
3745:     DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
3746:     return(0);
3747:   }
3748:   /* Get points */
3749:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3750:   if (!clPoints) {
3751:     PetscInt pStart, pEnd, q;

3753:     PetscSectionGetChart(section, &pStart, &pEnd);
3754:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3755:     /* Compress out points not in the section */
3756:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3757:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3758:         points[q*2]   = points[p];
3759:         points[q*2+1] = points[p+1];
3760:         ++q;
3761:       }
3762:     }
3763:     numPoints = q;
3764:   } else {
3765:     PetscInt dof, off;

3767:     PetscSectionGetDof(clSection, point, &dof);
3768:     PetscSectionGetOffset(clSection, point, &off);
3769:     ISGetIndices(clPoints, &clp);
3770:     numPoints = dof/2;
3771:     points    = (PetscInt *) &clp[off];
3772:   }
3773:   /* Get array */
3774:   VecGetArray(v, &array);
3775:   /* Get values */
3776:   if (numFields > 0) {
3777:     PetscInt offset = 0, fcomp, f;
3778:     for (f = 0; f < numFields; ++f) {
3779:       PetscSectionGetFieldComponents(section, f, &fcomp);
3780:       switch (mode) {
3781:       case INSERT_VALUES:
3782:         for (p = 0; p < numPoints*2; p += 2) {
3783:           const PetscInt point = points[p];
3784:           const PetscInt o     = points[p+1];
3785:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3786:         } break;
3787:       case INSERT_ALL_VALUES:
3788:         for (p = 0; p < numPoints*2; p += 2) {
3789:           const PetscInt point = points[p];
3790:           const PetscInt o     = points[p+1];
3791:           updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3792:         } break;
3793:       case INSERT_BC_VALUES:
3794:         for (p = 0; p < numPoints*2; p += 2) {
3795:           const PetscInt point = points[p];
3796:           const PetscInt o     = points[p+1];
3797:           updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3798:         } break;
3799:       case ADD_VALUES:
3800:         for (p = 0; p < numPoints*2; p += 2) {
3801:           const PetscInt point = points[p];
3802:           const PetscInt o     = points[p+1];
3803:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3804:         } break;
3805:       case ADD_ALL_VALUES:
3806:         for (p = 0; p < numPoints*2; p += 2) {
3807:           const PetscInt point = points[p];
3808:           const PetscInt o     = points[p+1];
3809:           updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3810:         } break;
3811:       default:
3812:         SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3813:       }
3814:     }
3815:   } else {
3816:     PetscInt dof, off;

3818:     switch (mode) {
3819:     case INSERT_VALUES:
3820:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3821:         PetscInt o = points[p+1];
3822:         PetscSectionGetDof(section, points[p], &dof);
3823:         updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
3824:       } break;
3825:     case INSERT_ALL_VALUES:
3826:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3827:         PetscInt o = points[p+1];
3828:         PetscSectionGetDof(section, points[p], &dof);
3829:         updatePoint_private(section, points[p], dof, insert, PETSC_TRUE,  o, &values[off], array);
3830:       } break;
3831:     case INSERT_BC_VALUES:
3832:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3833:         PetscInt o = points[p+1];
3834:         PetscSectionGetDof(section, points[p], &dof);
3835:         updatePointBC_private(section, points[p], dof, insert,  o, &values[off], array);
3836:       } break;
3837:     case ADD_VALUES:
3838:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3839:         PetscInt o = points[p+1];
3840:         PetscSectionGetDof(section, points[p], &dof);
3841:         updatePoint_private(section, points[p], dof, add,    PETSC_FALSE, o, &values[off], array);
3842:       } break;
3843:     case ADD_ALL_VALUES:
3844:       for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3845:         PetscInt o = points[p+1];
3846:         PetscSectionGetDof(section, points[p], &dof);
3847:         updatePoint_private(section, points[p], dof, add,    PETSC_TRUE,  o, &values[off], array);
3848:       } break;
3849:     default:
3850:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3851:     }
3852:   }
3853:   /* Cleanup points */
3854:   if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3855:   else           {ISRestoreIndices(clPoints, &clp);}
3856:   /* Cleanup array */
3857:   VecRestoreArray(v, &array);
3858:   return(0);
3859: }

3863: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
3864: {
3865:   PetscSection    clSection;
3866:   IS              clPoints;
3867:   PetscScalar    *array;
3868:   PetscInt       *points = NULL;
3869:   const PetscInt *clp;
3870:   PetscInt        numFields, numPoints, p;
3871:   PetscInt        offset = 0, fcomp, f;
3872:   PetscErrorCode  ierr;

3876:   if (!section) {DMGetDefaultSection(dm, &section);}
3879:   PetscSectionGetNumFields(section, &numFields);
3880:   /* Get points */
3881:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3882:   if (!clPoints) {
3883:     PetscInt pStart, pEnd, q;

3885:     PetscSectionGetChart(section, &pStart, &pEnd);
3886:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3887:     /* Compress out points not in the section */
3888:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
3889:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
3890:         points[q*2]   = points[p];
3891:         points[q*2+1] = points[p+1];
3892:         ++q;
3893:       }
3894:     }
3895:     numPoints = q;
3896:   } else {
3897:     PetscInt dof, off;

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

3963: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
3964: {
3965:   PetscMPIInt    rank;
3966:   PetscInt       i, j;

3970:   MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
3971:   PetscViewerASCIIPrintf(viewer, "[%D]mat for sieve point %D\n", rank, point);
3972:   for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
3973:   for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
3974:   numCIndices = numCIndices ? numCIndices : numRIndices;
3975:   for (i = 0; i < numRIndices; i++) {
3976:     PetscViewerASCIIPrintf(viewer, "[%d]", rank);
3977:     for (j = 0; j < numCIndices; j++) {
3978: #if defined(PETSC_USE_COMPLEX)
3979:       PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
3980: #else
3981:       PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
3982: #endif
3983:     }
3984:     PetscViewerASCIIPrintf(viewer, "\n");
3985:   }
3986:   return(0);
3987: }

3991: /* . off - The global offset of this point */
3992: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
3993: {
3994:   PetscInt        dof;    /* The number of unknowns on this point */
3995:   PetscInt        cdof;   /* The number of constraints on this point */
3996:   const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3997:   PetscInt        cind = 0, k;
3998:   PetscErrorCode  ierr;

4001:   PetscSectionGetDof(section, point, &dof);
4002:   PetscSectionGetConstraintDof(section, point, &cdof);
4003:   if (!cdof || setBC) {
4004:     if (orientation >= 0) {
4005:       for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
4006:     } else {
4007:       for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
4008:     }
4009:   } else {
4010:     PetscSectionGetConstraintIndices(section, point, &cdofs);
4011:     if (orientation >= 0) {
4012:       for (k = 0; k < dof; ++k) {
4013:         if ((cind < cdof) && (k == cdofs[cind])) {
4014:           /* Insert check for returning constrained indices */
4015:           indices[*loff+k] = -(off+k+1);
4016:           ++cind;
4017:         } else {
4018:           indices[*loff+k] = off+k-cind;
4019:         }
4020:       }
4021:     } else {
4022:       for (k = 0; k < dof; ++k) {
4023:         if ((cind < cdof) && (k == cdofs[cind])) {
4024:           /* Insert check for returning constrained indices */
4025:           indices[*loff+dof-k-1] = -(off+k+1);
4026:           ++cind;
4027:         } else {
4028:           indices[*loff+dof-k-1] = off+k-cind;
4029:         }
4030:       }
4031:     }
4032:   }
4033:   *loff += dof;
4034:   return(0);
4035: }

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

4046:   PetscSectionGetNumFields(section, &numFields);
4047:   for (f = 0, foff = 0; f < numFields; ++f) {
4048:     PetscInt        fdof, fcomp, cfdof;
4049:     const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4050:     PetscInt        cind = 0, k, c;

4052:     PetscSectionGetFieldComponents(section, f, &fcomp);
4053:     PetscSectionGetFieldDof(section, point, f, &fdof);
4054:     PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
4055:     if (!cfdof || setBC) {
4056:       if (orientation >= 0) {
4057:         for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
4058:       } else {
4059:         for (k = fdof/fcomp-1; k >= 0; --k) {
4060:           for (c = 0; c < fcomp; ++c) {
4061:             indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
4062:           }
4063:         }
4064:       }
4065:     } else {
4066:       PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
4067:       if (orientation >= 0) {
4068:         for (k = 0; k < fdof; ++k) {
4069:           if ((cind < cfdof) && (k == fcdofs[cind])) {
4070:             indices[foffs[f]+k] = -(off+foff+k+1);
4071:             ++cind;
4072:           } else {
4073:             indices[foffs[f]+k] = off+foff+k-cind;
4074:           }
4075:         }
4076:       } else {
4077:         for (k = fdof/fcomp-1; k >= 0; --k) {
4078:           for (c = 0; c < fcomp; ++c) {
4079:             if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
4080:               indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
4081:               ++cind;
4082:             } else {
4083:               indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
4084:             }
4085:           }
4086:         }
4087:       }
4088:     }
4089:     foff     += fdof - cfdof;
4090:     foffs[f] += fdof;
4091:   }
4092:   return(0);
4093: }

4097: static PetscErrorCode DMPlexAnchorsModifyMat(DM dm, PetscSection section, PetscInt numPoints, PetscInt numIndices, const PetscInt points[], const PetscScalar values[], PetscInt *outNumPoints, PetscInt *outNumIndices, PetscInt *outPoints[], PetscScalar *outValues[], PetscInt offsets[])
4098: {
4099:   Mat             cMat;
4100:   PetscSection    aSec, cSec;
4101:   IS              aIS;
4102:   PetscInt        aStart = -1, aEnd = -1;
4103:   const PetscInt  *anchors;
4104:   PetscInt        numFields, f, p, q, newP = 0;
4105:   PetscInt        newNumPoints = 0, newNumIndices = 0;
4106:   PetscInt        *newPoints, *indices, *newIndices;
4107:   PetscInt        maxAnchor, maxDof;
4108:   PetscInt        newOffsets[32];
4109:   PetscInt        *pointMatOffsets[32];
4110:   PetscInt        *newPointOffsets[32];
4111:   PetscScalar     *pointMat[32];
4112:   PetscScalar     *newValues,*tmpValues;
4113:   PetscBool       anyConstrained = PETSC_FALSE;
4114:   PetscErrorCode  ierr;

4119:   PetscSectionGetNumFields(section, &numFields);

4121:   DMPlexGetAnchors(dm,&aSec,&aIS);
4122:   /* if there are point-to-point constraints */
4123:   if (aSec) {
4124:     PetscMemzero(newOffsets, 32 * sizeof(PetscInt));
4125:     ISGetIndices(aIS,&anchors);
4126:     PetscSectionGetChart(aSec,&aStart,&aEnd);
4127:     /* figure out how many points are going to be in the new element matrix
4128:      * (we allow double counting, because it's all just going to be summed
4129:      * into the global matrix anyway) */
4130:     for (p = 0; p < 2*numPoints; p+=2) {
4131:       PetscInt b    = points[p];
4132:       PetscInt bDof = 0;

4134:       if (b >= aStart && b < aEnd) {
4135:         PetscSectionGetDof(aSec,b,&bDof);
4136:       }
4137:       if (bDof) {
4138:         /* this point is constrained */
4139:         /* it is going to be replaced by its anchors */
4140:         PetscInt bOff, q;

4142:         anyConstrained = PETSC_TRUE;
4143:         newNumPoints  += bDof;
4144:         PetscSectionGetOffset(aSec,b,&bOff);
4145:         for (q = 0; q < bDof; q++) {
4146:           PetscInt a = anchors[bOff + q];
4147:           PetscInt aDof;

4149:           PetscSectionGetDof(section,a,&aDof);
4150:           newNumIndices += aDof;
4151:           for (f = 0; f < numFields; ++f) {
4152:             PetscInt fDof;

4154:             PetscSectionGetFieldDof(section, a, f, &fDof);
4155:             newOffsets[f+1] += fDof;
4156:           }
4157:         }
4158:       }
4159:       else {
4160:         /* this point is not constrained */
4161:         newNumPoints++;
4162:         PetscSectionGetDof(section,b,&bDof);
4163:         newNumIndices += bDof;
4164:         for (f = 0; f < numFields; ++f) {
4165:           PetscInt fDof;

4167:           PetscSectionGetFieldDof(section, b, f, &fDof);
4168:           newOffsets[f+1] += fDof;
4169:         }
4170:       }
4171:     }
4172:   }
4173:   if (!anyConstrained) {
4174:     *outNumPoints  = 0;
4175:     *outNumIndices = 0;
4176:     *outPoints     = NULL;
4177:     *outValues     = NULL;
4178:     if (aSec) {
4179:       ISRestoreIndices(aIS,&anchors);
4180:     }
4181:     return(0);
4182:   }

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

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

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

4190:   /* output arrays */
4191:   DMGetWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4192:   DMGetWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);

4194:   /* workspaces */
4195:   DMGetWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4196:   if (numFields) {
4197:     for (f = 0; f < numFields; f++) {
4198:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4199:       DMGetWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4200:     }
4201:   }
4202:   else {
4203:     DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4204:     DMGetWorkArray(dm,numPoints,PETSC_INT,&newPointOffsets[0]);
4205:   }

4207:   /* get workspaces for the point-to-point matrices */
4208:   if (numFields) {
4209:     for (p = 0; p < numPoints; p++) {
4210:       PetscInt b    = points[2*p];
4211:       PetscInt bDof = 0;

4213:       if (b >= aStart && b < aEnd) {
4214:         PetscSectionGetDof(aSec, b, &bDof);
4215:       }
4216:       if (bDof) {
4217:         for (f = 0; f < numFields; f++) {
4218:           PetscInt fDof, q, bOff, allFDof = 0;

4220:           PetscSectionGetFieldDof(section, b, f, &fDof);
4221:           PetscSectionGetOffset(aSec, b, &bOff);
4222:           for (q = 0; q < bDof; q++) {
4223:             PetscInt a = anchors[bOff + q];
4224:             PetscInt aFDof;

4226:             PetscSectionGetFieldDof(section, a, f, &aFDof);
4227:             allFDof += aFDof;
4228:           }
4229:           newPointOffsets[f][p+1] = allFDof;
4230:           pointMatOffsets[f][p+1] = fDof * allFDof;
4231:         }
4232:       }
4233:       else {
4234:         for (f = 0; f < numFields; f++) {
4235:           PetscInt fDof;

4237:           PetscSectionGetFieldDof(section, b, f, &fDof);
4238:           newPointOffsets[f][p+1] = fDof;
4239:           pointMatOffsets[f][p+1] = 0;
4240:         }
4241:       }
4242:     }
4243:     for (f = 0; f < numFields; f++) {
4244:       newPointOffsets[f][0] = 0;
4245:       pointMatOffsets[f][0] = 0;
4246:       for (p = 0; p < numPoints; p++) {
4247:         newPointOffsets[f][p+1] += newPointOffsets[f][p];
4248:         pointMatOffsets[f][p+1] += pointMatOffsets[f][p];
4249:       }
4250:       DMGetWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4251:     }
4252:   }
4253:   else {
4254:     for (p = 0; p < numPoints; p++) {
4255:       PetscInt b    = points[2*p];
4256:       PetscInt bDof = 0;

4258:       if (b >= aStart && b < aEnd) {
4259:         PetscSectionGetDof(aSec, b, &bDof);
4260:       }
4261:       if (bDof) {
4262:         PetscInt dof, bOff, q, allDof = 0;

4264:         PetscSectionGetDof(section, b, &dof);
4265:         PetscSectionGetOffset(aSec, b, &bOff);
4266:         for (q = 0; q < bDof; q++) {
4267:           PetscInt a = anchors[bOff + q], aDof;

4269:           PetscSectionGetDof(section, a, &aDof);
4270:           allDof += aDof;
4271:         }
4272:         newPointOffsets[0][p+1] = allDof;
4273:         pointMatOffsets[0][p+1] = dof * allDof;
4274:       }
4275:       else {
4276:         PetscInt dof;

4278:         PetscSectionGetDof(section, b, &dof);
4279:         newPointOffsets[0][p+1] = dof;
4280:         pointMatOffsets[0][p+1] = 0;
4281:       }
4282:     }
4283:     newPointOffsets[0][0] = 0;
4284:     pointMatOffsets[0][0] = 0;
4285:     for (p = 0; p < numPoints; p++) {
4286:       newPointOffsets[0][p+1] += newPointOffsets[0][p];
4287:       pointMatOffsets[0][p+1] += pointMatOffsets[0][p];
4288:     }
4289:     DMGetWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4290:   }

4292:   /* get the point-to-point matrices; construct newPoints */
4293:   PetscSectionGetMaxDof(aSec, &maxAnchor);
4294:   PetscSectionGetMaxDof(section, &maxDof);
4295:   DMGetWorkArray(dm,maxDof,PETSC_INT,&indices);
4296:   DMGetWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4297:   if (numFields) {
4298:     for (p = 0, newP = 0; p < numPoints; p++) {
4299:       PetscInt b    = points[2*p];
4300:       PetscInt o    = points[2*p+1];
4301:       PetscInt bDof = 0;

4303:       if (b >= aStart && b < aEnd) {
4304:         PetscSectionGetDof(aSec, b, &bDof);
4305:       }
4306:       if (bDof) {
4307:         PetscInt fStart[32], fEnd[32], fAnchorStart[32], fAnchorEnd[32], bOff, q;

4309:         fStart[0] = 0;
4310:         fEnd[0]   = 0;
4311:         for (f = 0; f < numFields; f++) {
4312:           PetscInt fDof;

4314:           PetscSectionGetFieldDof(cSec, b, f, &fDof);
4315:           fStart[f+1] = fStart[f] + fDof;
4316:           fEnd[f+1]   = fStart[f+1];
4317:         }
4318:         PetscSectionGetOffset(cSec, b, &bOff);
4319:         indicesPointFields_private(cSec, b, bOff, fEnd, PETSC_TRUE, o, indices);

4321:         fAnchorStart[0] = 0;
4322:         fAnchorEnd[0]   = 0;
4323:         for (f = 0; f < numFields; f++) {
4324:           PetscInt fDof = newPointOffsets[f][p + 1] - newPointOffsets[f][p];

4326:           fAnchorStart[f+1] = fAnchorStart[f] + fDof;
4327:           fAnchorEnd[f+1]   = fAnchorStart[f + 1];
4328:         }
4329:         PetscSectionGetOffset (aSec, b, &bOff);
4330:         for (q = 0; q < bDof; q++) {
4331:           PetscInt a = anchors[bOff + q], aOff;

4333:           /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4334:           newPoints[2*(newP + q)]     = a;
4335:           newPoints[2*(newP + q) + 1] = 0;
4336:           PetscSectionGetOffset(section, a, &aOff);
4337:           indicesPointFields_private(section, a, aOff, fAnchorEnd, PETSC_TRUE, 0, newIndices);
4338:         }
4339:         newP += bDof;

4341:         /* get the point-to-point submatrix */
4342:         for (f = 0; f < numFields; f++) {
4343:           MatGetValues(cMat,fEnd[f]-fStart[f],indices + fStart[f],fAnchorEnd[f] - fAnchorStart[f],newIndices + fAnchorStart[f],pointMat[f] + pointMatOffsets[f][p]);
4344:         }
4345:       }
4346:       else {
4347:         newPoints[2 * newP]     = b;
4348:         newPoints[2 * newP + 1] = o;
4349:         newP++;
4350:       }
4351:     }
4352:   } else {
4353:     for (p = 0; p < numPoints; p++) {
4354:       PetscInt b    = points[2*p];
4355:       PetscInt o    = points[2*p+1];
4356:       PetscInt bDof = 0;

4358:       if (b >= aStart && b < aEnd) {
4359:         PetscSectionGetDof(aSec, b, &bDof);
4360:       }
4361:       if (bDof) {
4362:         PetscInt bEnd = 0, bAnchorEnd = 0, bOff;

4364:         PetscSectionGetOffset(cSec, b, &bOff);
4365:         indicesPoint_private(cSec, b, bOff, &bEnd, PETSC_TRUE, o, indices);

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

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

4373:           newPoints[2*(newP + q)]     = a;
4374:           newPoints[2*(newP + q) + 1] = 0;
4375:           PetscSectionGetOffset(section, a, &aOff);
4376:           indicesPoint_private(section, a, aOff, &bAnchorEnd, PETSC_TRUE, 0, newIndices);
4377:         }
4378:         newP += bDof;

4380:         /* get the point-to-point submatrix */
4381:         MatGetValues(cMat,bEnd,indices,bAnchorEnd,newIndices,pointMat[0] + pointMatOffsets[0][p]);
4382:       }
4383:       else {
4384:         newPoints[2 * newP]     = b;
4385:         newPoints[2 * newP + 1] = o;
4386:         newP++;
4387:       }
4388:     }
4389:   }

4391:   PetscMemzero(tmpValues,newNumIndices*numIndices*sizeof(*tmpValues));
4392:   /* multiply constraints on the right */
4393:   if (numFields) {
4394:     for (f = 0; f < numFields; f++) {
4395:       PetscInt oldOff = offsets[f];

4397:       for (p = 0; p < numPoints; p++) {
4398:         PetscInt cStart = newPointOffsets[f][p];
4399:         PetscInt b      = points[2 * p];
4400:         PetscInt c, r, k;
4401:         PetscInt dof;

4403:         PetscSectionGetFieldDof(section,b,f,&dof);
4404:         if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4405:           PetscInt nCols         = newPointOffsets[f][p+1]-cStart;
4406:           const PetscScalar *mat = pointMat[f] + pointMatOffsets[f][p];

4408:           for (r = 0; r < numIndices; r++) {
4409:             for (c = 0; c < nCols; c++) {
4410:               for (k = 0; k < dof; k++) {
4411:                 tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4412:               }
4413:             }
4414:           }
4415:         }
4416:         else {
4417:           /* copy this column as is */
4418:           for (r = 0; r < numIndices; r++) {
4419:             for (c = 0; c < dof; c++) {
4420:               tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4421:             }
4422:           }
4423:         }
4424:         oldOff += dof;
4425:       }
4426:     }
4427:   }
4428:   else {
4429:     PetscInt oldOff = 0;
4430:     for (p = 0; p < numPoints; p++) {
4431:       PetscInt cStart = newPointOffsets[0][p];
4432:       PetscInt b      = points[2 * p];
4433:       PetscInt c, r, k;
4434:       PetscInt dof;

4436:       PetscSectionGetDof(section,b,&dof);
4437:       if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4438:         PetscInt nCols         = newPointOffsets[0][p+1]-cStart;
4439:         const PetscScalar *mat = pointMat[0] + pointMatOffsets[0][p];

4441:         for (r = 0; r < numIndices; r++) {
4442:           for (c = 0; c < nCols; c++) {
4443:             for (k = 0; k < dof; k++) {
4444:               tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4445:             }
4446:           }
4447:         }
4448:       }
4449:       else {
4450:         /* copy this column as is */
4451:         for (r = 0; r < numIndices; r++) {
4452:           for (c = 0; c < dof; c++) {
4453:             tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4454:           }
4455:         }
4456:       }
4457:       oldOff += dof;
4458:     }
4459:   }

4461:   PetscMemzero(newValues,newNumIndices*newNumIndices*sizeof(*newValues));
4462:   /* multiply constraints transpose on the left */
4463:   if (numFields) {
4464:     for (f = 0; f < numFields; f++) {
4465:       PetscInt oldOff = offsets[f];

4467:       for (p = 0; p < numPoints; p++) {
4468:         PetscInt rStart = newPointOffsets[f][p];
4469:         PetscInt b      = points[2 * p];
4470:         PetscInt c, r, k;
4471:         PetscInt dof;

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

4478:           for (r = 0; r < nRows; r++) {
4479:             for (c = 0; c < newNumIndices; c++) {
4480:               for (k = 0; k < dof; k++) {
4481:                 newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4482:               }
4483:             }
4484:           }
4485:         }
4486:         else {
4487:           /* copy this row as is */
4488:           for (r = 0; r < dof; r++) {
4489:             for (c = 0; c < newNumIndices; c++) {
4490:               newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4491:             }
4492:           }
4493:         }
4494:         oldOff += dof;
4495:       }
4496:     }
4497:   }
4498:   else {
4499:     PetscInt oldOff = 0;

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

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

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

4532:   /* clean up */
4533:   DMRestoreWorkArray(dm,maxDof,PETSC_INT,&indices);
4534:   DMRestoreWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4535:   if (numFields) {
4536:     for (f = 0; f < numFields; f++) {
4537:       DMRestoreWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4538:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4539:       DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4540:     }
4541:   }
4542:   else {
4543:     DMRestoreWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4544:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4545:     DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[0]);
4546:   }
4547:   DMRestoreWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4548:   ISRestoreIndices(aIS,&anchors);

4550:   /* output */
4551:   *outNumPoints  = newNumPoints;
4552:   *outNumIndices = newNumIndices;
4553:   *outPoints     = newPoints;
4554:   *outValues     = newValues;
4555:   for (f = 0; f < numFields; f++) {
4556:     offsets[f] = newOffsets[f];
4557:   }
4558:   return(0);
4559: }

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

4566:   Not collective

4568:   Input Parameters:
4569: + dm - The DM
4570: . section - The section describing the layout in v, or NULL to use the default section
4571: . globalSection - The section describing the layout in v, or NULL to use the default global section
4572: . A - The matrix
4573: . point - The sieve point in the DM
4574: . values - The array of values
4575: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions

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

4580:   Level: intermediate

4582: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
4583: @*/
4584: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4585: {
4586:   DM_Plex        *mesh   = (DM_Plex*) dm->data;
4587:   PetscSection    clSection;
4588:   IS              clPoints;
4589:   PetscInt       *points = NULL, *newPoints;
4590:   const PetscInt *clp;
4591:   PetscInt       *indices;
4592:   PetscInt        offsets[32];
4593:   PetscInt        numFields, numPoints, newNumPoints, numIndices, newNumIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
4594:   PetscScalar    *newValues;
4595:   PetscErrorCode  ierr;

4599:   if (!section) {DMGetDefaultSection(dm, &section);}
4601:   if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
4604:   PetscSectionGetNumFields(section, &numFields);
4605:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4606:   PetscMemzero(offsets, 32 * sizeof(PetscInt));
4607:   PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4608:   if (!clPoints) {
4609:     DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4610:     /* Compress out points not in the section */
4611:     PetscSectionGetChart(section, &pStart, &pEnd);
4612:     for (p = 0, q = 0; p < numPoints*2; p += 2) {
4613:       if ((points[p] >= pStart) && (points[p] < pEnd)) {
4614:         points[q*2]   = points[p];
4615:         points[q*2+1] = points[p+1];
4616:         ++q;
4617:       }
4618:     }
4619:     numPoints = q;
4620:   } else {
4621:     PetscInt dof, off;

4623:     PetscSectionGetDof(clSection, point, &dof);
4624:     numPoints = dof/2;
4625:     PetscSectionGetOffset(clSection, point, &off);
4626:     ISGetIndices(clPoints, &clp);
4627:     points = (PetscInt *) &clp[off];
4628:   }
4629:   for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
4630:     PetscInt fdof;

4632:     PetscSectionGetDof(section, points[p], &dof);
4633:     for (f = 0; f < numFields; ++f) {
4634:       PetscSectionGetFieldDof(section, points[p], f, &fdof);
4635:       offsets[f+1] += fdof;
4636:     }
4637:     numIndices += dof;
4638:   }
4639:   for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];

4641:   if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
4642:   DMPlexAnchorsModifyMat(dm,section,numPoints,numIndices,points,values,&newNumPoints,&newNumIndices,&newPoints,&newValues,offsets);
4643:   if (newNumPoints) {
4644:     if (!clPoints) {
4645:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4646:     } else {
4647:       ISRestoreIndices(clPoints, &clp);
4648:     }
4649:     numPoints  = newNumPoints;
4650:     numIndices = newNumIndices;
4651:     points     = newPoints;
4652:     values     = newValues;
4653:   }
4654:   DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
4655:   if (numFields) {
4656:     for (p = 0; p < numPoints*2; p += 2) {
4657:       PetscInt o = points[p+1];
4658:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4659:       indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
4660:     }
4661:   } else {
4662:     for (p = 0, off = 0; p < numPoints*2; p += 2) {
4663:       PetscInt o = points[p+1];
4664:       PetscSectionGetOffset(globalSection, points[p], &globalOff);
4665:       indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
4666:     }
4667:   }
4668:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
4669:   MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
4670:   if (mesh->printFEM > 1) {
4671:     PetscInt i;
4672:     PetscPrintf(PETSC_COMM_SELF, "  Indices:");
4673:     for (i = 0; i < numIndices; ++i) {PetscPrintf(PETSC_COMM_SELF, " %d", indices[i]);}
4674:     PetscPrintf(PETSC_COMM_SELF, "\n");
4675:   }
4676:   if (ierr) {
4677:     PetscMPIInt    rank;
4678:     PetscErrorCode ierr2;

4680:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4681:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4682:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
4683:     ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
4684: 
4685:   }
4686:   if (newNumPoints) {
4687:     DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4688:     DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4689:   }
4690:   else {
4691:     if (!clPoints) {
4692:       DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4693:     } else {
4694:       ISRestoreIndices(clPoints, &clp);
4695:     }
4696:   }
4697:   DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
4698:   return(0);
4699: }

4703: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4704: {
4705:   DM_Plex        *mesh   = (DM_Plex*) dmf->data;
4706:   PetscInt       *fpoints = NULL, *ftotpoints = NULL;
4707:   PetscInt       *cpoints = NULL;
4708:   PetscInt       *findices, *cindices;
4709:   PetscInt        foffsets[32], coffsets[32];
4710:   CellRefiner     cellRefiner;
4711:   PetscInt        numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4712:   PetscErrorCode  ierr;

4717:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4719:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
4721:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4723:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4726:   PetscSectionGetNumFields(fsection, &numFields);
4727:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4728:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4729:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4730:   /* Column indices */
4731:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4732:   maxFPoints = numCPoints;
4733:   /* Compress out points not in the section */
4734:   /*   TODO: Squeeze out points with 0 dof as well */
4735:   PetscSectionGetChart(csection, &pStart, &pEnd);
4736:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4737:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4738:       cpoints[q*2]   = cpoints[p];
4739:       cpoints[q*2+1] = cpoints[p+1];
4740:       ++q;
4741:     }
4742:   }
4743:   numCPoints = q;
4744:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4745:     PetscInt fdof;

4747:     PetscSectionGetDof(csection, cpoints[p], &dof);
4748:     if (!dof) continue;
4749:     for (f = 0; f < numFields; ++f) {
4750:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4751:       coffsets[f+1] += fdof;
4752:     }
4753:     numCIndices += dof;
4754:   }
4755:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4756:   /* Row indices */
4757:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4758:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4759:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4760:   for (r = 0, q = 0; r < numSubcells; ++r) {
4761:     /* TODO Map from coarse to fine cells */
4762:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4763:     /* Compress out points not in the section */
4764:     PetscSectionGetChart(fsection, &pStart, &pEnd);
4765:     for (p = 0; p < numFPoints*2; p += 2) {
4766:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4767:         PetscSectionGetDof(fsection, fpoints[p], &dof);
4768:         if (!dof) continue;
4769:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4770:         if (s < q) continue;
4771:         ftotpoints[q*2]   = fpoints[p];
4772:         ftotpoints[q*2+1] = fpoints[p+1];
4773:         ++q;
4774:       }
4775:     }
4776:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4777:   }
4778:   numFPoints = q;
4779:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4780:     PetscInt fdof;

4782:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4783:     if (!dof) continue;
4784:     for (f = 0; f < numFields; ++f) {
4785:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4786:       foffsets[f+1] += fdof;
4787:     }
4788:     numFIndices += dof;
4789:   }
4790:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

4792:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4793:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4794:   DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4795:   DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4796:   if (numFields) {
4797:     for (p = 0; p < numFPoints*2; p += 2) {
4798:       PetscInt o = ftotpoints[p+1];
4799:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4800:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4801:     }
4802:     for (p = 0; p < numCPoints*2; p += 2) {
4803:       PetscInt o = cpoints[p+1];
4804:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4805:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4806:     }
4807:   } else {
4808:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4809:       PetscInt o = ftotpoints[p+1];
4810:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4811:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4812:     }
4813:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4814:       PetscInt o = cpoints[p+1];
4815:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4816:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4817:     }
4818:   }
4819:   if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
4820:   MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
4821:   if (ierr) {
4822:     PetscMPIInt    rank;
4823:     PetscErrorCode ierr2;

4825:     ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4826:     ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4827:     ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
4828:     ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
4829:     ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
4830: 
4831:   }
4832:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4833:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4834:   DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4835:   DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4836:   return(0);
4837: }

4841: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
4842: {
4843:   PetscInt      *fpoints = NULL, *ftotpoints = NULL;
4844:   PetscInt      *cpoints = NULL;
4845:   PetscInt       foffsets[32], coffsets[32];
4846:   CellRefiner    cellRefiner;
4847:   PetscInt       numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;

4853:   if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4855:   if (!csection) {DMGetDefaultSection(dmc, &csection);}
4857:   if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4859:   if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4861:   PetscSectionGetNumFields(fsection, &numFields);
4862:   if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4863:   PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4864:   PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4865:   /* Column indices */
4866:   DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4867:   maxFPoints = numCPoints;
4868:   /* Compress out points not in the section */
4869:   /*   TODO: Squeeze out points with 0 dof as well */
4870:   PetscSectionGetChart(csection, &pStart, &pEnd);
4871:   for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4872:     if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4873:       cpoints[q*2]   = cpoints[p];
4874:       cpoints[q*2+1] = cpoints[p+1];
4875:       ++q;
4876:     }
4877:   }
4878:   numCPoints = q;
4879:   for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4880:     PetscInt fdof;

4882:     PetscSectionGetDof(csection, cpoints[p], &dof);
4883:     if (!dof) continue;
4884:     for (f = 0; f < numFields; ++f) {
4885:       PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4886:       coffsets[f+1] += fdof;
4887:     }
4888:     numCIndices += dof;
4889:   }
4890:   for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4891:   /* Row indices */
4892:   DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4893:   CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4894:   DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4895:   for (r = 0, q = 0; r < numSubcells; ++r) {
4896:     /* TODO Map from coarse to fine cells */
4897:     DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4898:     /* Compress out points not in the section */
4899:     PetscSectionGetChart(fsection, &pStart, &pEnd);
4900:     for (p = 0; p < numFPoints*2; p += 2) {
4901:       if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4902:         PetscSectionGetDof(fsection, fpoints[p], &dof);
4903:         if (!dof) continue;
4904:         for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4905:         if (s < q) continue;
4906:         ftotpoints[q*2]   = fpoints[p];
4907:         ftotpoints[q*2+1] = fpoints[p+1];
4908:         ++q;
4909:       }
4910:     }
4911:     DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4912:   }
4913:   numFPoints = q;
4914:   for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4915:     PetscInt fdof;

4917:     PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4918:     if (!dof) continue;
4919:     for (f = 0; f < numFields; ++f) {
4920:       PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4921:       foffsets[f+1] += fdof;
4922:     }
4923:     numFIndices += dof;
4924:   }
4925:   for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];

4927:   if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4928:   if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4929:   if (numFields) {
4930:     for (p = 0; p < numFPoints*2; p += 2) {
4931:       PetscInt o = ftotpoints[p+1];
4932:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4933:       indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4934:     }
4935:     for (p = 0; p < numCPoints*2; p += 2) {
4936:       PetscInt o = cpoints[p+1];
4937:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4938:       indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4939:     }
4940:   } else {
4941:     for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4942:       PetscInt o = ftotpoints[p+1];
4943:       PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4944:       indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4945:     }
4946:     for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4947:       PetscInt o = cpoints[p+1];
4948:       PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4949:       indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4950:     }
4951:   }
4952:   DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4953:   DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4954:   return(0);
4955: }

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

4962:   Input Parameter:
4963: . dm - The DMPlex object

4965:   Output Parameters:
4966: + cMax - The first hybrid cell
4967: . fMax - The first hybrid face
4968: . eMax - The first hybrid edge
4969: - vMax - The first hybrid vertex

4971:   Level: developer

4973: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
4974: @*/
4975: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
4976: {
4977:   DM_Plex       *mesh = (DM_Plex*) dm->data;
4978:   PetscInt       dim;

4983:   DMGetDimension(dm, &dim);
4984:   if (cMax) *cMax = mesh->hybridPointMax[dim];
4985:   if (fMax) *fMax = mesh->hybridPointMax[dim-1];
4986:   if (eMax) *eMax = mesh->hybridPointMax[1];
4987:   if (vMax) *vMax = mesh->hybridPointMax[0];
4988:   return(0);
4989: }

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

4996:   Input Parameters:
4997: . dm   - The DMPlex object
4998: . cMax - The first hybrid cell
4999: . fMax - The first hybrid face
5000: . eMax - The first hybrid edge
5001: - vMax - The first hybrid vertex

5003:   Level: developer

5005: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5006: @*/
5007: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5008: {
5009:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5010:   PetscInt       dim;

5015:   DMGetDimension(dm, &dim);
5016:   if (cMax >= 0) mesh->hybridPointMax[dim]   = cMax;
5017:   if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5018:   if (eMax >= 0) mesh->hybridPointMax[1]     = eMax;
5019:   if (vMax >= 0) mesh->hybridPointMax[0]     = vMax;
5020:   return(0);
5021: }

5025: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5026: {
5027:   DM_Plex *mesh = (DM_Plex*) dm->data;

5032:   *cellHeight = mesh->vtkCellHeight;
5033:   return(0);
5034: }

5038: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5039: {
5040:   DM_Plex *mesh = (DM_Plex*) dm->data;

5044:   mesh->vtkCellHeight = cellHeight;
5045:   return(0);
5046: }

5050: /* We can easily have a form that takes an IS instead */
5051: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5052: {
5053:   PetscSection   section, globalSection;
5054:   PetscInt      *numbers, p;

5058:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
5059:   PetscSectionSetChart(section, pStart, pEnd);
5060:   for (p = pStart; p < pEnd; ++p) {
5061:     PetscSectionSetDof(section, p, 1);
5062:   }
5063:   PetscSectionSetUp(section);
5064:   PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, &globalSection);
5065:   PetscMalloc1(pEnd - pStart, &numbers);
5066:   for (p = pStart; p < pEnd; ++p) {
5067:     PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5068:     if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5069:     else                       numbers[p-pStart] += shift;
5070:   }
5071:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5072:   if (globalSize) {
5073:     PetscLayout layout;
5074:     PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5075:     PetscLayoutGetSize(layout, globalSize);
5076:     PetscLayoutDestroy(&layout);
5077:   }
5078:   PetscSectionDestroy(&section);
5079:   PetscSectionDestroy(&globalSection);
5080:   return(0);
5081: }

5085: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5086: {
5087:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5088:   PetscInt       cellHeight, cStart, cEnd, cMax;

5093:   if (!mesh->globalCellNumbers) {
5094:     DMPlexGetVTKCellHeight(dm, &cellHeight);
5095:     DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5096:     DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5097:     if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
5098:     DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, &mesh->globalCellNumbers);
5099:   }
5100:   *globalCellNumbers = mesh->globalCellNumbers;
5101:   return(0);
5102: }

5106: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5107: {
5108:   DM_Plex       *mesh = (DM_Plex*) dm->data;
5109:   PetscInt       vStart, vEnd, vMax;

5114:   if (!mesh->globalVertexNumbers) {
5115:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5116:     DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5117:     if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
5118:     DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, &mesh->globalVertexNumbers);
5119:   }
5120:   *globalVertexNumbers = mesh->globalVertexNumbers;
5121:   return(0);
5122: }

5126: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5127: {
5128:   IS             nums[4];
5129:   PetscInt       depths[4];
5130:   PetscInt       depth, d, shift = 0;

5135:   DMPlexGetDepth(dm, &depth);
5136:   /* For unstratified meshes use dim instead of depth */
5137:   if (depth < 0) {DMGetDimension(dm, &depth);}
5138:   depths[0] = depth; depths[1] = 0;
5139:   for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5140:   for (d = 0; d <= depth; ++d) {
5141:     PetscInt pStart, pEnd, gsize;

5143:     DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5144:     DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5145:     shift += gsize;
5146:   }
5147:   ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5148:   for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5149:   return(0);
5150: }


5155: /*@C
5156:   PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
5157:   the local section and an SF describing the section point overlap.

5159:   Input Parameters:
5160:   + s - The PetscSection for the local field layout
5161:   . sf - The SF describing parallel layout of the section points
5162:   . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
5163:   . label - The label specifying the points
5164:   - labelValue - The label stratum specifying the points

5166:   Output Parameter:
5167:   . gsection - The PetscSection for the global field layout

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

5171:   Level: developer

5173: .seealso: PetscSectionCreate()
5174: @*/
5175: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
5176: {
5177:   PetscInt      *neg = NULL, *tmpOff = NULL;
5178:   PetscInt       pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;

5182:   PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
5183:   PetscSectionGetChart(s, &pStart, &pEnd);
5184:   PetscSectionSetChart(*gsection, pStart, pEnd);
5185:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
5186:   if (nroots >= 0) {
5187:     if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
5188:     PetscCalloc1(nroots, &neg);
5189:     if (nroots > pEnd-pStart) {
5190:       PetscCalloc1(nroots, &tmpOff);
5191:     } else {
5192:       tmpOff = &(*gsection)->atlasDof[-pStart];
5193:     }
5194:   }
5195:   /* Mark ghost points with negative dof */
5196:   for (p = pStart; p < pEnd; ++p) {
5197:     PetscInt value;

5199:     DMLabelGetValue(label, p, &value);
5200:     if (value != labelValue) continue;
5201:     PetscSectionGetDof(s, p, &dof);
5202:     PetscSectionSetDof(*gsection, p, dof);
5203:     PetscSectionGetConstraintDof(s, p, &cdof);
5204:     if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
5205:     if (neg) neg[p] = -(dof+1);
5206:   }
5207:   PetscSectionSetUpBC(*gsection);
5208:   if (nroots >= 0) {
5209:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5210:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5211:     if (nroots > pEnd-pStart) {
5212:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
5213:     }
5214:   }
5215:   /* Calculate new sizes, get proccess offset, and calculate point offsets */
5216:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5217:     cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
5218:     (*gsection)->atlasOff[p] = off;
5219:     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
5220:   }
5221:   MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
5222:   globalOff -= off;
5223:   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5224:     (*gsection)->atlasOff[p] += globalOff;
5225:     if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
5226:   }
5227:   /* Put in negative offsets for ghost points */
5228:   if (nroots >= 0) {
5229:     PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5230:     PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5231:     if (nroots > pEnd-pStart) {
5232:       for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
5233:     }
5234:   }
5235:   if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
5236:   PetscFree(neg);
5237:   return(0);
5238: }

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

5245:   Input Parameters:
5246:   + dm - The DMPlex object

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

5250:   Level: developer

5252: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5253: @*/
5254: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5255: {
5256:   PetscSection    coneSection, supportSection;
5257:   const PetscInt *cone, *support;
5258:   PetscInt        coneSize, c, supportSize, s;
5259:   PetscInt        pStart, pEnd, p, csize, ssize;
5260:   PetscErrorCode  ierr;

5264:   DMPlexGetConeSection(dm, &coneSection);
5265:   DMPlexGetSupportSection(dm, &supportSection);
5266:   /* Check that point p is found in the support of its cone points, and vice versa */
5267:   DMPlexGetChart(dm, &pStart, &pEnd);
5268:   for (p = pStart; p < pEnd; ++p) {
5269:     DMPlexGetConeSize(dm, p, &coneSize);
5270:     DMPlexGetCone(dm, p, &cone);
5271:     for (c = 0; c < coneSize; ++c) {
5272:       PetscBool dup = PETSC_FALSE;
5273:       PetscInt  d;
5274:       for (d = c-1; d >= 0; --d) {
5275:         if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5276:       }
5277:       DMPlexGetSupportSize(dm, cone[c], &supportSize);
5278:       DMPlexGetSupport(dm, cone[c], &support);
5279:       for (s = 0; s < supportSize; ++s) {
5280:         if (support[s] == p) break;
5281:       }
5282:       if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5283:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5284:         for (s = 0; s < coneSize; ++s) {
5285:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
5286:         }
5287:         PetscPrintf(PETSC_COMM_SELF, "\n");
5288:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
5289:         for (s = 0; s < supportSize; ++s) {
5290:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
5291:         }
5292:         PetscPrintf(PETSC_COMM_SELF, "\n");
5293:         if (dup) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
5294:         else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
5295:       }
5296:     }
5297:     DMPlexGetSupportSize(dm, p, &supportSize);
5298:     DMPlexGetSupport(dm, p, &support);
5299:     for (s = 0; s < supportSize; ++s) {
5300:       DMPlexGetConeSize(dm, support[s], &coneSize);
5301:       DMPlexGetCone(dm, support[s], &cone);
5302:       for (c = 0; c < coneSize; ++c) {
5303:         if (cone[c] == p) break;
5304:       }
5305:       if (c >= coneSize) {
5306:         PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
5307:         for (c = 0; c < supportSize; ++c) {
5308:           PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
5309:         }
5310:         PetscPrintf(PETSC_COMM_SELF, "\n");
5311:         PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
5312:         for (c = 0; c < coneSize; ++c) {
5313:           PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
5314:         }
5315:         PetscPrintf(PETSC_COMM_SELF, "\n");
5316:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
5317:       }
5318:     }
5319:   }
5320:   PetscSectionGetStorageSize(coneSection, &csize);
5321:   PetscSectionGetStorageSize(supportSection, &ssize);
5322:   if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
5323:   return(0);
5324: }

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

5331:   Input Parameters:
5332: + dm - The DMPlex object
5333: . isSimplex - Are the cells simplices or tensor products
5334: - cellHeight - Normally 0

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

5338:   Level: developer

5340: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
5341: @*/
5342: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5343: {
5344:   PetscInt       dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;

5349:   DMGetDimension(dm, &dim);
5350:   switch (dim) {
5351:   case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
5352:   case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
5353:   case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
5354:   default:
5355:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
5356:   }
5357:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5358:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5359:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5360:   cMax = cMax >= 0 ? cMax : cEnd;
5361:   for (c = cStart; c < cMax; ++c) {
5362:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5364:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5365:     for (cl = 0; cl < closureSize*2; cl += 2) {
5366:       const PetscInt p = closure[cl];
5367:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5368:     }
5369:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5370:     if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has  %d vertices != %d", c, coneSize, numCorners);
5371:   }
5372:   for (c = cMax; c < cEnd; ++c) {
5373:     PetscInt *closure = NULL, closureSize, cl, coneSize = 0;

5375:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5376:     for (cl = 0; cl < closureSize*2; cl += 2) {
5377:       const PetscInt p = closure[cl];
5378:       if ((p >= vStart) && (p < vEnd)) ++coneSize;
5379:     }
5380:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5381:     if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has  %d vertices > %d", c, coneSize, numHybridCorners);
5382:   }
5383:   return(0);
5384: }

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

5391:   Input Parameters:
5392: + dm - The DMPlex object
5393: . isSimplex - Are the cells simplices or tensor products
5394: - cellHeight - Normally 0

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

5398:   Level: developer

5400: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
5401: @*/
5402: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5403: {
5404:   PetscInt       pMax[4];
5405:   PetscInt       dim, vStart, vEnd, cStart, cEnd, c, h;

5410:   DMGetDimension(dm, &dim);
5411:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5412:   DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
5413:   for (h = cellHeight; h < dim; ++h) {
5414:     DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
5415:     for (c = cStart; c < cEnd; ++c) {
5416:       const PetscInt *cone, *ornt, *faces;
5417:       PetscInt        numFaces, faceSize, coneSize,f;
5418:       PetscInt       *closure = NULL, closureSize, cl, numCorners = 0;

5420:       if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
5421:       DMPlexGetConeSize(dm, c, &coneSize);
5422:       DMPlexGetCone(dm, c, &cone);
5423:       DMPlexGetConeOrientation(dm, c, &ornt);
5424:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5425:       for (cl = 0; cl < closureSize*2; cl += 2) {
5426:         const PetscInt p = closure[cl];
5427:         if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
5428:       }
5429:       DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
5430:       if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
5431:       for (f = 0; f < numFaces; ++f) {
5432:         PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;

5434:         DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
5435:         for (cl = 0; cl < fclosureSize*2; cl += 2) {
5436:           const PetscInt p = fclosure[cl];
5437:           if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
5438:         }
5439:         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);
5440:         for (v = 0; v < fnumCorners; ++v) {
5441:           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]);
5442:         }
5443:         DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
5444:       }
5445:       DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
5446:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5447:     }
5448:   }
5449:   return(0);
5450: }

5454: /* Pointwise interpolation
5455:      Just code FEM for now
5456:      u^f = I u^c
5457:      sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
5458:      u^f_i = sum_j psi^f_i I phi^c_j u^c_j
5459:      I_{ij} = psi^f_i phi^c_j
5460: */
5461: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
5462: {
5463:   PetscSection   gsc, gsf;
5464:   PetscInt       m, n;
5465:   void          *ctx;

5469:   /*
5470:   Loop over coarse cells
5471:     Loop over coarse basis functions
5472:       Loop over fine cells in coarse cell
5473:         Loop over fine dual basis functions
5474:           Evaluate coarse basis on fine dual basis quad points
5475:           Sum
5476:           Update local element matrix
5477:     Accumulate to interpolation matrix

5479:    Can extend PetscFEIntegrateJacobian_Basic() to do a specialized cell loop
5480:   */
5481:   DMGetDefaultGlobalSection(dmFine, &gsf);
5482:   PetscSectionGetConstrainedStorageSize(gsf, &m);
5483:   DMGetDefaultGlobalSection(dmCoarse, &gsc);
5484:   PetscSectionGetConstrainedStorageSize(gsc, &n);
5485:   /* We need to preallocate properly */
5486:   MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
5487:   MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
5488:   MatSetType(*interpolation, dmCoarse->mattype);
5489:   DMGetApplicationContext(dmFine, &ctx);
5490:   DMPlexComputeInterpolatorFEM(dmCoarse, dmFine, *interpolation, ctx);
5491:   /* Use naive scaling */
5492:   DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
5493:   return(0);
5494: }

5498: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat)
5499: {
5501:   VecScatter     ctx;

5504:   DMPlexComputeInjectorFEM(dmCoarse, dmFine, &ctx, NULL);
5505:   MatCreateScatter(PetscObjectComm((PetscObject)ctx), ctx, mat);
5506:   VecScatterDestroy(&ctx);
5507:   return(0);
5508: }

5512: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
5513: {
5514:   PetscSection   section;
5515:   IS            *bcPoints, *bcComps;
5516:   PetscBool     *isFE;
5517:   PetscInt      *bcFields, *numComp, *numDof;
5518:   PetscInt       depth, dim, numBd, numBC = 0, numFields, bd, bc = 0, f;
5519:   PetscInt       cStart, cEnd, cEndInterior;

5523:   DMGetNumFields(dm, &numFields);
5524:   /* FE and FV boundary conditions are handled slightly differently */
5525:   PetscMalloc1(numFields, &isFE);
5526:   for (f = 0; f < numFields; ++f) {
5527:     PetscObject  obj;
5528:     PetscClassId id;

5530:     DMGetField(dm, f, &obj);
5531:     PetscObjectGetClassId(obj, &id);
5532:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
5533:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
5534:     else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
5535:   }
5536:   /* Allocate boundary point storage for FEM boundaries */
5537:   DMPlexGetDepth(dm, &depth);
5538:   DMGetDimension(dm, &dim);
5539:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
5540:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
5541:   DMPlexGetNumBoundary(dm, &numBd);
5542:   for (bd = 0; bd < numBd; ++bd) {
5543:     PetscInt  field;
5544:     PetscBool isEssential;

5546:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL, NULL);
5547:     if (isFE[field] && isEssential) ++numBC;
5548:   }
5549:   /* Add ghost cell boundaries for FVM */
5550:   for (f = 0; f < numFields; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
5551:   PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);
5552:   /* Constrain ghost cells for FV */
5553:   for (f = 0; f < numFields; ++f) {
5554:     PetscInt *newidx, c;

5556:     if (isFE[f] || cEndInterior < 0) continue;
5557:     PetscMalloc1(cEnd-cEndInterior,&newidx);
5558:     for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
5559:     bcFields[bc] = f;
5560:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5561:   }
5562:   /* Handle FEM Dirichlet boundaries */
5563:   for (bd = 0; bd < numBd; ++bd) {
5564:     const char     *bdLabel;
5565:     DMLabel         label;
5566:     const PetscInt *comps;
5567:     const PetscInt *values;
5568:     PetscInt        bd2, field, numComps, numValues;
5569:     PetscBool       isEssential, duplicate = PETSC_FALSE;

5571:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
5572:     if (!isFE[field]) continue;
5573:     DMPlexGetLabel(dm, bdLabel, &label);
5574:     /* Only want to modify label once */
5575:     for (bd2 = 0; bd2 < bd; ++bd2) {
5576:       const char *bdname;
5577:       DMPlexGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
5578:       PetscStrcmp(bdname, bdLabel, &duplicate);
5579:       if (duplicate) break;
5580:     }
5581:     if (!duplicate && (isFE[field])) {
5582:       DMPlexLabelComplete(dm, label);
5583:       DMPlexLabelAddCells(dm, label);
5584:     }
5585:     /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
5586:     if (isEssential) {
5587:       PetscInt       *newidx;
5588:       PetscInt        n, newn = 0, p, v;

5590:       bcFields[bc] = field;
5591:       if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
5592:       for (v = 0; v < numValues; ++v) {
5593:         IS              tmp;
5594:         const PetscInt *idx;

5596:         DMPlexGetStratumIS(dm, bdLabel, values[v], &tmp);
5597:         if (!tmp) continue;
5598:         ISGetLocalSize(tmp, &n);
5599:         ISGetIndices(tmp, &idx);
5600:         if (isFE[field]) {
5601:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
5602:         } else {
5603:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
5604:         }
5605:         ISRestoreIndices(tmp, &idx);
5606:         ISDestroy(&tmp);
5607:       }
5608:       PetscMalloc1(newn,&newidx);
5609:       newn = 0;
5610:       for (v = 0; v < numValues; ++v) {
5611:         IS              tmp;
5612:         const PetscInt *idx;

5614:         DMPlexGetStratumIS(dm, bdLabel, values[v], &tmp);
5615:         if (!tmp) continue;
5616:         ISGetLocalSize(tmp, &n);
5617:         ISGetIndices(tmp, &idx);
5618:         if (isFE[field]) {
5619:           for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
5620:         } else {
5621:           for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
5622:         }
5623:         ISRestoreIndices(tmp, &idx);
5624:         ISDestroy(&tmp);
5625:       }
5626:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5627:     }
5628:   }
5629:   /* Handle discretization */
5630:   PetscCalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
5631:   for (f = 0; f < numFields; ++f) {
5632:     PetscObject obj;

5634:     DMGetField(dm, f, &obj);
5635:     if (isFE[f]) {
5636:       PetscFE         fe = (PetscFE) obj;
5637:       const PetscInt *numFieldDof;
5638:       PetscInt        d;

5640:       PetscFEGetNumComponents(fe, &numComp[f]);
5641:       PetscFEGetNumDof(fe, &numFieldDof);
5642:       for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
5643:     } else {
5644:       PetscFV fv = (PetscFV) obj;

5646:       PetscFVGetNumComponents(fv, &numComp[f]);
5647:       numDof[f*(dim+1)+dim] = numComp[f];
5648:     }
5649:   }
5650:   for (f = 0; f < numFields; ++f) {
5651:     PetscInt d;
5652:     for (d = 1; d < dim; ++d) {
5653:       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.");
5654:     }
5655:   }
5656:   DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, &section);
5657:   for (f = 0; f < numFields; ++f) {
5658:     PetscFE     fe;
5659:     const char *name;

5661:     DMGetField(dm, f, (PetscObject *) &fe);
5662:     PetscObjectGetName((PetscObject) fe, &name);
5663:     PetscSectionSetFieldName(section, f, name);
5664:   }
5665:   DMSetDefaultSection(dm, section);
5666:   PetscSectionDestroy(&section);
5667:   for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
5668:   PetscFree3(bcFields,bcPoints,bcComps);
5669:   PetscFree2(numComp,numDof);
5670:   PetscFree(isFE);
5671:   return(0);
5672: }

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

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

5682:   Output Parameter:
5683: . cdm - The coarse DM

5685:   Level: intermediate

5687: .seealso: DMPlexSetCoarseDM()
5688: @*/
5689: PetscErrorCode DMPlexGetCoarseDM(DM dm, DM *cdm)
5690: {
5694:   *cdm = ((DM_Plex *) dm->data)->coarseMesh;
5695:   return(0);
5696: }

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

5703:   Input Parameters:
5704: + dm - The DMPlex object
5705: - cdm - The coarse DM

5707:   Level: intermediate

5709: .seealso: DMPlexGetCoarseDM()
5710: @*/
5711: PetscErrorCode DMPlexSetCoarseDM(DM dm, DM cdm)
5712: {
5713:   DM_Plex       *mesh;

5719:   mesh = (DM_Plex *) dm->data;
5720:   DMDestroy(&mesh->coarseMesh);
5721:   mesh->coarseMesh = cdm;
5722:   PetscObjectReference((PetscObject) mesh->coarseMesh);
5723:   return(0);
5724: }

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

5733:   not collective

5735:   Input Parameters:
5736: . dm - The DMPlex object

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


5743:   Level: intermediate

5745: .seealso: DMPlexSetAnchors(), DMGetConstraints(), DMSetConstraints()
5746: @*/
5747: PetscErrorCode DMPlexGetAnchors(DM dm, PetscSection *anchorSection, IS *anchorIS)
5748: {
5749:   DM_Plex *plex = (DM_Plex *)dm->data;

5754:   if (!plex->anchorSection && !plex->anchorIS && plex->createanchors) {(*plex->createanchors)(dm);}
5755:   if (anchorSection) *anchorSection = plex->anchorSection;
5756:   if (anchorIS) *anchorIS = plex->anchorIS;
5757:   return(0);
5758: }

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

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

5770:   collective on dm

5772:   Input Parameters:
5773: + dm - The DMPlex object
5774: . anchorSection - The section that describes the mapping from constrained points to the anchor points listed in anchorIS.  Must have a local communicator (PETSC_COMM_SELF or derivative).
5775: - anchorIS - The list of all anchor points.  Must have a local communicator (PETSC_COMM_SELF or derivative).

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

5779:   Level: intermediate

5781: .seealso: DMPlexGetAnchors(), DMGetConstraints(), DMSetConstraints()
5782: @*/
5783: PetscErrorCode DMPlexSetAnchors(DM dm, PetscSection anchorSection, IS anchorIS)
5784: {
5785:   DM_Plex        *plex = (DM_Plex *)dm->data;
5786:   PetscMPIInt    result;

5791:   if (anchorSection) {
5793:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorSection),&result);
5794:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor section must have local communicator");
5795:   }
5796:   if (anchorIS) {
5798:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorIS),&result);
5799:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor IS must have local communicator");
5800:   }

5802:   PetscObjectReference((PetscObject)anchorSection);
5803:   PetscSectionDestroy(&plex->anchorSection);
5804:   plex->anchorSection = anchorSection;

5806:   PetscObjectReference((PetscObject)anchorIS);
5807:   ISDestroy(&plex->anchorIS);
5808:   plex->anchorIS = anchorIS;

5810: #if defined(PETSC_USE_DEBUG)
5811:   if (anchorIS && anchorSection) {
5812:     PetscInt size, a, pStart, pEnd;
5813:     const PetscInt *anchors;

5815:     PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5816:     ISGetLocalSize(anchorIS,&size);
5817:     ISGetIndices(anchorIS,&anchors);
5818:     for (a = 0; a < size; a++) {
5819:       PetscInt p;

5821:       p = anchors[a];
5822:       if (p >= pStart && p < pEnd) {
5823:         PetscInt dof;

5825:         PetscSectionGetDof(anchorSection,p,&dof);
5826:         if (dof) {
5827:           PetscErrorCode ierr2;

5829:           ierr2 = ISRestoreIndices(anchorIS,&anchors);CHKERRQ(ierr2);
5830:           SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point %d cannot be constrained and an anchor",p);
5831:         }
5832:       }
5833:     }
5834:     ISRestoreIndices(anchorIS,&anchors);
5835:   }
5836: #endif
5837:   /* reset the generic constraints */
5838:   DMSetDefaultConstraints(dm,NULL,NULL);
5839:   return(0);
5840: }

5844: static PetscErrorCode DMPlexCreateConstraintSection_Anchors(DM dm, PetscSection section, PetscSection *cSec)
5845: {
5846:   PetscSection anchorSection;
5847:   PetscInt pStart, pEnd, p, dof, numFields, f;

5852:   DMPlexGetAnchors(dm,&anchorSection,NULL);
5853:   PetscSectionCreate(PETSC_COMM_SELF,cSec);
5854:   PetscSectionGetNumFields(section,&numFields);
5855:   PetscSectionSetNumFields(*cSec,numFields);
5856:   PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5857:   PetscSectionSetChart(*cSec,pStart,pEnd);
5858:   for (p = pStart; p < pEnd; p++) {
5859:     PetscSectionGetDof(anchorSection,p,&dof);
5860:     if (dof) {
5861:       PetscSectionGetDof(section,p,&dof);
5862:       PetscSectionSetDof(*cSec,p,dof);
5863:       for (f = 0; f < numFields; f++) {
5864:         PetscSectionGetFieldDof(section,p,f,&dof);
5865:         PetscSectionSetFieldDof(*cSec,p,f,dof);
5866:       }
5867:     }
5868:   }
5869:   PetscSectionSetUp(*cSec);
5870:   return(0);
5871: }

5875: static PetscErrorCode DMPlexCreateConstraintMatrix_Anchors(DM dm, PetscSection section, PetscSection cSec, Mat *cMat)
5876: {
5877:   PetscSection aSec;
5878:   PetscInt pStart, pEnd, p, dof, aDof, aOff, off, nnz, annz, m, n, q, a, offset, *i, *j;
5879:   const PetscInt *anchors;
5880:   PetscInt numFields, f;
5881:   IS aIS;

5886:   PetscSectionGetStorageSize(cSec, &m);
5887:   PetscSectionGetStorageSize(section, &n);
5888:   MatCreate(PETSC_COMM_SELF,cMat);
5889:   MatSetSizes(*cMat,m,n,m,n);
5890:   MatSetType(*cMat,MATSEQAIJ);
5891:   DMPlexGetAnchors(dm,&aSec,&aIS);
5892:   ISGetIndices(aIS,&anchors);
5893:   PetscSectionGetChart(aSec,&pStart,&pEnd);
5894:   PetscMalloc1(m+1,&i);
5895:   i[0] = 0;
5896:   PetscSectionGetNumFields(section,&numFields);
5897:   for (p = pStart; p < pEnd; p++) {
5898:     PetscSectionGetDof(aSec,p,&dof);
5899:     if (!dof) continue;
5900:     PetscSectionGetOffset(aSec,p,&off);
5901:     if (numFields) {
5902:       for (f = 0; f < numFields; f++) {
5903:         annz = 0;
5904:         for (q = 0; q < dof; q++) {
5905:           a = anchors[off + q];
5906:           PetscSectionGetFieldDof(section,a,f,&aDof);
5907:           annz += aDof;
5908:         }
5909:         PetscSectionGetFieldDof(cSec,p,f,&dof);
5910:         PetscSectionGetFieldOffset(cSec,p,f,&off);
5911:         for (q = 0; q < dof; q++) {
5912:           i[off + q + 1] = i[off + q] + annz;
5913:         }
5914:       }
5915:     }
5916:     else {
5917:       annz = 0;
5918:       for (q = 0; q < dof; q++) {
5919:         a = anchors[off + q];
5920:         PetscSectionGetDof(section,a,&aDof);
5921:         annz += aDof;
5922:       }
5923:       PetscSectionGetDof(cSec,p,&dof);
5924:       PetscSectionGetOffset(cSec,p,&off);
5925:       for (q = 0; q < dof; q++) {
5926:         i[off + q + 1] = i[off + q] + annz;
5927:       }
5928:     }
5929:   }
5930:   nnz = i[m];
5931:   PetscMalloc1(nnz,&j);
5932:   offset = 0;
5933:   for (p = pStart; p < pEnd; p++) {
5934:     if (numFields) {
5935:       for (f = 0; f < numFields; f++) {
5936:         PetscSectionGetFieldDof(cSec,p,f,&dof);
5937:         for (q = 0; q < dof; q++) {
5938:           PetscInt rDof, rOff, r;
5939:           PetscSectionGetDof(aSec,p,&rDof);
5940:           PetscSectionGetOffset(aSec,p,&rOff);
5941:           for (r = 0; r < rDof; r++) {
5942:             PetscInt s;

5944:             a = anchors[rOff + r];

5946:             PetscSectionGetFieldDof(section,a,f,&aDof);
5947:             PetscSectionGetFieldOffset(section,a,f,&aOff);
5948:             for (s = 0; s < aDof; s++) {
5949:               j[offset++] = aOff + s;
5950:             }
5951:           }
5952:         }
5953:       }
5954:     }
5955:     else {
5956:       PetscSectionGetDof(cSec,p,&dof);
5957:       for (q = 0; q < dof; q++) {
5958:         PetscInt rDof, rOff, r;
5959:         PetscSectionGetDof(aSec,p,&rDof);
5960:         PetscSectionGetOffset(aSec,p,&rOff);
5961:         for (r = 0; r < rDof; r++) {
5962:           PetscInt s;

5964:           a = anchors[rOff + r];

5966:           PetscSectionGetDof(section,a,&aDof);
5967:           PetscSectionGetOffset(section,a,&aOff);
5968:           for (s = 0; s < aDof; s++) {
5969:             j[offset++] = aOff + s;
5970:           }
5971:         }
5972:       }
5973:     }
5974:   }
5975:   MatSeqAIJSetPreallocationCSR(*cMat,i,j,NULL);
5976:   PetscFree(i);
5977:   PetscFree(j);
5978:   ISRestoreIndices(aIS,&anchors);
5979:   return(0);
5980: }

5984: PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm)
5985: {
5986:   DM_Plex        *plex = (DM_Plex *)dm->data;
5987:   PetscSection   anchorSection, section, cSec;
5988:   Mat            cMat;

5993:   DMPlexGetAnchors(dm,&anchorSection,NULL);
5994:   if (anchorSection) {
5995:     PetscDS  ds;
5996:     PetscInt nf;

5998:     DMGetDefaultSection(dm,&section);
5999:     DMPlexCreateConstraintSection_Anchors(dm,section,&cSec);
6000:     DMPlexCreateConstraintMatrix_Anchors(dm,section,cSec,&cMat);
6001:     DMGetDS(dm,&ds);
6002:     PetscDSGetNumFields(ds,&nf);
6003:     if (nf && plex->computeanchormatrix) {(*plex->computeanchormatrix)(dm,section,cSec,cMat);}
6004:     DMSetDefaultConstraints(dm,cSec,cMat);
6005:     PetscSectionDestroy(&cSec);
6006:     MatDestroy(&cMat);
6007:   }
6008:   return(0);
6009: }