Actual source code: plex.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petsc/private/isimpl.h>
3: #include <petsc/private/vecimpl.h>
4: #include <petscsf.h>
5: #include <petscds.h>
7: /* Logging support */
8: PetscLogEvent DMPLEX_Interpolate, PETSCPARTITIONER_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeOverlap, DMPLEX_DistributeField, DMPLEX_DistributeData, DMPLEX_Migrate, DMPLEX_GlobalToNaturalBegin, DMPLEX_GlobalToNaturalEnd, DMPLEX_NaturalToGlobalBegin, DMPLEX_NaturalToGlobalEnd, DMPLEX_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM, DMPLEX_InterpolatorFEM, DMPLEX_InjectorFEM, DMPLEX_IntegralFEM, DMPLEX_CreateGmsh;
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, PETSC_TRUE, v, 0.0, NULL, NULL, NULL);}
75: }
76: if (isvtk) {
77: PetscSection section;
78: PetscViewerVTKFieldType ft;
79: PetscInt pStart, pEnd;
81: DMGetDefaultSection(dm, §ion);
82: DMPlexGetFieldType_Internal(dm, section, PETSC_DETERMINE, &pStart, &pEnd, &ft);
83: PetscObjectReference((PetscObject) v); /* viewer drops reference */
84: PetscViewerVTKAddField(viewer, (PetscObject) dm, DMPlexVTKWriteAll, ft, (PetscObject) v);
85: } else if (ishdf5) {
86: #if defined(PETSC_HAVE_HDF5)
87: VecView_Plex_Local_HDF5(v, viewer);
88: #else
89: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
90: #endif
91: } else {
92: if (isseq) {VecView_Seq(v, viewer);}
93: else {VecView_MPI(v, viewer);}
94: }
95: return(0);
96: }
100: PetscErrorCode VecView_Plex(Vec v, PetscViewer viewer)
101: {
102: DM dm;
103: PetscBool isvtk, ishdf5, isseq;
107: VecGetDM(v, &dm);
108: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
109: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
110: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
111: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
112: if (isvtk) {
113: Vec locv;
114: const char *name;
116: DMGetLocalVector(dm, &locv);
117: PetscObjectGetName((PetscObject) v, &name);
118: PetscObjectSetName((PetscObject) locv, name);
119: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
120: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
121: VecView_Plex_Local(locv, viewer);
122: DMRestoreLocalVector(dm, &locv);
123: } else if (ishdf5) {
124: #if defined(PETSC_HAVE_HDF5)
125: VecView_Plex_HDF5(v, viewer);
126: #else
127: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
128: #endif
129: } else {
130: if (isseq) {VecView_Seq(v, viewer);}
131: else {VecView_MPI(v, viewer);}
132: }
133: return(0);
134: }
138: PetscErrorCode VecView_Plex_Native(Vec originalv, PetscViewer viewer)
139: {
140: DM dm;
141: MPI_Comm comm;
142: PetscViewerFormat format;
143: Vec v;
144: PetscBool isvtk, ishdf5;
145: PetscErrorCode ierr;
148: VecGetDM(originalv, &dm);
149: PetscObjectGetComm((PetscObject) originalv, &comm);
150: if (!dm) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
151: PetscViewerGetFormat(viewer, &format);
152: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
153: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
154: if (format == PETSC_VIEWER_NATIVE) {
155: const char *vecname;
156: PetscInt n, nroots;
158: if (dm->sfNatural) {
159: VecGetLocalSize(originalv, &n);
160: PetscSFGetGraph(dm->sfNatural, &nroots, NULL, NULL, NULL);
161: if (n == nroots) {
162: DMGetGlobalVector(dm, &v);
163: DMPlexGlobalToNaturalBegin(dm, originalv, v);
164: DMPlexGlobalToNaturalEnd(dm, originalv, v);
165: PetscObjectGetName((PetscObject) originalv, &vecname);
166: PetscObjectSetName((PetscObject) v, vecname);
167: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "DM global to natural SF only handles global vectors");
168: } else SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created");
169: } else {
170: /* we are viewing a natural DMPlex vec. */
171: v = originalv;
172: }
173: if (ishdf5) {
174: #if defined(PETSC_HAVE_HDF5)
175: VecView_Plex_HDF5_Native(v, viewer);
176: #else
177: SETERRQ(comm, PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
178: #endif
179: } else if (isvtk) {
180: SETERRQ(comm, PETSC_ERR_SUP, "VTK format does not support viewing in natural order. Please switch to HDF5.");
181: } else {
182: PetscBool isseq;
184: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
185: if (isseq) {VecView_Seq(v, viewer);}
186: else {VecView_MPI(v, viewer);}
187: }
188: if (format == PETSC_VIEWER_NATIVE) {DMRestoreGlobalVector(dm, &v);}
189: return(0);
190: }
194: PetscErrorCode VecLoad_Plex_Local(Vec v, PetscViewer viewer)
195: {
196: DM dm;
197: PetscBool ishdf5;
201: VecGetDM(v, &dm);
202: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
203: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
204: if (ishdf5) {
205: DM dmBC;
206: Vec gv;
207: const char *name;
209: DMGetOutputDM(dm, &dmBC);
210: DMGetGlobalVector(dmBC, &gv);
211: PetscObjectGetName((PetscObject) v, &name);
212: PetscObjectSetName((PetscObject) gv, name);
213: VecLoad_Default(gv, viewer);
214: DMGlobalToLocalBegin(dmBC, gv, INSERT_VALUES, v);
215: DMGlobalToLocalEnd(dmBC, gv, INSERT_VALUES, v);
216: DMRestoreGlobalVector(dmBC, &gv);
217: } else {
218: VecLoad_Default(v, viewer);
219: }
220: return(0);
221: }
225: PetscErrorCode VecLoad_Plex(Vec v, PetscViewer viewer)
226: {
227: DM dm;
228: PetscBool ishdf5;
232: VecGetDM(v, &dm);
233: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
234: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
235: if (ishdf5) {
236: #if defined(PETSC_HAVE_HDF5)
237: VecLoad_Plex_HDF5(v, viewer);
238: #else
239: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
240: #endif
241: } else {
242: VecLoad_Default(v, viewer);
243: }
244: return(0);
245: }
249: PetscErrorCode VecLoad_Plex_Native(Vec originalv, PetscViewer viewer)
250: {
251: DM dm;
252: PetscViewerFormat format;
253: PetscBool ishdf5;
254: PetscErrorCode ierr;
257: VecGetDM(originalv, &dm);
258: if (!dm) SETERRQ(PetscObjectComm((PetscObject) originalv), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
259: PetscViewerGetFormat(viewer, &format);
260: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
261: if (format == PETSC_VIEWER_NATIVE) {
262: if (dm->sfNatural) {
263: if (ishdf5) {
264: #if defined(PETSC_HAVE_HDF5)
265: Vec v;
266: const char *vecname;
268: DMGetGlobalVector(dm, &v);
269: PetscObjectGetName((PetscObject) originalv, &vecname);
270: PetscObjectSetName((PetscObject) v, vecname);
271: VecLoad_Plex_HDF5_Native(v, viewer);
272: DMPlexNaturalToGlobalBegin(dm, v, originalv);
273: DMPlexNaturalToGlobalEnd(dm, v, originalv);
274: DMRestoreGlobalVector(dm, &v);
275: #else
276: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
277: #endif
278: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Reading in natural order is not supported for anything but HDF5.");
279: }
280: }
281: return(0);
282: }
286: PetscErrorCode DMPlexView_Ascii_Geometry(DM dm, PetscViewer viewer)
287: {
288: PetscSection coordSection;
289: Vec coordinates;
290: DMLabel depthLabel;
291: const char *name[4];
292: const PetscScalar *a;
293: PetscInt dim, pStart, pEnd, cStart, cEnd, c;
294: PetscErrorCode ierr;
297: DMGetDimension(dm, &dim);
298: DMGetCoordinatesLocal(dm, &coordinates);
299: DMGetCoordinateSection(dm, &coordSection);
300: DMPlexGetDepthLabel(dm, &depthLabel);
301: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
302: PetscSectionGetChart(coordSection, &pStart, &pEnd);
303: VecGetArrayRead(coordinates, &a);
304: name[0] = "vertex";
305: name[1] = "edge";
306: name[dim-1] = "face";
307: name[dim] = "cell";
308: for (c = cStart; c < cEnd; ++c) {
309: PetscInt *closure = NULL;
310: PetscInt closureSize, cl;
312: PetscViewerASCIIPrintf(viewer, "Geometry for cell %D:\n", c);
313: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
314: PetscViewerASCIIPushTab(viewer);
315: for (cl = 0; cl < closureSize*2; cl += 2) {
316: PetscInt point = closure[cl], depth, dof, off, d, p;
318: if ((point < pStart) || (point >= pEnd)) continue;
319: PetscSectionGetDof(coordSection, point, &dof);
320: if (!dof) continue;
321: DMLabelGetValue(depthLabel, point, &depth);
322: PetscSectionGetOffset(coordSection, point, &off);
323: PetscViewerASCIIPrintf(viewer, "%s %D coords:", name[depth], point);
324: for (p = 0; p < dof/dim; ++p) {
325: PetscViewerASCIIPrintf(viewer, " (");
326: for (d = 0; d < dim; ++d) {
327: if (d > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
328: PetscViewerASCIIPrintf(viewer, "%g", PetscRealPart(a[off+p*dim+d]));
329: }
330: PetscViewerASCIIPrintf(viewer, ")");
331: }
332: PetscViewerASCIIPrintf(viewer, "\n");
333: }
334: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
335: PetscViewerASCIIPopTab(viewer);
336: }
337: VecRestoreArrayRead(coordinates, &a);
338: return(0);
339: }
343: PetscErrorCode DMPlexView_Ascii(DM dm, PetscViewer viewer)
344: {
345: DM_Plex *mesh = (DM_Plex*) dm->data;
346: DM cdm;
347: DMLabel markers;
348: PetscSection coordSection;
349: Vec coordinates;
350: PetscViewerFormat format;
351: PetscErrorCode ierr;
354: DMGetCoordinateDM(dm, &cdm);
355: DMGetDefaultSection(cdm, &coordSection);
356: DMGetCoordinatesLocal(dm, &coordinates);
357: PetscViewerGetFormat(viewer, &format);
358: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
359: const char *name;
360: PetscInt maxConeSize, maxSupportSize;
361: PetscInt pStart, pEnd, p;
362: PetscMPIInt rank, size;
364: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
365: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
366: PetscObjectGetName((PetscObject) dm, &name);
367: DMPlexGetChart(dm, &pStart, &pEnd);
368: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
369: PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
370: PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
371: PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
372: PetscViewerASCIIPushSynchronized(viewer);
373: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Max sizes cone: %D support: %D\n", rank,maxConeSize, maxSupportSize);
374: for (p = pStart; p < pEnd; ++p) {
375: PetscInt dof, off, s;
377: PetscSectionGetDof(mesh->supportSection, p, &dof);
378: PetscSectionGetOffset(mesh->supportSection, p, &off);
379: for (s = off; s < off+dof; ++s) {
380: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D ----> %D\n", rank, p, mesh->supports[s]);
381: }
382: }
383: PetscViewerFlush(viewer);
384: PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
385: for (p = pStart; p < pEnd; ++p) {
386: PetscInt dof, off, c;
388: PetscSectionGetDof(mesh->coneSection, p, &dof);
389: PetscSectionGetOffset(mesh->coneSection, p, &off);
390: for (c = off; c < off+dof; ++c) {
391: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D <---- %D (%D)\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]);
392: }
393: }
394: PetscViewerFlush(viewer);
395: PetscViewerASCIIPopSynchronized(viewer);
396: PetscSectionGetChart(coordSection, &pStart, NULL);
397: if (pStart >= 0) {PetscSectionVecView(coordSection, coordinates, viewer);}
398: DMGetLabel(dm, "marker", &markers);
399: DMLabelView(markers,viewer);
400: if (size > 1) {
401: PetscSF sf;
403: DMGetPointSF(dm, &sf);
404: PetscSFView(sf, viewer);
405: }
406: PetscViewerFlush(viewer);
407: } else if (format == PETSC_VIEWER_ASCII_LATEX) {
408: const char *name, *color;
409: const char *defcolors[3] = {"gray", "orange", "green"};
410: const char *deflcolors[4] = {"blue", "cyan", "red", "magenta"};
411: PetscReal scale = 2.0;
412: PetscBool useNumbers = PETSC_TRUE, useLabels, useColors;
413: double tcoords[3];
414: PetscScalar *coords;
415: PetscInt numLabels, l, numColors, numLColors, dim, depth, cStart, cEnd, c, vStart, vEnd, v, eStart = 0, eEnd = 0, e, p;
416: PetscMPIInt rank, size;
417: char **names, **colors, **lcolors;
419: DMGetDimension(dm, &dim);
420: DMPlexGetDepth(dm, &depth);
421: DMGetNumLabels(dm, &numLabels);
422: numLabels = PetscMax(numLabels, 10);
423: numColors = 10;
424: numLColors = 10;
425: PetscCalloc3(numLabels, &names, numColors, &colors, numLColors, &lcolors);
426: PetscOptionsGetReal(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_scale", &scale, NULL);
427: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_numbers", &useNumbers, NULL);
428: PetscOptionsGetStringArray(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_labels", names, &numLabels, &useLabels);
429: if (!useLabels) numLabels = 0;
430: PetscOptionsGetStringArray(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_colors", colors, &numColors, &useColors);
431: if (!useColors) {
432: numColors = 3;
433: for (c = 0; c < numColors; ++c) {PetscStrallocpy(defcolors[c], &colors[c]);}
434: }
435: PetscOptionsGetStringArray(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_lcolors", lcolors, &numLColors, &useColors);
436: if (!useColors) {
437: numLColors = 4;
438: for (c = 0; c < numLColors; ++c) {PetscStrallocpy(deflcolors[c], &lcolors[c]);}
439: }
440: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
441: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
442: PetscObjectGetName((PetscObject) dm, &name);
443: PetscViewerASCIIPrintf(viewer, "\
444: \\documentclass[tikz]{standalone}\n\n\
445: \\usepackage{pgflibraryshapes}\n\
446: \\usetikzlibrary{backgrounds}\n\
447: \\usetikzlibrary{arrows}\n\
448: \\begin{document}\n");
449: if (size > 1) {
450: PetscViewerASCIIPrintf(viewer, "%s for process ", name);
451: for (p = 0; p < size; ++p) {
452: if (p > 0 && p == size-1) {
453: PetscViewerASCIIPrintf(viewer, ", and ", colors[p%numColors], p);
454: } else if (p > 0) {
455: PetscViewerASCIIPrintf(viewer, ", ", colors[p%numColors], p);
456: }
457: PetscViewerASCIIPrintf(viewer, "{\\textcolor{%s}%D}", colors[p%numColors], p);
458: }
459: PetscViewerASCIIPrintf(viewer, ".\n\n\n");
460: }
461: PetscViewerASCIIPrintf(viewer, "\\begin{tikzpicture}[scale = %g,font=\\fontsize{8}{8}\\selectfont]\n", 1.0);
462: /* Plot vertices */
463: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
464: VecGetArray(coordinates, &coords);
465: PetscViewerASCIIPushSynchronized(viewer);
466: for (v = vStart; v < vEnd; ++v) {
467: PetscInt off, dof, d;
468: PetscBool isLabeled = PETSC_FALSE;
470: PetscSectionGetDof(coordSection, v, &dof);
471: PetscSectionGetOffset(coordSection, v, &off);
472: PetscViewerASCIISynchronizedPrintf(viewer, "\\path (");
473: if (PetscUnlikely(dof > 3)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"coordSection vertex %D has dof %D > 3",v,dof);
474: for (d = 0; d < dof; ++d) {
475: tcoords[d] = (double) (scale*PetscRealPart(coords[off+d]));
476: tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
477: }
478: /* Rotate coordinates since PGF makes z point out of the page instead of up */
479: if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
480: for (d = 0; d < dof; ++d) {
481: if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
482: PetscViewerASCIISynchronizedPrintf(viewer, "%g", tcoords[d]);
483: }
484: color = colors[rank%numColors];
485: for (l = 0; l < numLabels; ++l) {
486: PetscInt val;
487: DMGetLabelValue(dm, names[l], v, &val);
488: if (val >= 0) {color = lcolors[l%numLColors]; isLabeled = PETSC_TRUE; break;}
489: }
490: if (useNumbers) {
491: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [draw,shape=circle,color=%s] {%D};\n", v, rank, color, v);
492: } else {
493: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [fill,inner sep=%dpt,shape=circle,color=%s] {};\n", v, rank, !isLabeled ? 1 : 2, color);
494: }
495: }
496: VecRestoreArray(coordinates, &coords);
497: PetscViewerFlush(viewer);
498: /* Plot edges */
499: if (depth > 1) {DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);}
500: if (dim < 3 && useNumbers) {
501: VecGetArray(coordinates, &coords);
502: PetscViewerASCIIPrintf(viewer, "\\path\n");
503: for (e = eStart; e < eEnd; ++e) {
504: const PetscInt *cone;
505: PetscInt coneSize, offA, offB, dof, d;
507: DMPlexGetConeSize(dm, e, &coneSize);
508: if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %D cone should have two vertices, not %D", e, coneSize);
509: DMPlexGetCone(dm, e, &cone);
510: PetscSectionGetDof(coordSection, cone[0], &dof);
511: PetscSectionGetOffset(coordSection, cone[0], &offA);
512: PetscSectionGetOffset(coordSection, cone[1], &offB);
513: PetscViewerASCIISynchronizedPrintf(viewer, "(");
514: for (d = 0; d < dof; ++d) {
515: tcoords[d] = (double) (scale*PetscRealPart(coords[offA+d]+coords[offB+d]));
516: tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
517: }
518: /* Rotate coordinates since PGF makes z point out of the page instead of up */
519: if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
520: for (d = 0; d < dof; ++d) {
521: if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
522: PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)tcoords[d]);
523: }
524: color = colors[rank%numColors];
525: for (l = 0; l < numLabels; ++l) {
526: PetscInt val;
527: DMGetLabelValue(dm, names[l], v, &val);
528: if (val >= 0) {color = lcolors[l%numLColors]; break;}
529: }
530: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [draw,shape=circle,color=%s] {%D} --\n", e, rank, color, e);
531: }
532: VecRestoreArray(coordinates, &coords);
533: PetscViewerFlush(viewer);
534: PetscViewerASCIIPrintf(viewer, "(0,0);\n");
535: }
536: /* Plot cells */
537: if (dim == 3 || !useNumbers) {
538: for (e = eStart; e < eEnd; ++e) {
539: const PetscInt *cone;
541: color = colors[rank%numColors];
542: for (l = 0; l < numLabels; ++l) {
543: PetscInt val;
544: DMGetLabelValue(dm, names[l], e, &val);
545: if (val >= 0) {color = lcolors[l%numLColors]; break;}
546: }
547: DMPlexGetCone(dm, e, &cone);
548: PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] (%D_%d) -- (%D_%d);\n", color, cone[0], rank, cone[1], rank);
549: }
550: } else {
551: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
552: for (c = cStart; c < cEnd; ++c) {
553: PetscInt *closure = NULL;
554: PetscInt closureSize, firstPoint = -1;
556: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
557: PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] ", colors[rank%numColors]);
558: for (p = 0; p < closureSize*2; p += 2) {
559: const PetscInt point = closure[p];
561: if ((point < vStart) || (point >= vEnd)) continue;
562: if (firstPoint >= 0) {PetscViewerASCIISynchronizedPrintf(viewer, " -- ");}
563: PetscViewerASCIISynchronizedPrintf(viewer, "(%D_%d)", point, rank);
564: if (firstPoint < 0) firstPoint = point;
565: }
566: /* Why doesn't this work? PetscViewerASCIISynchronizedPrintf(viewer, " -- cycle;\n"); */
567: PetscViewerASCIISynchronizedPrintf(viewer, " -- (%D_%d);\n", firstPoint, rank);
568: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
569: }
570: }
571: PetscViewerFlush(viewer);
572: PetscViewerASCIIPopSynchronized(viewer);
573: PetscViewerASCIIPrintf(viewer, "\\end{tikzpicture}\n");
574: PetscViewerASCIIPrintf(viewer, "\\end{document}\n", name);
575: for (l = 0; l < numLabels; ++l) {PetscFree(names[l]);}
576: for (c = 0; c < numColors; ++c) {PetscFree(colors[c]);}
577: for (c = 0; c < numLColors; ++c) {PetscFree(lcolors[c]);}
578: PetscFree3(names, colors, lcolors);
579: } else {
580: MPI_Comm comm;
581: PetscInt *sizes, *hybsizes;
582: PetscInt locDepth, depth, dim, d, pMax[4];
583: PetscInt pStart, pEnd, p;
584: PetscInt numLabels, l;
585: const char *name;
586: PetscMPIInt size;
588: PetscObjectGetComm((PetscObject)dm,&comm);
589: MPI_Comm_size(comm, &size);
590: DMGetDimension(dm, &dim);
591: PetscObjectGetName((PetscObject) dm, &name);
592: if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dim);}
593: else {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dim);}
594: DMPlexGetDepth(dm, &locDepth);
595: MPIU_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
596: DMPlexGetHybridBounds(dm, &pMax[depth], depth > 0 ? &pMax[depth-1] : NULL, &pMax[1], &pMax[0]);
597: PetscMalloc2(size,&sizes,size,&hybsizes);
598: if (depth == 1) {
599: DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd);
600: pEnd = pEnd - pStart;
601: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
602: PetscViewerASCIIPrintf(viewer, " %d-cells:", 0);
603: for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
604: PetscViewerASCIIPrintf(viewer, "\n");
605: DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd);
606: pEnd = pEnd - pStart;
607: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
608: PetscViewerASCIIPrintf(viewer, " %D-cells:", dim);
609: for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
610: PetscViewerASCIIPrintf(viewer, "\n");
611: } else {
612: for (d = 0; d <= dim; d++) {
613: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
614: pEnd -= pStart;
615: pMax[d] -= pStart;
616: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
617: MPI_Gather(&pMax[d], 1, MPIU_INT, hybsizes, 1, MPIU_INT, 0, comm);
618: PetscViewerASCIIPrintf(viewer, " %D-cells:", d);
619: for (p = 0; p < size; ++p) {
620: if (hybsizes[p] >= 0) {PetscViewerASCIIPrintf(viewer, " %D (%D)", sizes[p], sizes[p] - hybsizes[p]);}
621: else {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
622: }
623: PetscViewerASCIIPrintf(viewer, "\n");
624: }
625: }
626: PetscFree2(sizes,hybsizes);
627: DMGetNumLabels(dm, &numLabels);
628: if (numLabels) {PetscViewerASCIIPrintf(viewer, "Labels:\n");}
629: for (l = 0; l < numLabels; ++l) {
630: DMLabel label;
631: const char *name;
632: IS valueIS;
633: const PetscInt *values;
634: PetscInt numValues, v;
636: DMGetLabelName(dm, l, &name);
637: DMGetLabel(dm, name, &label);
638: DMLabelGetNumValues(label, &numValues);
639: PetscViewerASCIIPrintf(viewer, " %s: %D strata of sizes (", name, numValues);
640: DMLabelGetValueIS(label, &valueIS);
641: ISGetIndices(valueIS, &values);
642: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
643: for (v = 0; v < numValues; ++v) {
644: PetscInt size;
646: DMLabelGetStratumSize(label, values[v], &size);
647: if (v > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
648: PetscViewerASCIIPrintf(viewer, "%D", size);
649: }
650: PetscViewerASCIIPrintf(viewer, ")\n");
651: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
652: ISRestoreIndices(valueIS, &values);
653: ISDestroy(&valueIS);
654: }
655: DMGetCoarseDM(dm, &cdm);
656: if (cdm) {
657: PetscViewerASCIIPushTab(viewer);
658: DMPlexView_Ascii(cdm, viewer);
659: PetscViewerASCIIPopTab(viewer);
660: }
661: }
662: return(0);
663: }
667: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
668: {
669: PetscBool iascii, ishdf5, isvtk;
675: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
676: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
677: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
678: if (iascii) {
679: DMPlexView_Ascii(dm, viewer);
680: } else if (ishdf5) {
681: #if defined(PETSC_HAVE_HDF5)
682: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
683: DMPlexView_HDF5(dm, viewer);
684: PetscViewerPopFormat(viewer);
685: #else
686: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
687: #endif
688: }
689: else if (isvtk) {
690: DMPlexVTKWriteAll((PetscObject) dm,viewer);
691: }
692: return(0);
693: }
697: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
698: {
699: PetscBool isbinary, ishdf5;
705: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
706: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
707: if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");}
708: else if (ishdf5) {
709: #if defined(PETSC_HAVE_HDF5)
710: DMPlexLoad_HDF5(dm, viewer);
711: #else
712: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
713: #endif
714: }
715: return(0);
716: }
721: PetscErrorCode DMDestroy_Plex(DM dm)
722: {
723: DM_Plex *mesh = (DM_Plex*) dm->data;
727: if (--mesh->refct > 0) return(0);
728: PetscSectionDestroy(&mesh->coneSection);
729: PetscFree(mesh->cones);
730: PetscFree(mesh->coneOrientations);
731: PetscSectionDestroy(&mesh->supportSection);
732: PetscFree(mesh->supports);
733: PetscFree(mesh->facesTmp);
734: PetscFree(mesh->tetgenOpts);
735: PetscFree(mesh->triangleOpts);
736: PetscPartitionerDestroy(&mesh->partitioner);
737: DMLabelDestroy(&mesh->subpointMap);
738: ISDestroy(&mesh->globalVertexNumbers);
739: ISDestroy(&mesh->globalCellNumbers);
740: PetscSectionDestroy(&mesh->anchorSection);
741: ISDestroy(&mesh->anchorIS);
742: PetscSectionDestroy(&mesh->parentSection);
743: PetscFree(mesh->parents);
744: PetscFree(mesh->childIDs);
745: PetscSectionDestroy(&mesh->childSection);
746: PetscFree(mesh->children);
747: DMDestroy(&mesh->referenceTree);
748: PetscGridHashDestroy(&mesh->lbox);
749: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
750: PetscFree(mesh);
751: return(0);
752: }
756: PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
757: {
758: PetscSection sectionGlobal;
759: PetscInt bs = -1;
760: PetscInt localSize;
761: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock;
763: MatType mtype;
764: ISLocalToGlobalMapping ltog;
767: MatInitializePackage();
768: mtype = dm->mattype;
769: DMGetDefaultGlobalSection(dm, §ionGlobal);
770: /* PetscSectionGetStorageSize(sectionGlobal, &localSize); */
771: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
772: MatCreate(PetscObjectComm((PetscObject)dm), J);
773: MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
774: MatSetType(*J, mtype);
775: MatSetFromOptions(*J);
776: PetscStrcmp(mtype, MATSHELL, &isShell);
777: PetscStrcmp(mtype, MATBAIJ, &isBlock);
778: PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
779: PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
780: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
781: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
782: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
783: if (!isShell) {
784: PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
785: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin;
787: if (bs < 0) {
788: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
789: PetscInt pStart, pEnd, p, dof, cdof;
791: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
792: for (p = pStart; p < pEnd; ++p) {
793: PetscInt bdof;
795: PetscSectionGetDof(sectionGlobal, p, &dof);
796: PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
797: bdof = PetscAbs(dof) - cdof;
798: if (bdof) {
799: if (bs < 0) {
800: bs = bdof;
801: } else if (bs != bdof) {
802: /* Layout does not admit a pointwise block size */
803: bs = 1;
804: break;
805: }
806: }
807: }
808: /* Must have same blocksize on all procs (some might have no points) */
809: bsLocal = bs;
810: MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
811: bsLocal = bs < 0 ? bsMax : bs;
812: MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
813: if (bsMin != bsMax) {
814: bs = 1;
815: } else {
816: bs = bsMax;
817: }
818: } else {
819: bs = 1;
820: }
821: }
822: PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
823: DMPlexPreallocateOperator(dm, bs, dnz, onz, dnzu, onzu, *J, fillMatrix);
824: PetscFree4(dnz, onz, dnzu, onzu);
826: /* Set localtoglobalmapping on the matrix for MatSetValuesLocal() to work */
827: DMGetLocalToGlobalMapping(dm,<og);
828: MatSetLocalToGlobalMapping(*J,ltog,ltog);
829: }
830: return(0);
831: }
835: /*@
836: DMPlexGetChart - Return the interval for all mesh points [pStart, pEnd)
838: Not collective
840: Input Parameter:
841: . mesh - The DMPlex
843: Output Parameters:
844: + pStart - The first mesh point
845: - pEnd - The upper bound for mesh points
847: Level: beginner
849: .seealso: DMPlexCreate(), DMPlexSetChart()
850: @*/
851: PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
852: {
853: DM_Plex *mesh = (DM_Plex*) dm->data;
858: PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
859: return(0);
860: }
864: /*@
865: DMPlexSetChart - Set the interval for all mesh points [pStart, pEnd)
867: Not collective
869: Input Parameters:
870: + mesh - The DMPlex
871: . pStart - The first mesh point
872: - pEnd - The upper bound for mesh points
874: Output Parameters:
876: Level: beginner
878: .seealso: DMPlexCreate(), DMPlexGetChart()
879: @*/
880: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
881: {
882: DM_Plex *mesh = (DM_Plex*) dm->data;
887: PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
888: PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
889: return(0);
890: }
894: /*@
895: DMPlexGetConeSize - Return the number of in-edges for this point in the Sieve DAG
897: Not collective
899: Input Parameters:
900: + mesh - The DMPlex
901: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
903: Output Parameter:
904: . size - The cone size for point p
906: Level: beginner
908: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
909: @*/
910: PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size)
911: {
912: DM_Plex *mesh = (DM_Plex*) dm->data;
918: PetscSectionGetDof(mesh->coneSection, p, size);
919: return(0);
920: }
924: /*@
925: DMPlexSetConeSize - Set the number of in-edges for this point in the Sieve DAG
927: Not collective
929: Input Parameters:
930: + mesh - The DMPlex
931: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
932: - size - The cone size for point p
934: Output Parameter:
936: Note:
937: This should be called after DMPlexSetChart().
939: Level: beginner
941: .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart()
942: @*/
943: PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size)
944: {
945: DM_Plex *mesh = (DM_Plex*) dm->data;
950: PetscSectionSetDof(mesh->coneSection, p, size);
952: mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
953: return(0);
954: }
958: /*@
959: DMPlexAddConeSize - Add the given number of in-edges to this point in the Sieve DAG
961: Not collective
963: Input Parameters:
964: + mesh - The DMPlex
965: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
966: - size - The additional cone size for point p
968: Output Parameter:
970: Note:
971: This should be called after DMPlexSetChart().
973: Level: beginner
975: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart()
976: @*/
977: PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size)
978: {
979: DM_Plex *mesh = (DM_Plex*) dm->data;
980: PetscInt csize;
985: PetscSectionAddDof(mesh->coneSection, p, size);
986: PetscSectionGetDof(mesh->coneSection, p, &csize);
988: mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
989: return(0);
990: }
994: /*@C
995: DMPlexGetCone - Return the points on the in-edges for this point in the Sieve DAG
997: Not collective
999: Input Parameters:
1000: + mesh - The DMPlex
1001: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1003: Output Parameter:
1004: . cone - An array of points which are on the in-edges for point p
1006: Level: beginner
1008: Fortran Notes:
1009: Since it returns an array, this routine is only available in Fortran 90, and you must
1010: include petsc.h90 in your code.
1012: You must also call DMPlexRestoreCone() after you finish using the returned array.
1014: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart()
1015: @*/
1016: PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
1017: {
1018: DM_Plex *mesh = (DM_Plex*) dm->data;
1019: PetscInt off;
1025: PetscSectionGetOffset(mesh->coneSection, p, &off);
1026: *cone = &mesh->cones[off];
1027: return(0);
1028: }
1032: /*@
1033: DMPlexSetCone - Set the points on the in-edges for this point in the Sieve DAG
1035: Not collective
1037: Input Parameters:
1038: + mesh - The DMPlex
1039: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1040: - cone - An array of points which are on the in-edges for point p
1042: Output Parameter:
1044: Note:
1045: This should be called after all calls to DMPlexSetConeSize() and DMSetUp().
1047: Level: beginner
1049: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1050: @*/
1051: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
1052: {
1053: DM_Plex *mesh = (DM_Plex*) dm->data;
1054: PetscInt pStart, pEnd;
1055: PetscInt dof, off, c;
1060: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1061: PetscSectionGetDof(mesh->coneSection, p, &dof);
1063: PetscSectionGetOffset(mesh->coneSection, p, &off);
1064: 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);
1065: for (c = 0; c < dof; ++c) {
1066: 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);
1067: mesh->cones[off+c] = cone[c];
1068: }
1069: return(0);
1070: }
1074: /*@C
1075: DMPlexGetConeOrientation - Return the orientations on the in-edges for this point in the Sieve DAG
1077: Not collective
1079: Input Parameters:
1080: + mesh - The DMPlex
1081: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1083: Output Parameter:
1084: . coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1085: integer giving the prescription for cone traversal. If it is negative, the cone is
1086: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1087: the index of the cone point on which to start.
1089: Level: beginner
1091: Fortran Notes:
1092: Since it returns an array, this routine is only available in Fortran 90, and you must
1093: include petsc.h90 in your code.
1095: You must also call DMPlexRestoreConeOrientation() after you finish using the returned array.
1097: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
1098: @*/
1099: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
1100: {
1101: DM_Plex *mesh = (DM_Plex*) dm->data;
1102: PetscInt off;
1107: #if defined(PETSC_USE_DEBUG)
1108: {
1109: PetscInt dof;
1110: PetscSectionGetDof(mesh->coneSection, p, &dof);
1112: }
1113: #endif
1114: PetscSectionGetOffset(mesh->coneSection, p, &off);
1116: *coneOrientation = &mesh->coneOrientations[off];
1117: return(0);
1118: }
1122: /*@
1123: DMPlexSetConeOrientation - Set the orientations on the in-edges for this point in the Sieve DAG
1125: Not collective
1127: Input Parameters:
1128: + mesh - The DMPlex
1129: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1130: - coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1131: integer giving the prescription for cone traversal. If it is negative, the cone is
1132: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1133: the index of the cone point on which to start.
1135: Output Parameter:
1137: Note:
1138: This should be called after all calls to DMPlexSetConeSize() and DMSetUp().
1140: Level: beginner
1142: .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1143: @*/
1144: PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
1145: {
1146: DM_Plex *mesh = (DM_Plex*) dm->data;
1147: PetscInt pStart, pEnd;
1148: PetscInt dof, off, c;
1153: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1154: PetscSectionGetDof(mesh->coneSection, p, &dof);
1156: PetscSectionGetOffset(mesh->coneSection, p, &off);
1157: 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);
1158: for (c = 0; c < dof; ++c) {
1159: PetscInt cdof, o = coneOrientation[c];
1161: PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1162: 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);
1163: mesh->coneOrientations[off+c] = o;
1164: }
1165: return(0);
1166: }
1170: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
1171: {
1172: DM_Plex *mesh = (DM_Plex*) dm->data;
1173: PetscInt pStart, pEnd;
1174: PetscInt dof, off;
1179: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1180: 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);
1181: 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);
1182: PetscSectionGetDof(mesh->coneSection, p, &dof);
1183: PetscSectionGetOffset(mesh->coneSection, p, &off);
1184: 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);
1185: mesh->cones[off+conePos] = conePoint;
1186: return(0);
1187: }
1191: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1192: {
1193: DM_Plex *mesh = (DM_Plex*) dm->data;
1194: PetscInt pStart, pEnd;
1195: PetscInt dof, off;
1200: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1201: 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);
1202: PetscSectionGetDof(mesh->coneSection, p, &dof);
1203: PetscSectionGetOffset(mesh->coneSection, p, &off);
1204: 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);
1205: mesh->coneOrientations[off+conePos] = coneOrientation;
1206: return(0);
1207: }
1211: /*@
1212: DMPlexGetSupportSize - Return the number of out-edges for this point in the Sieve DAG
1214: Not collective
1216: Input Parameters:
1217: + mesh - The DMPlex
1218: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1220: Output Parameter:
1221: . size - The support size for point p
1223: Level: beginner
1225: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1226: @*/
1227: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1228: {
1229: DM_Plex *mesh = (DM_Plex*) dm->data;
1235: PetscSectionGetDof(mesh->supportSection, p, size);
1236: return(0);
1237: }
1241: /*@
1242: DMPlexSetSupportSize - Set the number of out-edges for this point in the Sieve DAG
1244: Not collective
1246: Input Parameters:
1247: + mesh - The DMPlex
1248: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1249: - size - The support size for point p
1251: Output Parameter:
1253: Note:
1254: This should be called after DMPlexSetChart().
1256: Level: beginner
1258: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1259: @*/
1260: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1261: {
1262: DM_Plex *mesh = (DM_Plex*) dm->data;
1267: PetscSectionSetDof(mesh->supportSection, p, size);
1269: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1270: return(0);
1271: }
1275: /*@C
1276: DMPlexGetSupport - Return the points on the out-edges for this point in the Sieve DAG
1278: Not collective
1280: Input Parameters:
1281: + mesh - The DMPlex
1282: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1284: Output Parameter:
1285: . support - An array of points which are on the out-edges for point p
1287: Level: beginner
1289: Fortran Notes:
1290: Since it returns an array, this routine is only available in Fortran 90, and you must
1291: include petsc.h90 in your code.
1293: You must also call DMPlexRestoreSupport() after you finish using the returned array.
1295: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1296: @*/
1297: PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1298: {
1299: DM_Plex *mesh = (DM_Plex*) dm->data;
1300: PetscInt off;
1306: PetscSectionGetOffset(mesh->supportSection, p, &off);
1307: *support = &mesh->supports[off];
1308: return(0);
1309: }
1313: /*@
1314: DMPlexSetSupport - Set the points on the out-edges for this point in the Sieve DAG
1316: Not collective
1318: Input Parameters:
1319: + mesh - The DMPlex
1320: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1321: - support - An array of points which are on the in-edges for point p
1323: Output Parameter:
1325: Note:
1326: This should be called after all calls to DMPlexSetSupportSize() and DMSetUp().
1328: Level: beginner
1330: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1331: @*/
1332: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1333: {
1334: DM_Plex *mesh = (DM_Plex*) dm->data;
1335: PetscInt pStart, pEnd;
1336: PetscInt dof, off, c;
1341: PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1342: PetscSectionGetDof(mesh->supportSection, p, &dof);
1344: PetscSectionGetOffset(mesh->supportSection, p, &off);
1345: 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);
1346: for (c = 0; c < dof; ++c) {
1347: 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);
1348: mesh->supports[off+c] = support[c];
1349: }
1350: return(0);
1351: }
1355: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1356: {
1357: DM_Plex *mesh = (DM_Plex*) dm->data;
1358: PetscInt pStart, pEnd;
1359: PetscInt dof, off;
1364: PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1365: PetscSectionGetDof(mesh->supportSection, p, &dof);
1366: PetscSectionGetOffset(mesh->supportSection, p, &off);
1367: 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);
1368: 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);
1369: 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);
1370: mesh->supports[off+supportPos] = supportPoint;
1371: return(0);
1372: }
1376: /*@C
1377: DMPlexGetTransitiveClosure - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1379: Not collective
1381: Input Parameters:
1382: + mesh - The DMPlex
1383: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1384: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1385: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
1387: Output Parameters:
1388: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1389: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
1391: Note:
1392: If using internal storage (points is NULL on input), each call overwrites the last output.
1394: Fortran Notes:
1395: Since it returns an array, this routine is only available in Fortran 90, and you must
1396: include petsc.h90 in your code.
1398: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1400: Level: beginner
1402: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1403: @*/
1404: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1405: {
1406: DM_Plex *mesh = (DM_Plex*) dm->data;
1407: PetscInt *closure, *fifo;
1408: const PetscInt *tmp = NULL, *tmpO = NULL;
1409: PetscInt tmpSize, t;
1410: PetscInt depth = 0, maxSize;
1411: PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0;
1412: PetscErrorCode ierr;
1416: DMPlexGetDepth(dm, &depth);
1417: /* This is only 1-level */
1418: if (useCone) {
1419: DMPlexGetConeSize(dm, p, &tmpSize);
1420: DMPlexGetCone(dm, p, &tmp);
1421: DMPlexGetConeOrientation(dm, p, &tmpO);
1422: } else {
1423: DMPlexGetSupportSize(dm, p, &tmpSize);
1424: DMPlexGetSupport(dm, p, &tmp);
1425: }
1426: if (depth == 1) {
1427: if (*points) {
1428: closure = *points;
1429: } else {
1430: maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1431: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1432: }
1433: closure[0] = p; closure[1] = 0;
1434: for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1435: closure[closureSize] = tmp[t];
1436: closure[closureSize+1] = tmpO ? tmpO[t] : 0;
1437: }
1438: if (numPoints) *numPoints = closureSize/2;
1439: if (points) *points = closure;
1440: return(0);
1441: }
1442: {
1443: PetscInt c, coneSeries, s,supportSeries;
1445: c = mesh->maxConeSize;
1446: coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1447: s = mesh->maxSupportSize;
1448: supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1449: maxSize = 2*PetscMax(coneSeries,supportSeries);
1450: }
1451: DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1452: if (*points) {
1453: closure = *points;
1454: } else {
1455: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1456: }
1457: closure[0] = p; closure[1] = 0;
1458: for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1459: const PetscInt cp = tmp[t];
1460: const PetscInt co = tmpO ? tmpO[t] : 0;
1462: closure[closureSize] = cp;
1463: closure[closureSize+1] = co;
1464: fifo[fifoSize] = cp;
1465: fifo[fifoSize+1] = co;
1466: }
1467: /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1468: while (fifoSize - fifoStart) {
1469: const PetscInt q = fifo[fifoStart];
1470: const PetscInt o = fifo[fifoStart+1];
1471: const PetscInt rev = o >= 0 ? 0 : 1;
1472: const PetscInt off = rev ? -(o+1) : o;
1474: if (useCone) {
1475: DMPlexGetConeSize(dm, q, &tmpSize);
1476: DMPlexGetCone(dm, q, &tmp);
1477: DMPlexGetConeOrientation(dm, q, &tmpO);
1478: } else {
1479: DMPlexGetSupportSize(dm, q, &tmpSize);
1480: DMPlexGetSupport(dm, q, &tmp);
1481: tmpO = NULL;
1482: }
1483: for (t = 0; t < tmpSize; ++t) {
1484: const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize;
1485: const PetscInt cp = tmp[i];
1486: /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1487: /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1488: const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1489: PetscInt co = tmpO ? tmpO[i] : 0;
1490: PetscInt c;
1492: if (rev) {
1493: PetscInt childSize, coff;
1494: DMPlexGetConeSize(dm, cp, &childSize);
1495: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1496: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1497: }
1498: /* Check for duplicate */
1499: for (c = 0; c < closureSize; c += 2) {
1500: if (closure[c] == cp) break;
1501: }
1502: if (c == closureSize) {
1503: closure[closureSize] = cp;
1504: closure[closureSize+1] = co;
1505: fifo[fifoSize] = cp;
1506: fifo[fifoSize+1] = co;
1507: closureSize += 2;
1508: fifoSize += 2;
1509: }
1510: }
1511: fifoStart += 2;
1512: }
1513: if (numPoints) *numPoints = closureSize/2;
1514: if (points) *points = closure;
1515: DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1516: return(0);
1517: }
1521: /*@C
1522: 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
1524: Not collective
1526: Input Parameters:
1527: + mesh - The DMPlex
1528: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1529: . orientation - The orientation of the point
1530: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1531: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
1533: Output Parameters:
1534: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1535: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
1537: Note:
1538: If using internal storage (points is NULL on input), each call overwrites the last output.
1540: Fortran Notes:
1541: Since it returns an array, this routine is only available in Fortran 90, and you must
1542: include petsc.h90 in your code.
1544: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1546: Level: beginner
1548: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1549: @*/
1550: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1551: {
1552: DM_Plex *mesh = (DM_Plex*) dm->data;
1553: PetscInt *closure, *fifo;
1554: const PetscInt *tmp = NULL, *tmpO = NULL;
1555: PetscInt tmpSize, t;
1556: PetscInt depth = 0, maxSize;
1557: PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0;
1558: PetscErrorCode ierr;
1562: DMPlexGetDepth(dm, &depth);
1563: /* This is only 1-level */
1564: if (useCone) {
1565: DMPlexGetConeSize(dm, p, &tmpSize);
1566: DMPlexGetCone(dm, p, &tmp);
1567: DMPlexGetConeOrientation(dm, p, &tmpO);
1568: } else {
1569: DMPlexGetSupportSize(dm, p, &tmpSize);
1570: DMPlexGetSupport(dm, p, &tmp);
1571: }
1572: if (depth == 1) {
1573: if (*points) {
1574: closure = *points;
1575: } else {
1576: maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1577: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1578: }
1579: closure[0] = p; closure[1] = ornt;
1580: for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1581: const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1582: closure[closureSize] = tmp[i];
1583: closure[closureSize+1] = tmpO ? tmpO[i] : 0;
1584: }
1585: if (numPoints) *numPoints = closureSize/2;
1586: if (points) *points = closure;
1587: return(0);
1588: }
1589: {
1590: PetscInt c, coneSeries, s,supportSeries;
1592: c = mesh->maxConeSize;
1593: coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1594: s = mesh->maxSupportSize;
1595: supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1596: maxSize = 2*PetscMax(coneSeries,supportSeries);
1597: }
1598: DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1599: if (*points) {
1600: closure = *points;
1601: } else {
1602: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1603: }
1604: closure[0] = p; closure[1] = ornt;
1605: for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1606: const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1607: const PetscInt cp = tmp[i];
1608: PetscInt co = tmpO ? tmpO[i] : 0;
1610: if (ornt < 0) {
1611: PetscInt childSize, coff;
1612: DMPlexGetConeSize(dm, cp, &childSize);
1613: coff = co < 0 ? -(tmpO[i]+1) : tmpO[i];
1614: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1615: }
1616: closure[closureSize] = cp;
1617: closure[closureSize+1] = co;
1618: fifo[fifoSize] = cp;
1619: fifo[fifoSize+1] = co;
1620: }
1621: /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1622: while (fifoSize - fifoStart) {
1623: const PetscInt q = fifo[fifoStart];
1624: const PetscInt o = fifo[fifoStart+1];
1625: const PetscInt rev = o >= 0 ? 0 : 1;
1626: const PetscInt off = rev ? -(o+1) : o;
1628: if (useCone) {
1629: DMPlexGetConeSize(dm, q, &tmpSize);
1630: DMPlexGetCone(dm, q, &tmp);
1631: DMPlexGetConeOrientation(dm, q, &tmpO);
1632: } else {
1633: DMPlexGetSupportSize(dm, q, &tmpSize);
1634: DMPlexGetSupport(dm, q, &tmp);
1635: tmpO = NULL;
1636: }
1637: for (t = 0; t < tmpSize; ++t) {
1638: const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize;
1639: const PetscInt cp = tmp[i];
1640: /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1641: /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1642: const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1643: PetscInt co = tmpO ? tmpO[i] : 0;
1644: PetscInt c;
1646: if (rev) {
1647: PetscInt childSize, coff;
1648: DMPlexGetConeSize(dm, cp, &childSize);
1649: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1650: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1651: }
1652: /* Check for duplicate */
1653: for (c = 0; c < closureSize; c += 2) {
1654: if (closure[c] == cp) break;
1655: }
1656: if (c == closureSize) {
1657: closure[closureSize] = cp;
1658: closure[closureSize+1] = co;
1659: fifo[fifoSize] = cp;
1660: fifo[fifoSize+1] = co;
1661: closureSize += 2;
1662: fifoSize += 2;
1663: }
1664: }
1665: fifoStart += 2;
1666: }
1667: if (numPoints) *numPoints = closureSize/2;
1668: if (points) *points = closure;
1669: DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1670: return(0);
1671: }
1675: /*@C
1676: DMPlexRestoreTransitiveClosure - Restore the array of points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1678: Not collective
1680: Input Parameters:
1681: + mesh - The DMPlex
1682: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1683: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1684: . numPoints - The number of points in the closure, so points[] is of size 2*numPoints, zeroed on exit
1685: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...], zeroed on exit
1687: Note:
1688: If not using internal storage (points is not NULL on input), this call is unnecessary
1690: Fortran Notes:
1691: Since it returns an array, this routine is only available in Fortran 90, and you must
1692: include petsc.h90 in your code.
1694: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1696: Level: beginner
1698: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1699: @*/
1700: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1701: {
1708: DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1709: if (numPoints) *numPoints = 0;
1710: return(0);
1711: }
1715: /*@
1716: DMPlexGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG
1718: Not collective
1720: Input Parameter:
1721: . mesh - The DMPlex
1723: Output Parameters:
1724: + maxConeSize - The maximum number of in-edges
1725: - maxSupportSize - The maximum number of out-edges
1727: Level: beginner
1729: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1730: @*/
1731: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1732: {
1733: DM_Plex *mesh = (DM_Plex*) dm->data;
1737: if (maxConeSize) *maxConeSize = mesh->maxConeSize;
1738: if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1739: return(0);
1740: }
1744: PetscErrorCode DMSetUp_Plex(DM dm)
1745: {
1746: DM_Plex *mesh = (DM_Plex*) dm->data;
1747: PetscInt size;
1752: PetscSectionSetUp(mesh->coneSection);
1753: PetscSectionGetStorageSize(mesh->coneSection, &size);
1754: PetscMalloc1(size, &mesh->cones);
1755: PetscCalloc1(size, &mesh->coneOrientations);
1756: if (mesh->maxSupportSize) {
1757: PetscSectionSetUp(mesh->supportSection);
1758: PetscSectionGetStorageSize(mesh->supportSection, &size);
1759: PetscMalloc1(size, &mesh->supports);
1760: }
1761: return(0);
1762: }
1766: PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1767: {
1771: if (subdm) {DMClone(dm, subdm);}
1772: DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1773: return(0);
1774: }
1778: /*@
1779: DMPlexSymmetrize - Creates support (out-edge) information from cone (in-edge) inoformation
1781: Not collective
1783: Input Parameter:
1784: . mesh - The DMPlex
1786: Output Parameter:
1788: Note:
1789: This should be called after all calls to DMPlexSetCone()
1791: Level: beginner
1793: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1794: @*/
1795: PetscErrorCode DMPlexSymmetrize(DM dm)
1796: {
1797: DM_Plex *mesh = (DM_Plex*) dm->data;
1798: PetscInt *offsets;
1799: PetscInt supportSize;
1800: PetscInt pStart, pEnd, p;
1805: if (mesh->supports) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Supports were already setup in this DMPlex");
1806: /* Calculate support sizes */
1807: DMPlexGetChart(dm, &pStart, &pEnd);
1808: for (p = pStart; p < pEnd; ++p) {
1809: PetscInt dof, off, c;
1811: PetscSectionGetDof(mesh->coneSection, p, &dof);
1812: PetscSectionGetOffset(mesh->coneSection, p, &off);
1813: for (c = off; c < off+dof; ++c) {
1814: PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1815: }
1816: }
1817: for (p = pStart; p < pEnd; ++p) {
1818: PetscInt dof;
1820: PetscSectionGetDof(mesh->supportSection, p, &dof);
1822: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1823: }
1824: PetscSectionSetUp(mesh->supportSection);
1825: /* Calculate supports */
1826: PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1827: PetscMalloc1(supportSize, &mesh->supports);
1828: PetscCalloc1(pEnd - pStart, &offsets);
1829: for (p = pStart; p < pEnd; ++p) {
1830: PetscInt dof, off, c;
1832: PetscSectionGetDof(mesh->coneSection, p, &dof);
1833: PetscSectionGetOffset(mesh->coneSection, p, &off);
1834: for (c = off; c < off+dof; ++c) {
1835: const PetscInt q = mesh->cones[c];
1836: PetscInt offS;
1838: PetscSectionGetOffset(mesh->supportSection, q, &offS);
1840: mesh->supports[offS+offsets[q]] = p;
1841: ++offsets[q];
1842: }
1843: }
1844: PetscFree(offsets);
1845: return(0);
1846: }
1850: /*@
1851: DMPlexStratify - The Sieve DAG for most topologies is a graded poset (http://en.wikipedia.org/wiki/Graded_poset), and
1852: can be illustrated by Hasse Diagram (a http://en.wikipedia.org/wiki/Hasse_diagram). The strata group all points of the
1853: same grade, and this function calculates the strata. This grade can be seen as the height (or depth) of the point in
1854: the DAG.
1856: Collective on dm
1858: Input Parameter:
1859: . mesh - The DMPlex
1861: Output Parameter:
1863: Notes:
1864: Concretely, DMPlexStratify() creates a new label named "depth" containing the dimension of each element: 0 for vertices,
1865: 1 for edges, and so on. The depth label can be accessed through DMPlexGetDepthLabel() or DMPlexGetDepthStratum(), or
1866: manually via DMGetLabel(). The height is defined implicitly by height = maxDimension - depth, and can be accessed
1867: via DMPlexGetHeightStratum(). For example, cells have height 0 and faces have height 1.
1869: DMPlexStratify() should be called after all calls to DMPlexSymmetrize()
1871: Level: beginner
1873: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1874: @*/
1875: PetscErrorCode DMPlexStratify(DM dm)
1876: {
1877: DM_Plex *mesh = (DM_Plex*) dm->data;
1878: DMLabel label;
1879: PetscInt pStart, pEnd, p;
1880: PetscInt numRoots = 0, numLeaves = 0;
1885: PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1886: /* Calculate depth */
1887: DMPlexGetChart(dm, &pStart, &pEnd);
1888: DMCreateLabel(dm, "depth");
1889: DMPlexGetDepthLabel(dm, &label);
1890: /* Initialize roots and count leaves */
1891: for (p = pStart; p < pEnd; ++p) {
1892: PetscInt coneSize, supportSize;
1894: DMPlexGetConeSize(dm, p, &coneSize);
1895: DMPlexGetSupportSize(dm, p, &supportSize);
1896: if (!coneSize && supportSize) {
1897: ++numRoots;
1898: DMLabelSetValue(label, p, 0);
1899: } else if (!supportSize && coneSize) {
1900: ++numLeaves;
1901: } else if (!supportSize && !coneSize) {
1902: /* Isolated points */
1903: DMLabelSetValue(label, p, 0);
1904: }
1905: }
1906: if (numRoots + numLeaves == (pEnd - pStart)) {
1907: for (p = pStart; p < pEnd; ++p) {
1908: PetscInt coneSize, supportSize;
1910: DMPlexGetConeSize(dm, p, &coneSize);
1911: DMPlexGetSupportSize(dm, p, &supportSize);
1912: if (!supportSize && coneSize) {
1913: DMLabelSetValue(label, p, 1);
1914: }
1915: }
1916: } else {
1917: IS pointIS;
1918: PetscInt numPoints = 0, level = 0;
1920: DMLabelGetStratumIS(label, level, &pointIS);
1921: if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1922: while (numPoints) {
1923: const PetscInt *points;
1924: const PetscInt newLevel = level+1;
1926: ISGetIndices(pointIS, &points);
1927: for (p = 0; p < numPoints; ++p) {
1928: const PetscInt point = points[p];
1929: const PetscInt *support;
1930: PetscInt supportSize, s;
1932: DMPlexGetSupportSize(dm, point, &supportSize);
1933: DMPlexGetSupport(dm, point, &support);
1934: for (s = 0; s < supportSize; ++s) {
1935: DMLabelSetValue(label, support[s], newLevel);
1936: }
1937: }
1938: ISRestoreIndices(pointIS, &points);
1939: ++level;
1940: ISDestroy(&pointIS);
1941: DMLabelGetStratumIS(label, level, &pointIS);
1942: if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1943: else {numPoints = 0;}
1944: }
1945: ISDestroy(&pointIS);
1946: }
1947: { /* just in case there is an empty process */
1948: PetscInt numValues, maxValues = 0, v;
1950: DMLabelGetNumValues(label,&numValues);
1951: MPI_Allreduce(&numValues,&maxValues,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
1952: for (v = numValues; v < maxValues; v++) {
1953: DMLabelAddStratum(label,v);
1954: }
1955: }
1957: DMLabelGetState(label, &mesh->depthState);
1958: PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1959: return(0);
1960: }
1964: /*@C
1965: DMPlexGetJoin - Get an array for the join of the set of points
1967: Not Collective
1969: Input Parameters:
1970: + dm - The DMPlex object
1971: . numPoints - The number of input points for the join
1972: - points - The input points
1974: Output Parameters:
1975: + numCoveredPoints - The number of points in the join
1976: - coveredPoints - The points in the join
1978: Level: intermediate
1980: Note: Currently, this is restricted to a single level join
1982: Fortran Notes:
1983: Since it returns an array, this routine is only available in Fortran 90, and you must
1984: include petsc.h90 in your code.
1986: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1988: .keywords: mesh
1989: .seealso: DMPlexRestoreJoin(), DMPlexGetMeet()
1990: @*/
1991: PetscErrorCode DMPlexGetJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1992: {
1993: DM_Plex *mesh = (DM_Plex*) dm->data;
1994: PetscInt *join[2];
1995: PetscInt joinSize, i = 0;
1996: PetscInt dof, off, p, c, m;
2004: DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[0]);
2005: DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1]);
2006: /* Copy in support of first point */
2007: PetscSectionGetDof(mesh->supportSection, points[0], &dof);
2008: PetscSectionGetOffset(mesh->supportSection, points[0], &off);
2009: for (joinSize = 0; joinSize < dof; ++joinSize) {
2010: join[i][joinSize] = mesh->supports[off+joinSize];
2011: }
2012: /* Check each successive support */
2013: for (p = 1; p < numPoints; ++p) {
2014: PetscInt newJoinSize = 0;
2016: PetscSectionGetDof(mesh->supportSection, points[p], &dof);
2017: PetscSectionGetOffset(mesh->supportSection, points[p], &off);
2018: for (c = 0; c < dof; ++c) {
2019: const PetscInt point = mesh->supports[off+c];
2021: for (m = 0; m < joinSize; ++m) {
2022: if (point == join[i][m]) {
2023: join[1-i][newJoinSize++] = point;
2024: break;
2025: }
2026: }
2027: }
2028: joinSize = newJoinSize;
2029: i = 1-i;
2030: }
2031: *numCoveredPoints = joinSize;
2032: *coveredPoints = join[i];
2033: DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2034: return(0);
2035: }
2039: /*@C
2040: DMPlexRestoreJoin - Restore an array for the join of the set of points
2042: Not Collective
2044: Input Parameters:
2045: + dm - The DMPlex object
2046: . numPoints - The number of input points for the join
2047: - points - The input points
2049: Output Parameters:
2050: + numCoveredPoints - The number of points in the join
2051: - coveredPoints - The points in the join
2053: Fortran Notes:
2054: Since it returns an array, this routine is only available in Fortran 90, and you must
2055: include petsc.h90 in your code.
2057: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2059: Level: intermediate
2061: .keywords: mesh
2062: .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet()
2063: @*/
2064: PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2065: {
2073: DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2074: if (numCoveredPoints) *numCoveredPoints = 0;
2075: return(0);
2076: }
2080: /*@C
2081: DMPlexGetFullJoin - Get an array for the join of the set of points
2083: Not Collective
2085: Input Parameters:
2086: + dm - The DMPlex object
2087: . numPoints - The number of input points for the join
2088: - points - The input points
2090: Output Parameters:
2091: + numCoveredPoints - The number of points in the join
2092: - coveredPoints - The points in the join
2094: Fortran Notes:
2095: Since it returns an array, this routine is only available in Fortran 90, and you must
2096: include petsc.h90 in your code.
2098: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2100: Level: intermediate
2102: .keywords: mesh
2103: .seealso: DMPlexGetJoin(), DMPlexRestoreJoin(), DMPlexGetMeet()
2104: @*/
2105: PetscErrorCode DMPlexGetFullJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2106: {
2107: DM_Plex *mesh = (DM_Plex*) dm->data;
2108: PetscInt *offsets, **closures;
2109: PetscInt *join[2];
2110: PetscInt depth = 0, maxSize, joinSize = 0, i = 0;
2111: PetscInt p, d, c, m, ms;
2120: DMPlexGetDepth(dm, &depth);
2121: PetscCalloc1(numPoints, &closures);
2122: DMGetWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2123: ms = mesh->maxSupportSize;
2124: maxSize = (ms > 1) ? ((PetscPowInt(ms,depth+1)-1)/(ms-1)) : depth + 1;
2125: DMGetWorkArray(dm, maxSize, PETSC_INT, &join[0]);
2126: DMGetWorkArray(dm, maxSize, PETSC_INT, &join[1]);
2128: for (p = 0; p < numPoints; ++p) {
2129: PetscInt closureSize;
2131: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closures[p]);
2133: offsets[p*(depth+2)+0] = 0;
2134: for (d = 0; d < depth+1; ++d) {
2135: PetscInt pStart, pEnd, i;
2137: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
2138: for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) {
2139: if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2140: offsets[p*(depth+2)+d+1] = i;
2141: break;
2142: }
2143: }
2144: if (i == closureSize) offsets[p*(depth+2)+d+1] = i;
2145: }
2146: 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);
2147: }
2148: for (d = 0; d < depth+1; ++d) {
2149: PetscInt dof;
2151: /* Copy in support of first point */
2152: dof = offsets[d+1] - offsets[d];
2153: for (joinSize = 0; joinSize < dof; ++joinSize) {
2154: join[i][joinSize] = closures[0][(offsets[d]+joinSize)*2];
2155: }
2156: /* Check each successive cone */
2157: for (p = 1; p < numPoints && joinSize; ++p) {
2158: PetscInt newJoinSize = 0;
2160: dof = offsets[p*(depth+2)+d+1] - offsets[p*(depth+2)+d];
2161: for (c = 0; c < dof; ++c) {
2162: const PetscInt point = closures[p][(offsets[p*(depth+2)+d]+c)*2];
2164: for (m = 0; m < joinSize; ++m) {
2165: if (point == join[i][m]) {
2166: join[1-i][newJoinSize++] = point;
2167: break;
2168: }
2169: }
2170: }
2171: joinSize = newJoinSize;
2172: i = 1-i;
2173: }
2174: if (joinSize) break;
2175: }
2176: *numCoveredPoints = joinSize;
2177: *coveredPoints = join[i];
2178: for (p = 0; p < numPoints; ++p) {
2179: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, NULL, &closures[p]);
2180: }
2181: PetscFree(closures);
2182: DMRestoreWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2183: DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2184: return(0);
2185: }
2189: /*@C
2190: DMPlexGetMeet - Get an array for the meet of the set of points
2192: Not Collective
2194: Input Parameters:
2195: + dm - The DMPlex object
2196: . numPoints - The number of input points for the meet
2197: - points - The input points
2199: Output Parameters:
2200: + numCoveredPoints - The number of points in the meet
2201: - coveredPoints - The points in the meet
2203: Level: intermediate
2205: Note: Currently, this is restricted to a single level meet
2207: Fortran Notes:
2208: Since it returns an array, this routine is only available in Fortran 90, and you must
2209: include petsc.h90 in your code.
2211: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2213: .keywords: mesh
2214: .seealso: DMPlexRestoreMeet(), DMPlexGetJoin()
2215: @*/
2216: PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2217: {
2218: DM_Plex *mesh = (DM_Plex*) dm->data;
2219: PetscInt *meet[2];
2220: PetscInt meetSize, i = 0;
2221: PetscInt dof, off, p, c, m;
2229: DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[0]);
2230: DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1]);
2231: /* Copy in cone of first point */
2232: PetscSectionGetDof(mesh->coneSection, points[0], &dof);
2233: PetscSectionGetOffset(mesh->coneSection, points[0], &off);
2234: for (meetSize = 0; meetSize < dof; ++meetSize) {
2235: meet[i][meetSize] = mesh->cones[off+meetSize];
2236: }
2237: /* Check each successive cone */
2238: for (p = 1; p < numPoints; ++p) {
2239: PetscInt newMeetSize = 0;
2241: PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2242: PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2243: for (c = 0; c < dof; ++c) {
2244: const PetscInt point = mesh->cones[off+c];
2246: for (m = 0; m < meetSize; ++m) {
2247: if (point == meet[i][m]) {
2248: meet[1-i][newMeetSize++] = point;
2249: break;
2250: }
2251: }
2252: }
2253: meetSize = newMeetSize;
2254: i = 1-i;
2255: }
2256: *numCoveringPoints = meetSize;
2257: *coveringPoints = meet[i];
2258: DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2259: return(0);
2260: }
2264: /*@C
2265: DMPlexRestoreMeet - Restore an array for the meet of the set of points
2267: Not Collective
2269: Input Parameters:
2270: + dm - The DMPlex object
2271: . numPoints - The number of input points for the meet
2272: - points - The input points
2274: Output Parameters:
2275: + numCoveredPoints - The number of points in the meet
2276: - coveredPoints - The points in the meet
2278: Level: intermediate
2280: Fortran Notes:
2281: Since it returns an array, this routine is only available in Fortran 90, and you must
2282: include petsc.h90 in your code.
2284: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2286: .keywords: mesh
2287: .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin()
2288: @*/
2289: PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2290: {
2298: DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2299: if (numCoveredPoints) *numCoveredPoints = 0;
2300: return(0);
2301: }
2305: /*@C
2306: DMPlexGetFullMeet - Get an array for the meet of the set of points
2308: Not Collective
2310: Input Parameters:
2311: + dm - The DMPlex object
2312: . numPoints - The number of input points for the meet
2313: - points - The input points
2315: Output Parameters:
2316: + numCoveredPoints - The number of points in the meet
2317: - coveredPoints - The points in the meet
2319: Level: intermediate
2321: Fortran Notes:
2322: Since it returns an array, this routine is only available in Fortran 90, and you must
2323: include petsc.h90 in your code.
2325: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2327: .keywords: mesh
2328: .seealso: DMPlexGetMeet(), DMPlexRestoreMeet(), DMPlexGetJoin()
2329: @*/
2330: PetscErrorCode DMPlexGetFullMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2331: {
2332: DM_Plex *mesh = (DM_Plex*) dm->data;
2333: PetscInt *offsets, **closures;
2334: PetscInt *meet[2];
2335: PetscInt height = 0, maxSize, meetSize = 0, i = 0;
2336: PetscInt p, h, c, m, mc;
2345: DMPlexGetDepth(dm, &height);
2346: PetscMalloc1(numPoints, &closures);
2347: DMGetWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2348: mc = mesh->maxConeSize;
2349: maxSize = (mc > 1) ? ((PetscPowInt(mc,height+1)-1)/(mc-1)) : height + 1;
2350: DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[0]);
2351: DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[1]);
2353: for (p = 0; p < numPoints; ++p) {
2354: PetscInt closureSize;
2356: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closures[p]);
2358: offsets[p*(height+2)+0] = 0;
2359: for (h = 0; h < height+1; ++h) {
2360: PetscInt pStart, pEnd, i;
2362: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
2363: for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) {
2364: if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2365: offsets[p*(height+2)+h+1] = i;
2366: break;
2367: }
2368: }
2369: if (i == closureSize) offsets[p*(height+2)+h+1] = i;
2370: }
2371: 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);
2372: }
2373: for (h = 0; h < height+1; ++h) {
2374: PetscInt dof;
2376: /* Copy in cone of first point */
2377: dof = offsets[h+1] - offsets[h];
2378: for (meetSize = 0; meetSize < dof; ++meetSize) {
2379: meet[i][meetSize] = closures[0][(offsets[h]+meetSize)*2];
2380: }
2381: /* Check each successive cone */
2382: for (p = 1; p < numPoints && meetSize; ++p) {
2383: PetscInt newMeetSize = 0;
2385: dof = offsets[p*(height+2)+h+1] - offsets[p*(height+2)+h];
2386: for (c = 0; c < dof; ++c) {
2387: const PetscInt point = closures[p][(offsets[p*(height+2)+h]+c)*2];
2389: for (m = 0; m < meetSize; ++m) {
2390: if (point == meet[i][m]) {
2391: meet[1-i][newMeetSize++] = point;
2392: break;
2393: }
2394: }
2395: }
2396: meetSize = newMeetSize;
2397: i = 1-i;
2398: }
2399: if (meetSize) break;
2400: }
2401: *numCoveredPoints = meetSize;
2402: *coveredPoints = meet[i];
2403: for (p = 0; p < numPoints; ++p) {
2404: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, NULL, &closures[p]);
2405: }
2406: PetscFree(closures);
2407: DMRestoreWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2408: DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2409: return(0);
2410: }
2414: /*@C
2415: DMPlexEqual - Determine if two DMs have the same topology
2417: Not Collective
2419: Input Parameters:
2420: + dmA - A DMPlex object
2421: - dmB - A DMPlex object
2423: Output Parameters:
2424: . equal - PETSC_TRUE if the topologies are identical
2426: Level: intermediate
2428: Notes:
2429: We are not solving graph isomorphism, so we do not permutation.
2431: .keywords: mesh
2432: .seealso: DMPlexGetCone()
2433: @*/
2434: PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal)
2435: {
2436: PetscInt depth, depthB, pStart, pEnd, pStartB, pEndB, p;
2444: *equal = PETSC_FALSE;
2445: DMPlexGetDepth(dmA, &depth);
2446: DMPlexGetDepth(dmB, &depthB);
2447: if (depth != depthB) return(0);
2448: DMPlexGetChart(dmA, &pStart, &pEnd);
2449: DMPlexGetChart(dmB, &pStartB, &pEndB);
2450: if ((pStart != pStartB) || (pEnd != pEndB)) return(0);
2451: for (p = pStart; p < pEnd; ++p) {
2452: const PetscInt *cone, *coneB, *ornt, *orntB, *support, *supportB;
2453: PetscInt coneSize, coneSizeB, c, supportSize, supportSizeB, s;
2455: DMPlexGetConeSize(dmA, p, &coneSize);
2456: DMPlexGetCone(dmA, p, &cone);
2457: DMPlexGetConeOrientation(dmA, p, &ornt);
2458: DMPlexGetConeSize(dmB, p, &coneSizeB);
2459: DMPlexGetCone(dmB, p, &coneB);
2460: DMPlexGetConeOrientation(dmB, p, &orntB);
2461: if (coneSize != coneSizeB) return(0);
2462: for (c = 0; c < coneSize; ++c) {
2463: if (cone[c] != coneB[c]) return(0);
2464: if (ornt[c] != orntB[c]) return(0);
2465: }
2466: DMPlexGetSupportSize(dmA, p, &supportSize);
2467: DMPlexGetSupport(dmA, p, &support);
2468: DMPlexGetSupportSize(dmB, p, &supportSizeB);
2469: DMPlexGetSupport(dmB, p, &supportB);
2470: if (supportSize != supportSizeB) return(0);
2471: for (s = 0; s < supportSize; ++s) {
2472: if (support[s] != supportB[s]) return(0);
2473: }
2474: }
2475: *equal = PETSC_TRUE;
2476: return(0);
2477: }
2481: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2482: {
2483: MPI_Comm comm;
2487: PetscObjectGetComm((PetscObject)dm,&comm);
2489: switch (cellDim) {
2490: case 0:
2491: *numFaceVertices = 0;
2492: break;
2493: case 1:
2494: *numFaceVertices = 1;
2495: break;
2496: case 2:
2497: switch (numCorners) {
2498: case 3: /* triangle */
2499: *numFaceVertices = 2; /* Edge has 2 vertices */
2500: break;
2501: case 4: /* quadrilateral */
2502: *numFaceVertices = 2; /* Edge has 2 vertices */
2503: break;
2504: case 6: /* quadratic triangle, tri and quad cohesive Lagrange cells */
2505: *numFaceVertices = 3; /* Edge has 3 vertices */
2506: break;
2507: case 9: /* quadratic quadrilateral, quadratic quad cohesive Lagrange cells */
2508: *numFaceVertices = 3; /* Edge has 3 vertices */
2509: break;
2510: default:
2511: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %D for dimension %D", numCorners, cellDim);
2512: }
2513: break;
2514: case 3:
2515: switch (numCorners) {
2516: case 4: /* tetradehdron */
2517: *numFaceVertices = 3; /* Face has 3 vertices */
2518: break;
2519: case 6: /* tet cohesive cells */
2520: *numFaceVertices = 4; /* Face has 4 vertices */
2521: break;
2522: case 8: /* hexahedron */
2523: *numFaceVertices = 4; /* Face has 4 vertices */
2524: break;
2525: case 9: /* tet cohesive Lagrange cells */
2526: *numFaceVertices = 6; /* Face has 6 vertices */
2527: break;
2528: case 10: /* quadratic tetrahedron */
2529: *numFaceVertices = 6; /* Face has 6 vertices */
2530: break;
2531: case 12: /* hex cohesive Lagrange cells */
2532: *numFaceVertices = 6; /* Face has 6 vertices */
2533: break;
2534: case 18: /* quadratic tet cohesive Lagrange cells */
2535: *numFaceVertices = 6; /* Face has 6 vertices */
2536: break;
2537: case 27: /* quadratic hexahedron, quadratic hex cohesive Lagrange cells */
2538: *numFaceVertices = 9; /* Face has 9 vertices */
2539: break;
2540: default:
2541: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %D for dimension %D", numCorners, cellDim);
2542: }
2543: break;
2544: default:
2545: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid cell dimension %D", cellDim);
2546: }
2547: return(0);
2548: }
2552: /*@
2553: DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point
2555: Not Collective
2557: Input Parameter:
2558: . dm - The DMPlex object
2560: Output Parameter:
2561: . depthLabel - The DMLabel recording point depth
2563: Level: developer
2565: .keywords: mesh, points
2566: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2567: @*/
2568: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
2569: {
2575: if (!dm->depthLabel) {DMGetLabel(dm, "depth", &dm->depthLabel);}
2576: *depthLabel = dm->depthLabel;
2577: return(0);
2578: }
2582: /*@
2583: DMPlexGetDepth - Get the depth of the DAG representing this mesh
2585: Not Collective
2587: Input Parameter:
2588: . dm - The DMPlex object
2590: Output Parameter:
2591: . depth - The number of strata (breadth first levels) in the DAG
2593: Level: developer
2595: .keywords: mesh, points
2596: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2597: @*/
2598: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
2599: {
2600: DMLabel label;
2601: PetscInt d = 0;
2607: DMPlexGetDepthLabel(dm, &label);
2608: if (label) {DMLabelGetNumValues(label, &d);}
2609: *depth = d-1;
2610: return(0);
2611: }
2615: /*@
2616: DMPlexGetDepthStratum - Get the bounds [start, end) for all points at a certain depth.
2618: Not Collective
2620: Input Parameters:
2621: + dm - The DMPlex object
2622: - stratumValue - The requested depth
2624: Output Parameters:
2625: + start - The first point at this depth
2626: - end - One beyond the last point at this depth
2628: Level: developer
2630: .keywords: mesh, points
2631: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
2632: @*/
2633: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2634: {
2635: DMLabel label;
2636: PetscInt pStart, pEnd;
2643: DMPlexGetChart(dm, &pStart, &pEnd);
2644: if (pStart == pEnd) return(0);
2645: if (stratumValue < 0) {
2646: if (start) *start = pStart;
2647: if (end) *end = pEnd;
2648: return(0);
2649: }
2650: DMPlexGetDepthLabel(dm, &label);
2651: if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2652: DMLabelGetStratumBounds(label, stratumValue, start, end);
2653: return(0);
2654: }
2658: /*@
2659: DMPlexGetHeightStratum - Get the bounds [start, end) for all points at a certain height.
2661: Not Collective
2663: Input Parameters:
2664: + dm - The DMPlex object
2665: - stratumValue - The requested height
2667: Output Parameters:
2668: + start - The first point at this height
2669: - end - One beyond the last point at this height
2671: Level: developer
2673: .keywords: mesh, points
2674: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
2675: @*/
2676: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2677: {
2678: DMLabel label;
2679: PetscInt depth, pStart, pEnd;
2686: DMPlexGetChart(dm, &pStart, &pEnd);
2687: if (pStart == pEnd) return(0);
2688: if (stratumValue < 0) {
2689: if (start) *start = pStart;
2690: if (end) *end = pEnd;
2691: return(0);
2692: }
2693: DMPlexGetDepthLabel(dm, &label);
2694: if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2695: DMLabelGetNumValues(label, &depth);
2696: DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
2697: return(0);
2698: }
2702: /* Set the number of dof on each point and separate by fields */
2703: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
2704: {
2705: PetscInt *pMax;
2706: PetscInt depth, pStart = 0, pEnd = 0;
2707: PetscInt Nf, p, d, dep, f;
2708: PetscBool *isFE;
2712: PetscMalloc1(numFields, &isFE);
2713: DMGetNumFields(dm, &Nf);
2714: for (f = 0; f < numFields; ++f) {
2715: PetscObject obj;
2716: PetscClassId id;
2718: isFE[f] = PETSC_FALSE;
2719: if (f >= Nf) continue;
2720: DMGetField(dm, f, &obj);
2721: PetscObjectGetClassId(obj, &id);
2722: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
2723: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
2724: }
2725: PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
2726: if (numFields > 0) {
2727: PetscSectionSetNumFields(*section, numFields);
2728: if (numComp) {
2729: for (f = 0; f < numFields; ++f) {
2730: PetscSectionSetFieldComponents(*section, f, numComp[f]);
2731: }
2732: }
2733: }
2734: DMPlexGetChart(dm, &pStart, &pEnd);
2735: PetscSectionSetChart(*section, pStart, pEnd);
2736: DMPlexGetDepth(dm, &depth);
2737: PetscMalloc1(depth+1,&pMax);
2738: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2739: for (dep = 0; dep <= depth; ++dep) {
2740: d = dim == depth ? dep : (!dep ? 0 : dim);
2741: DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
2742: pMax[dep] = pMax[dep] < 0 ? pEnd : pMax[dep];
2743: for (p = pStart; p < pEnd; ++p) {
2744: PetscInt tot = 0;
2746: for (f = 0; f < numFields; ++f) {
2747: if (isFE[f] && p >= pMax[dep]) continue;
2748: PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
2749: tot += numDof[f*(dim+1)+d];
2750: }
2751: PetscSectionSetDof(*section, p, tot);
2752: }
2753: }
2754: PetscFree(pMax);
2755: PetscFree(isFE);
2756: return(0);
2757: }
2761: /* Set the number of dof on each point and separate by fields
2762: If bcComps is NULL or the IS is NULL, constrain every dof on the point
2763: */
2764: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2765: {
2766: PetscInt numFields;
2767: PetscInt bc;
2768: PetscSection aSec;
2772: PetscSectionGetNumFields(section, &numFields);
2773: for (bc = 0; bc < numBC; ++bc) {
2774: PetscInt field = 0;
2775: const PetscInt *comp;
2776: const PetscInt *idx;
2777: PetscInt Nc = -1, n, i;
2779: if (numFields) field = bcField[bc];
2780: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2781: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2782: ISGetLocalSize(bcPoints[bc], &n);
2783: ISGetIndices(bcPoints[bc], &idx);
2784: for (i = 0; i < n; ++i) {
2785: const PetscInt p = idx[i];
2786: PetscInt numConst;
2788: if (numFields) {
2789: PetscSectionGetFieldDof(section, p, field, &numConst);
2790: } else {
2791: PetscSectionGetDof(section, p, &numConst);
2792: }
2793: /* If Nc < 0, constrain every dof on the point */
2794: if (Nc > 0) numConst = PetscMin(numConst, Nc);
2795: if (numFields) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
2796: PetscSectionAddConstraintDof(section, p, numConst);
2797: }
2798: ISRestoreIndices(bcPoints[bc], &idx);
2799: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2800: }
2801: DMPlexGetAnchors(dm, &aSec, NULL);
2802: if (aSec) {
2803: PetscInt aStart, aEnd, a;
2805: PetscSectionGetChart(aSec, &aStart, &aEnd);
2806: for (a = aStart; a < aEnd; a++) {
2807: PetscInt dof, f;
2809: PetscSectionGetDof(aSec, a, &dof);
2810: if (dof) {
2811: /* if there are point-to-point constraints, then all dofs are constrained */
2812: PetscSectionGetDof(section, a, &dof);
2813: PetscSectionSetConstraintDof(section, a, dof);
2814: for (f = 0; f < numFields; f++) {
2815: PetscSectionGetFieldDof(section, a, f, &dof);
2816: PetscSectionSetFieldConstraintDof(section, a, f, dof);
2817: }
2818: }
2819: }
2820: }
2821: return(0);
2822: }
2826: /* Set the constrained field indices on each point
2827: If bcComps is NULL or the IS is NULL, constrain every dof on the point
2828: */
2829: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2830: {
2831: PetscSection aSec;
2832: PetscInt *indices;
2833: PetscInt numFields, maxDof, pStart, pEnd, p, bc, f, d;
2837: PetscSectionGetNumFields(section, &numFields);
2838: if (!numFields) return(0);
2839: /* Initialize all field indices to -1 */
2840: PetscSectionGetChart(section, &pStart, &pEnd);
2841: PetscSectionGetMaxDof(section, &maxDof);
2842: PetscMalloc1(maxDof, &indices);
2843: for (d = 0; d < maxDof; ++d) indices[d] = -1;
2844: for (p = pStart; p < pEnd; ++p) for (f = 0; f < numFields; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
2845: /* Handle BC constraints */
2846: for (bc = 0; bc < numBC; ++bc) {
2847: const PetscInt field = bcField[bc];
2848: const PetscInt *comp, *idx;
2849: PetscInt Nc = -1, n, i;
2851: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2852: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2853: ISGetLocalSize(bcPoints[bc], &n);
2854: ISGetIndices(bcPoints[bc], &idx);
2855: for (i = 0; i < n; ++i) {
2856: const PetscInt p = idx[i];
2857: const PetscInt *find;
2858: PetscInt fcdof, c;
2860: PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
2861: if (Nc < 0) {
2862: for (d = 0; d < fcdof; ++d) indices[d] = d;
2863: } else {
2864: PetscSectionGetFieldConstraintIndices(section, p, field, &find);
2865: for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
2866: for (c = 0; c < Nc; ++c) indices[d+c] = comp[c];
2867: PetscSortInt(d+Nc, indices);
2868: for (c = d+Nc; c < fcdof; ++c) indices[c] = -1;
2869: }
2870: PetscSectionSetFieldConstraintIndices(section, p, field, indices);
2871: }
2872: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2873: ISRestoreIndices(bcPoints[bc], &idx);
2874: }
2875: /* Handle anchors */
2876: DMPlexGetAnchors(dm, &aSec, NULL);
2877: if (aSec) {
2878: PetscInt aStart, aEnd, a;
2880: for (d = 0; d < maxDof; ++d) indices[d] = d;
2881: PetscSectionGetChart(aSec, &aStart, &aEnd);
2882: for (a = aStart; a < aEnd; a++) {
2883: PetscInt dof, fdof, f;
2885: PetscSectionGetDof(aSec, a, &dof);
2886: if (dof) {
2887: /* if there are point-to-point constraints, then all dofs are constrained */
2888: for (f = 0; f < numFields; f++) {
2889: PetscSectionGetFieldDof(section, a, f, &fdof);
2890: PetscSectionSetFieldConstraintIndices(section, a, f, indices);
2891: }
2892: }
2893: }
2894: }
2895: PetscFree(indices);
2896: return(0);
2897: }
2901: /* Set the constrained indices on each point */
2902: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
2903: {
2904: PetscInt *indices;
2905: PetscInt numFields, maxDof, pStart, pEnd, p, f, d;
2909: PetscSectionGetNumFields(section, &numFields);
2910: PetscSectionGetMaxDof(section, &maxDof);
2911: PetscSectionGetChart(section, &pStart, &pEnd);
2912: PetscMalloc1(maxDof, &indices);
2913: for (d = 0; d < maxDof; ++d) indices[d] = -1;
2914: for (p = pStart; p < pEnd; ++p) {
2915: PetscInt cdof, d;
2917: PetscSectionGetConstraintDof(section, p, &cdof);
2918: if (cdof) {
2919: if (numFields) {
2920: PetscInt numConst = 0, foff = 0;
2922: for (f = 0; f < numFields; ++f) {
2923: const PetscInt *find;
2924: PetscInt fcdof, fdof;
2926: PetscSectionGetFieldDof(section, p, f, &fdof);
2927: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
2928: /* Change constraint numbering from field component to local dof number */
2929: PetscSectionGetFieldConstraintIndices(section, p, f, &find);
2930: for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
2931: numConst += fcdof;
2932: foff += fdof;
2933: }
2934: if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
2935: } else {
2936: for (d = 0; d < cdof; ++d) indices[d] = d;
2937: }
2938: PetscSectionSetConstraintIndices(section, p, indices);
2939: }
2940: }
2941: PetscFree(indices);
2942: return(0);
2943: }
2947: /*@C
2948: DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided.
2950: Not Collective
2952: Input Parameters:
2953: + dm - The DMPlex object
2954: . dim - The spatial dimension of the problem
2955: . numFields - The number of fields in the problem
2956: . numComp - An array of size numFields that holds the number of components for each field
2957: . numDof - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
2958: . numBC - The number of boundary conditions
2959: . bcField - An array of size numBC giving the field number for each boundry condition
2960: . bcComps - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
2961: . bcPoints - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
2962: - perm - Optional permutation of the chart, or NULL
2964: Output Parameter:
2965: . section - The PetscSection object
2967: 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
2968: number of dof for field 0 on each edge.
2970: The chart permutation is the same one set using PetscSectionSetPermutation()
2972: Level: developer
2974: Fortran Notes:
2975: A Fortran 90 version is available as DMPlexCreateSectionF90()
2977: .keywords: mesh, elements
2978: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
2979: @*/
2980: 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)
2981: {
2982: PetscSection aSec;
2986: DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
2987: DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
2988: if (perm) {PetscSectionSetPermutation(*section, perm);}
2989: PetscSectionSetUp(*section);
2990: DMPlexGetAnchors(dm,&aSec,NULL);
2991: if (numBC || aSec) {
2992: DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
2993: DMPlexCreateSectionBCIndices(dm, *section);
2994: }
2995: PetscSectionViewFromOptions(*section,NULL,"-section_view");
2996: return(0);
2997: }
3001: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
3002: {
3003: PetscSection section, s;
3004: Mat m;
3008: DMClone(dm, cdm);
3009: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
3010: DMSetDefaultSection(*cdm, section);
3011: PetscSectionDestroy(§ion);
3012: PetscSectionCreate(PETSC_COMM_SELF, &s);
3013: MatCreate(PETSC_COMM_SELF, &m);
3014: DMSetDefaultConstraints(*cdm, s, m);
3015: PetscSectionDestroy(&s);
3016: MatDestroy(&m);
3017: return(0);
3018: }
3022: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
3023: {
3024: DM_Plex *mesh = (DM_Plex*) dm->data;
3028: if (section) *section = mesh->coneSection;
3029: return(0);
3030: }
3034: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
3035: {
3036: DM_Plex *mesh = (DM_Plex*) dm->data;
3040: if (section) *section = mesh->supportSection;
3041: return(0);
3042: }
3046: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
3047: {
3048: DM_Plex *mesh = (DM_Plex*) dm->data;
3052: if (cones) *cones = mesh->cones;
3053: return(0);
3054: }
3058: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
3059: {
3060: DM_Plex *mesh = (DM_Plex*) dm->data;
3064: if (coneOrientations) *coneOrientations = mesh->coneOrientations;
3065: return(0);
3066: }
3068: /******************************** FEM Support **********************************/
3072: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3073: {
3074: PetscScalar *array, *vArray;
3075: const PetscInt *cone, *coneO;
3076: PetscInt pStart, pEnd, p, numPoints, size = 0, offset = 0;
3077: PetscErrorCode ierr;
3080: PetscSectionGetChart(section, &pStart, &pEnd);
3081: DMPlexGetConeSize(dm, point, &numPoints);
3082: DMPlexGetCone(dm, point, &cone);
3083: DMPlexGetConeOrientation(dm, point, &coneO);
3084: if (!values || !*values) {
3085: if ((point >= pStart) && (point < pEnd)) {
3086: PetscInt dof;
3088: PetscSectionGetDof(section, point, &dof);
3089: size += dof;
3090: }
3091: for (p = 0; p < numPoints; ++p) {
3092: const PetscInt cp = cone[p];
3093: PetscInt dof;
3095: if ((cp < pStart) || (cp >= pEnd)) continue;
3096: PetscSectionGetDof(section, cp, &dof);
3097: size += dof;
3098: }
3099: if (!values) {
3100: if (csize) *csize = size;
3101: return(0);
3102: }
3103: DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
3104: } else {
3105: array = *values;
3106: }
3107: size = 0;
3108: VecGetArray(v, &vArray);
3109: if ((point >= pStart) && (point < pEnd)) {
3110: PetscInt dof, off, d;
3111: PetscScalar *varr;
3113: PetscSectionGetDof(section, point, &dof);
3114: PetscSectionGetOffset(section, point, &off);
3115: varr = &vArray[off];
3116: for (d = 0; d < dof; ++d, ++offset) {
3117: array[offset] = varr[d];
3118: }
3119: size += dof;
3120: }
3121: for (p = 0; p < numPoints; ++p) {
3122: const PetscInt cp = cone[p];
3123: PetscInt o = coneO[p];
3124: PetscInt dof, off, d;
3125: PetscScalar *varr;
3127: if ((cp < pStart) || (cp >= pEnd)) continue;
3128: PetscSectionGetDof(section, cp, &dof);
3129: PetscSectionGetOffset(section, cp, &off);
3130: varr = &vArray[off];
3131: if (o >= 0) {
3132: for (d = 0; d < dof; ++d, ++offset) {
3133: array[offset] = varr[d];
3134: }
3135: } else {
3136: for (d = dof-1; d >= 0; --d, ++offset) {
3137: array[offset] = varr[d];
3138: }
3139: }
3140: size += dof;
3141: }
3142: VecRestoreArray(v, &vArray);
3143: if (!*values) {
3144: if (csize) *csize = size;
3145: *values = array;
3146: } else {
3147: if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3148: *csize = size;
3149: }
3150: return(0);
3151: }
3155: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3156: {
3157: PetscInt offset = 0, p;
3161: *size = 0;
3162: for (p = 0; p < numPoints*2; p += 2) {
3163: const PetscInt point = points[p];
3164: const PetscInt o = points[p+1];
3165: PetscInt dof, off, d;
3166: const PetscScalar *varr;
3168: PetscSectionGetDof(section, point, &dof);
3169: PetscSectionGetOffset(section, point, &off);
3170: varr = &vArray[off];
3171: if (o >= 0) {
3172: for (d = 0; d < dof; ++d, ++offset) array[offset] = varr[d];
3173: } else {
3174: for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
3175: }
3176: }
3177: *size = offset;
3178: return(0);
3179: }
3183: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Fields_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], PetscInt numFields, const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3184: {
3185: PetscInt offset = 0, f;
3189: *size = 0;
3190: for (f = 0; f < numFields; ++f) {
3191: PetscInt fcomp, p;
3193: PetscSectionGetFieldComponents(section, f, &fcomp);
3194: for (p = 0; p < numPoints*2; p += 2) {
3195: const PetscInt point = points[p];
3196: const PetscInt o = points[p+1];
3197: PetscInt fdof, foff, d, c;
3198: const PetscScalar *varr;
3200: PetscSectionGetFieldDof(section, point, f, &fdof);
3201: PetscSectionGetFieldOffset(section, point, f, &foff);
3202: varr = &vArray[foff];
3203: if (o >= 0) {
3204: for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
3205: } else {
3206: for (d = fdof/fcomp-1; d >= 0; --d) {
3207: for (c = 0; c < fcomp; ++c, ++offset) {
3208: array[offset] = varr[d*fcomp+c];
3209: }
3210: }
3211: }
3212: }
3213: }
3214: *size = offset;
3215: return(0);
3216: }
3220: /*@C
3221: DMPlexVecGetClosure - Get an array of the values on the closure of 'point'
3223: Not collective
3225: Input Parameters:
3226: + dm - The DM
3227: . section - The section describing the layout in v, or NULL to use the default section
3228: . v - The local vector
3229: - point - The sieve point in the DM
3231: Output Parameters:
3232: + csize - The number of values in the closure, or NULL
3233: - values - The array of values, which is a borrowed array and should not be freed
3235: Fortran Notes:
3236: Since it returns an array, this routine is only available in Fortran 90, and you must
3237: include petsc.h90 in your code.
3239: The csize argument is not present in the Fortran 90 binding since it is internal to the array.
3241: Level: intermediate
3243: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3244: @*/
3245: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3246: {
3247: PetscSection clSection;
3248: IS clPoints;
3249: PetscScalar *array, *vArray;
3250: PetscInt *points = NULL;
3251: const PetscInt *clp;
3252: PetscInt depth, numFields, numPoints, size;
3253: PetscErrorCode ierr;
3257: if (!section) {DMGetDefaultSection(dm, §ion);}
3260: DMPlexGetDepth(dm, &depth);
3261: PetscSectionGetNumFields(section, &numFields);
3262: if (depth == 1 && numFields < 2) {
3263: DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
3264: return(0);
3265: }
3266: /* Get points */
3267: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3268: if (!clPoints) {
3269: PetscInt pStart, pEnd, p, q;
3271: PetscSectionGetChart(section, &pStart, &pEnd);
3272: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3273: /* Compress out points not in the section */
3274: for (p = 0, q = 0; p < numPoints*2; p += 2) {
3275: if ((points[p] >= pStart) && (points[p] < pEnd)) {
3276: points[q*2] = points[p];
3277: points[q*2+1] = points[p+1];
3278: ++q;
3279: }
3280: }
3281: numPoints = q;
3282: } else {
3283: PetscInt dof, off;
3285: PetscSectionGetDof(clSection, point, &dof);
3286: PetscSectionGetOffset(clSection, point, &off);
3287: ISGetIndices(clPoints, &clp);
3288: numPoints = dof/2;
3289: points = (PetscInt *) &clp[off];
3290: }
3291: /* Get array */
3292: if (!values || !*values) {
3293: PetscInt asize = 0, dof, p;
3295: for (p = 0; p < numPoints*2; p += 2) {
3296: PetscSectionGetDof(section, points[p], &dof);
3297: asize += dof;
3298: }
3299: if (!values) {
3300: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3301: else {ISRestoreIndices(clPoints, &clp);}
3302: if (csize) *csize = asize;
3303: return(0);
3304: }
3305: DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
3306: } else {
3307: array = *values;
3308: }
3309: VecGetArray(v, &vArray);
3310: /* Get values */
3311: if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);}
3312: else {DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);}
3313: /* Cleanup points */
3314: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3315: else {ISRestoreIndices(clPoints, &clp);}
3316: /* Cleanup array */
3317: VecRestoreArray(v, &vArray);
3318: if (!*values) {
3319: if (csize) *csize = size;
3320: *values = array;
3321: } else {
3322: if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %D < actual size %D", *csize, size);
3323: *csize = size;
3324: }
3325: return(0);
3326: }
3330: /*@C
3331: DMPlexVecRestoreClosure - Restore the array of the values on the closure of 'point'
3333: Not collective
3335: Input Parameters:
3336: + dm - The DM
3337: . section - The section describing the layout in v, or NULL to use the default section
3338: . v - The local vector
3339: . point - The sieve point in the DM
3340: . csize - The number of values in the closure, or NULL
3341: - values - The array of values, which is a borrowed array and should not be freed
3343: Fortran Notes:
3344: Since it returns an array, this routine is only available in Fortran 90, and you must
3345: include petsc.h90 in your code.
3347: The csize argument is not present in the Fortran 90 binding since it is internal to the array.
3349: Level: intermediate
3351: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3352: @*/
3353: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3354: {
3355: PetscInt size = 0;
3359: /* Should work without recalculating size */
3360: DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
3361: return(0);
3362: }
3364: PETSC_STATIC_INLINE void add (PetscScalar *x, PetscScalar y) {*x += y;}
3365: PETSC_STATIC_INLINE void insert(PetscScalar *x, PetscScalar y) {*x = y;}
3369: 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[])
3370: {
3371: PetscInt cdof; /* The number of constraints on this point */
3372: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3373: PetscScalar *a;
3374: PetscInt off, cind = 0, k;
3375: PetscErrorCode ierr;
3378: PetscSectionGetConstraintDof(section, point, &cdof);
3379: PetscSectionGetOffset(section, point, &off);
3380: a = &array[off];
3381: if (!cdof || setBC) {
3382: if (orientation >= 0) {
3383: for (k = 0; k < dof; ++k) {
3384: fuse(&a[k], values[k]);
3385: }
3386: } else {
3387: for (k = 0; k < dof; ++k) {
3388: fuse(&a[k], values[dof-k-1]);
3389: }
3390: }
3391: } else {
3392: PetscSectionGetConstraintIndices(section, point, &cdofs);
3393: if (orientation >= 0) {
3394: for (k = 0; k < dof; ++k) {
3395: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3396: fuse(&a[k], values[k]);
3397: }
3398: } else {
3399: for (k = 0; k < dof; ++k) {
3400: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3401: fuse(&a[k], values[dof-k-1]);
3402: }
3403: }
3404: }
3405: return(0);
3406: }
3410: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3411: {
3412: PetscInt cdof; /* The number of constraints on this point */
3413: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3414: PetscScalar *a;
3415: PetscInt off, cind = 0, k;
3416: PetscErrorCode ierr;
3419: PetscSectionGetConstraintDof(section, point, &cdof);
3420: PetscSectionGetOffset(section, point, &off);
3421: a = &array[off];
3422: if (cdof) {
3423: PetscSectionGetConstraintIndices(section, point, &cdofs);
3424: if (orientation >= 0) {
3425: for (k = 0; k < dof; ++k) {
3426: if ((cind < cdof) && (k == cdofs[cind])) {
3427: fuse(&a[k], values[k]);
3428: ++cind;
3429: }
3430: }
3431: } else {
3432: for (k = 0; k < dof; ++k) {
3433: if ((cind < cdof) && (k == cdofs[cind])) {
3434: fuse(&a[k], values[dof-k-1]);
3435: ++cind;
3436: }
3437: }
3438: }
3439: }
3440: return(0);
3441: }
3445: 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[])
3446: {
3447: PetscScalar *a;
3448: PetscInt fdof, foff, fcdof, foffset = *offset;
3449: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3450: PetscInt cind = 0, k, c;
3451: PetscErrorCode ierr;
3454: PetscSectionGetFieldDof(section, point, f, &fdof);
3455: PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3456: PetscSectionGetFieldOffset(section, point, f, &foff);
3457: a = &array[foff];
3458: if (!fcdof || setBC) {
3459: if (o >= 0) {
3460: for (k = 0; k < fdof; ++k) fuse(&a[k], values[foffset+k]);
3461: } else {
3462: for (k = fdof/fcomp-1; k >= 0; --k) {
3463: for (c = 0; c < fcomp; ++c) {
3464: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3465: }
3466: }
3467: }
3468: } else {
3469: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3470: if (o >= 0) {
3471: for (k = 0; k < fdof; ++k) {
3472: if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
3473: fuse(&a[k], values[foffset+k]);
3474: }
3475: } else {
3476: for (k = fdof/fcomp-1; k >= 0; --k) {
3477: for (c = 0; c < fcomp; ++c) {
3478: if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
3479: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3480: }
3481: }
3482: }
3483: }
3484: *offset += fdof;
3485: return(0);
3486: }
3490: 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[])
3491: {
3492: PetscScalar *a;
3493: PetscInt fdof, foff, fcdof, foffset = *offset;
3494: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3495: PetscInt cind = 0, k, c;
3496: PetscErrorCode ierr;
3499: PetscSectionGetFieldDof(section, point, f, &fdof);
3500: PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3501: PetscSectionGetFieldOffset(section, point, f, &foff);
3502: a = &array[foff];
3503: if (fcdof) {
3504: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3505: if (o >= 0) {
3506: for (k = 0; k < fdof; ++k) {
3507: if ((cind < fcdof) && (k == fcdofs[cind])) {
3508: fuse(&a[k], values[foffset+k]);
3509: ++cind;
3510: }
3511: }
3512: } else {
3513: for (k = fdof/fcomp-1; k >= 0; --k) {
3514: for (c = 0; c < fcomp; ++c) {
3515: if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {
3516: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3517: ++cind;
3518: }
3519: }
3520: }
3521: }
3522: }
3523: *offset += fdof;
3524: return(0);
3525: }
3529: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3530: {
3531: PetscScalar *array;
3532: const PetscInt *cone, *coneO;
3533: PetscInt pStart, pEnd, p, numPoints, off, dof;
3534: PetscErrorCode ierr;
3537: PetscSectionGetChart(section, &pStart, &pEnd);
3538: DMPlexGetConeSize(dm, point, &numPoints);
3539: DMPlexGetCone(dm, point, &cone);
3540: DMPlexGetConeOrientation(dm, point, &coneO);
3541: VecGetArray(v, &array);
3542: for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
3543: const PetscInt cp = !p ? point : cone[p-1];
3544: const PetscInt o = !p ? 0 : coneO[p-1];
3546: if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
3547: PetscSectionGetDof(section, cp, &dof);
3548: /* ADD_VALUES */
3549: {
3550: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3551: PetscScalar *a;
3552: PetscInt cdof, coff, cind = 0, k;
3554: PetscSectionGetConstraintDof(section, cp, &cdof);
3555: PetscSectionGetOffset(section, cp, &coff);
3556: a = &array[coff];
3557: if (!cdof) {
3558: if (o >= 0) {
3559: for (k = 0; k < dof; ++k) {
3560: a[k] += values[off+k];
3561: }
3562: } else {
3563: for (k = 0; k < dof; ++k) {
3564: a[k] += values[off+dof-k-1];
3565: }
3566: }
3567: } else {
3568: PetscSectionGetConstraintIndices(section, cp, &cdofs);
3569: if (o >= 0) {
3570: for (k = 0; k < dof; ++k) {
3571: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3572: a[k] += values[off+k];
3573: }
3574: } else {
3575: for (k = 0; k < dof; ++k) {
3576: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3577: a[k] += values[off+dof-k-1];
3578: }
3579: }
3580: }
3581: }
3582: }
3583: VecRestoreArray(v, &array);
3584: return(0);
3585: }
3589: /*@C
3590: DMPlexVecSetClosure - Set an array of the values on the closure of 'point'
3592: Not collective
3594: Input Parameters:
3595: + dm - The DM
3596: . section - The section describing the layout in v, or NULL to use the default section
3597: . v - The local vector
3598: . point - The sieve point in the DM
3599: . values - The array of values
3600: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
3602: Fortran Notes:
3603: This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
3605: Level: intermediate
3607: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
3608: @*/
3609: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3610: {
3611: PetscSection clSection;
3612: IS clPoints;
3613: PetscScalar *array;
3614: PetscInt *points = NULL;
3615: const PetscInt *clp;
3616: PetscInt depth, numFields, numPoints, p;
3617: PetscErrorCode ierr;
3621: if (!section) {DMGetDefaultSection(dm, §ion);}
3624: DMPlexGetDepth(dm, &depth);
3625: PetscSectionGetNumFields(section, &numFields);
3626: if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
3627: DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
3628: return(0);
3629: }
3630: /* Get points */
3631: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3632: if (!clPoints) {
3633: PetscInt pStart, pEnd, q;
3635: PetscSectionGetChart(section, &pStart, &pEnd);
3636: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3637: /* Compress out points not in the section */
3638: for (p = 0, q = 0; p < numPoints*2; p += 2) {
3639: if ((points[p] >= pStart) && (points[p] < pEnd)) {
3640: points[q*2] = points[p];
3641: points[q*2+1] = points[p+1];
3642: ++q;
3643: }
3644: }
3645: numPoints = q;
3646: } else {
3647: PetscInt dof, off;
3649: PetscSectionGetDof(clSection, point, &dof);
3650: PetscSectionGetOffset(clSection, point, &off);
3651: ISGetIndices(clPoints, &clp);
3652: numPoints = dof/2;
3653: points = (PetscInt *) &clp[off];
3654: }
3655: /* Get array */
3656: VecGetArray(v, &array);
3657: /* Get values */
3658: if (numFields > 0) {
3659: PetscInt offset = 0, fcomp, f;
3660: for (f = 0; f < numFields; ++f) {
3661: PetscSectionGetFieldComponents(section, f, &fcomp);
3662: switch (mode) {
3663: case INSERT_VALUES:
3664: for (p = 0; p < numPoints*2; p += 2) {
3665: const PetscInt point = points[p];
3666: const PetscInt o = points[p+1];
3667: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3668: } break;
3669: case INSERT_ALL_VALUES:
3670: for (p = 0; p < numPoints*2; p += 2) {
3671: const PetscInt point = points[p];
3672: const PetscInt o = points[p+1];
3673: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3674: } break;
3675: case INSERT_BC_VALUES:
3676: for (p = 0; p < numPoints*2; p += 2) {
3677: const PetscInt point = points[p];
3678: const PetscInt o = points[p+1];
3679: updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3680: } break;
3681: case ADD_VALUES:
3682: for (p = 0; p < numPoints*2; p += 2) {
3683: const PetscInt point = points[p];
3684: const PetscInt o = points[p+1];
3685: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3686: } break;
3687: case ADD_ALL_VALUES:
3688: for (p = 0; p < numPoints*2; p += 2) {
3689: const PetscInt point = points[p];
3690: const PetscInt o = points[p+1];
3691: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3692: } break;
3693: default:
3694: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3695: }
3696: }
3697: } else {
3698: PetscInt dof, off;
3700: switch (mode) {
3701: case INSERT_VALUES:
3702: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3703: PetscInt o = points[p+1];
3704: PetscSectionGetDof(section, points[p], &dof);
3705: updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
3706: } break;
3707: case INSERT_ALL_VALUES:
3708: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3709: PetscInt o = points[p+1];
3710: PetscSectionGetDof(section, points[p], &dof);
3711: updatePoint_private(section, points[p], dof, insert, PETSC_TRUE, o, &values[off], array);
3712: } break;
3713: case INSERT_BC_VALUES:
3714: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3715: PetscInt o = points[p+1];
3716: PetscSectionGetDof(section, points[p], &dof);
3717: updatePointBC_private(section, points[p], dof, insert, o, &values[off], array);
3718: } break;
3719: case ADD_VALUES:
3720: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3721: PetscInt o = points[p+1];
3722: PetscSectionGetDof(section, points[p], &dof);
3723: updatePoint_private(section, points[p], dof, add, PETSC_FALSE, o, &values[off], array);
3724: } break;
3725: case ADD_ALL_VALUES:
3726: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3727: PetscInt o = points[p+1];
3728: PetscSectionGetDof(section, points[p], &dof);
3729: updatePoint_private(section, points[p], dof, add, PETSC_TRUE, o, &values[off], array);
3730: } break;
3731: default:
3732: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3733: }
3734: }
3735: /* Cleanup points */
3736: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3737: else {ISRestoreIndices(clPoints, &clp);}
3738: /* Cleanup array */
3739: VecRestoreArray(v, &array);
3740: return(0);
3741: }
3745: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
3746: {
3747: PetscSection clSection;
3748: IS clPoints;
3749: PetscScalar *array;
3750: PetscInt *points = NULL;
3751: const PetscInt *clp;
3752: PetscInt numFields, numPoints, p;
3753: PetscInt offset = 0, fcomp, f;
3754: PetscErrorCode ierr;
3758: if (!section) {DMGetDefaultSection(dm, §ion);}
3761: PetscSectionGetNumFields(section, &numFields);
3762: /* Get points */
3763: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3764: if (!clPoints) {
3765: PetscInt pStart, pEnd, q;
3767: PetscSectionGetChart(section, &pStart, &pEnd);
3768: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3769: /* Compress out points not in the section */
3770: for (p = 0, q = 0; p < numPoints*2; p += 2) {
3771: if ((points[p] >= pStart) && (points[p] < pEnd)) {
3772: points[q*2] = points[p];
3773: points[q*2+1] = points[p+1];
3774: ++q;
3775: }
3776: }
3777: numPoints = q;
3778: } else {
3779: PetscInt dof, off;
3781: PetscSectionGetDof(clSection, point, &dof);
3782: PetscSectionGetOffset(clSection, point, &off);
3783: ISGetIndices(clPoints, &clp);
3784: numPoints = dof/2;
3785: points = (PetscInt *) &clp[off];
3786: }
3787: /* Get array */
3788: VecGetArray(v, &array);
3789: /* Get values */
3790: for (f = 0; f < numFields; ++f) {
3791: PetscSectionGetFieldComponents(section, f, &fcomp);
3792: if (!fieldActive[f]) {
3793: for (p = 0; p < numPoints*2; p += 2) {
3794: PetscInt fdof;
3795: PetscSectionGetFieldDof(section, points[p], f, &fdof);
3796: offset += fdof;
3797: }
3798: continue;
3799: }
3800: switch (mode) {
3801: case INSERT_VALUES:
3802: for (p = 0; p < numPoints*2; p += 2) {
3803: const PetscInt point = points[p];
3804: const PetscInt o = points[p+1];
3805: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3806: } break;
3807: case INSERT_ALL_VALUES:
3808: for (p = 0; p < numPoints*2; p += 2) {
3809: const PetscInt point = points[p];
3810: const PetscInt o = points[p+1];
3811: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3812: } break;
3813: case INSERT_BC_VALUES:
3814: for (p = 0; p < numPoints*2; p += 2) {
3815: const PetscInt point = points[p];
3816: const PetscInt o = points[p+1];
3817: updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3818: } break;
3819: case ADD_VALUES:
3820: for (p = 0; p < numPoints*2; p += 2) {
3821: const PetscInt point = points[p];
3822: const PetscInt o = points[p+1];
3823: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3824: } break;
3825: case ADD_ALL_VALUES:
3826: for (p = 0; p < numPoints*2; p += 2) {
3827: const PetscInt point = points[p];
3828: const PetscInt o = points[p+1];
3829: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3830: } break;
3831: default:
3832: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3833: }
3834: }
3835: /* Cleanup points */
3836: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3837: else {ISRestoreIndices(clPoints, &clp);}
3838: /* Cleanup array */
3839: VecRestoreArray(v, &array);
3840: return(0);
3841: }
3845: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
3846: {
3847: PetscMPIInt rank;
3848: PetscInt i, j;
3852: MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
3853: PetscViewerASCIIPrintf(viewer, "[%d]mat for sieve point %D\n", rank, point);
3854: for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
3855: for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
3856: numCIndices = numCIndices ? numCIndices : numRIndices;
3857: for (i = 0; i < numRIndices; i++) {
3858: PetscViewerASCIIPrintf(viewer, "[%d]", rank);
3859: for (j = 0; j < numCIndices; j++) {
3860: #if defined(PETSC_USE_COMPLEX)
3861: PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
3862: #else
3863: PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
3864: #endif
3865: }
3866: PetscViewerASCIIPrintf(viewer, "\n");
3867: }
3868: return(0);
3869: }
3873: /* . off - The global offset of this point */
3874: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
3875: {
3876: PetscInt dof; /* The number of unknowns on this point */
3877: PetscInt cdof; /* The number of constraints on this point */
3878: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3879: PetscInt cind = 0, k;
3880: PetscErrorCode ierr;
3883: PetscSectionGetDof(section, point, &dof);
3884: PetscSectionGetConstraintDof(section, point, &cdof);
3885: if (!cdof || setBC) {
3886: if (orientation >= 0) {
3887: for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
3888: } else {
3889: for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
3890: }
3891: } else {
3892: PetscSectionGetConstraintIndices(section, point, &cdofs);
3893: if (orientation >= 0) {
3894: for (k = 0; k < dof; ++k) {
3895: if ((cind < cdof) && (k == cdofs[cind])) {
3896: /* Insert check for returning constrained indices */
3897: indices[*loff+k] = -(off+k+1);
3898: ++cind;
3899: } else {
3900: indices[*loff+k] = off+k-cind;
3901: }
3902: }
3903: } else {
3904: for (k = 0; k < dof; ++k) {
3905: if ((cind < cdof) && (k == cdofs[cind])) {
3906: /* Insert check for returning constrained indices */
3907: indices[*loff+dof-k-1] = -(off+k+1);
3908: ++cind;
3909: } else {
3910: indices[*loff+dof-k-1] = off+k-cind;
3911: }
3912: }
3913: }
3914: }
3915: *loff += dof;
3916: return(0);
3917: }
3921: /* . off - The global offset of this point */
3922: PetscErrorCode indicesPointFields_private(PetscSection section, PetscInt point, PetscInt off, PetscInt foffs[], PetscBool setBC, PetscInt orientation, PetscInt indices[])
3923: {
3924: PetscInt numFields, foff, f;
3928: PetscSectionGetNumFields(section, &numFields);
3929: for (f = 0, foff = 0; f < numFields; ++f) {
3930: PetscInt fdof, fcomp, cfdof;
3931: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3932: PetscInt cind = 0, k, c;
3934: PetscSectionGetFieldComponents(section, f, &fcomp);
3935: PetscSectionGetFieldDof(section, point, f, &fdof);
3936: PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
3937: if (!cfdof || setBC) {
3938: if (orientation >= 0) {
3939: for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
3940: } else {
3941: for (k = fdof/fcomp-1; k >= 0; --k) {
3942: for (c = 0; c < fcomp; ++c) {
3943: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
3944: }
3945: }
3946: }
3947: } else {
3948: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3949: if (orientation >= 0) {
3950: for (k = 0; k < fdof; ++k) {
3951: if ((cind < cfdof) && (k == fcdofs[cind])) {
3952: indices[foffs[f]+k] = -(off+foff+k+1);
3953: ++cind;
3954: } else {
3955: indices[foffs[f]+k] = off+foff+k-cind;
3956: }
3957: }
3958: } else {
3959: for (k = fdof/fcomp-1; k >= 0; --k) {
3960: for (c = 0; c < fcomp; ++c) {
3961: if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
3962: indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
3963: ++cind;
3964: } else {
3965: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
3966: }
3967: }
3968: }
3969: }
3970: }
3971: foff += (setBC ? fdof : (fdof - cfdof));
3972: foffs[f] += fdof;
3973: }
3974: return(0);
3975: }
3979: 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[], PetscBool multiplyLeft)
3980: {
3981: Mat cMat;
3982: PetscSection aSec, cSec;
3983: IS aIS;
3984: PetscInt aStart = -1, aEnd = -1;
3985: const PetscInt *anchors;
3986: PetscInt numFields, f, p, q, newP = 0;
3987: PetscInt newNumPoints = 0, newNumIndices = 0;
3988: PetscInt *newPoints, *indices, *newIndices;
3989: PetscInt maxAnchor, maxDof;
3990: PetscInt newOffsets[32];
3991: PetscInt *pointMatOffsets[32];
3992: PetscInt *newPointOffsets[32];
3993: PetscScalar *pointMat[32];
3994: PetscScalar *newValues=NULL,*tmpValues;
3995: PetscBool anyConstrained = PETSC_FALSE;
3996: PetscErrorCode ierr;
4001: PetscSectionGetNumFields(section, &numFields);
4003: DMPlexGetAnchors(dm,&aSec,&aIS);
4004: /* if there are point-to-point constraints */
4005: if (aSec) {
4006: PetscMemzero(newOffsets, 32 * sizeof(PetscInt));
4007: ISGetIndices(aIS,&anchors);
4008: PetscSectionGetChart(aSec,&aStart,&aEnd);
4009: /* figure out how many points are going to be in the new element matrix
4010: * (we allow double counting, because it's all just going to be summed
4011: * into the global matrix anyway) */
4012: for (p = 0; p < 2*numPoints; p+=2) {
4013: PetscInt b = points[p];
4014: PetscInt bDof = 0, bSecDof;
4016: PetscSectionGetDof(section,b,&bSecDof);
4017: if (!bSecDof) {
4018: continue;
4019: }
4020: if (b >= aStart && b < aEnd) {
4021: PetscSectionGetDof(aSec,b,&bDof);
4022: }
4023: if (bDof) {
4024: /* this point is constrained */
4025: /* it is going to be replaced by its anchors */
4026: PetscInt bOff, q;
4028: anyConstrained = PETSC_TRUE;
4029: newNumPoints += bDof;
4030: PetscSectionGetOffset(aSec,b,&bOff);
4031: for (q = 0; q < bDof; q++) {
4032: PetscInt a = anchors[bOff + q];
4033: PetscInt aDof;
4035: PetscSectionGetDof(section,a,&aDof);
4036: newNumIndices += aDof;
4037: for (f = 0; f < numFields; ++f) {
4038: PetscInt fDof;
4040: PetscSectionGetFieldDof(section, a, f, &fDof);
4041: newOffsets[f+1] += fDof;
4042: }
4043: }
4044: }
4045: else {
4046: /* this point is not constrained */
4047: newNumPoints++;
4048: newNumIndices += bSecDof;
4049: for (f = 0; f < numFields; ++f) {
4050: PetscInt fDof;
4052: PetscSectionGetFieldDof(section, b, f, &fDof);
4053: newOffsets[f+1] += fDof;
4054: }
4055: }
4056: }
4057: }
4058: if (!anyConstrained) {
4059: if (outNumPoints) *outNumPoints = 0;
4060: if (outNumIndices) *outNumIndices = 0;
4061: if (outPoints) *outPoints = NULL;
4062: if (outValues) *outValues = NULL;
4063: if (aSec) {ISRestoreIndices(aIS,&anchors);}
4064: return(0);
4065: }
4067: if (outNumPoints) *outNumPoints = newNumPoints;
4068: if (outNumIndices) *outNumIndices = newNumIndices;
4070: for (f = 1; f < numFields; ++f) newOffsets[f+1] += newOffsets[f];
4072: if (!outPoints && !outValues) {
4073: if (offsets) {
4074: for (f = 0; f <= numFields; f++) {
4075: offsets[f] = newOffsets[f];
4076: }
4077: }
4078: if (aSec) {ISRestoreIndices(aIS,&anchors);}
4079: return(0);
4080: }
4082: if (numFields && newOffsets[numFields] != newNumIndices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid size for closure %D should be %D", newOffsets[numFields], newNumIndices);
4084: DMGetDefaultConstraints(dm, &cSec, &cMat);
4086: /* workspaces */
4087: if (numFields) {
4088: for (f = 0; f < numFields; f++) {
4089: DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4090: DMGetWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4091: }
4092: }
4093: else {
4094: DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4095: DMGetWorkArray(dm,numPoints,PETSC_INT,&newPointOffsets[0]);
4096: }
4098: /* get workspaces for the point-to-point matrices */
4099: if (numFields) {
4100: PetscInt totalOffset, totalMatOffset;
4102: for (p = 0; p < numPoints; p++) {
4103: PetscInt b = points[2*p];
4104: PetscInt bDof = 0, bSecDof;
4106: PetscSectionGetDof(section,b,&bSecDof);
4107: if (!bSecDof) {
4108: for (f = 0; f < numFields; f++) {
4109: newPointOffsets[f][p + 1] = 0;
4110: pointMatOffsets[f][p + 1] = 0;
4111: }
4112: continue;
4113: }
4114: if (b >= aStart && b < aEnd) {
4115: PetscSectionGetDof(aSec, b, &bDof);
4116: }
4117: if (bDof) {
4118: for (f = 0; f < numFields; f++) {
4119: PetscInt fDof, q, bOff, allFDof = 0;
4121: PetscSectionGetFieldDof(section, b, f, &fDof);
4122: PetscSectionGetOffset(aSec, b, &bOff);
4123: for (q = 0; q < bDof; q++) {
4124: PetscInt a = anchors[bOff + q];
4125: PetscInt aFDof;
4127: PetscSectionGetFieldDof(section, a, f, &aFDof);
4128: allFDof += aFDof;
4129: }
4130: newPointOffsets[f][p+1] = allFDof;
4131: pointMatOffsets[f][p+1] = fDof * allFDof;
4132: }
4133: }
4134: else {
4135: for (f = 0; f < numFields; f++) {
4136: PetscInt fDof;
4138: PetscSectionGetFieldDof(section, b, f, &fDof);
4139: newPointOffsets[f][p+1] = fDof;
4140: pointMatOffsets[f][p+1] = 0;
4141: }
4142: }
4143: }
4144: for (f = 0, totalOffset = 0, totalMatOffset = 0; f < numFields; f++) {
4145: newPointOffsets[f][0] = totalOffset;
4146: pointMatOffsets[f][0] = totalMatOffset;
4147: for (p = 0; p < numPoints; p++) {
4148: newPointOffsets[f][p+1] += newPointOffsets[f][p];
4149: pointMatOffsets[f][p+1] += pointMatOffsets[f][p];
4150: }
4151: totalOffset = newPointOffsets[f][numPoints];
4152: totalMatOffset = pointMatOffsets[f][numPoints];
4153: DMGetWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4154: }
4155: }
4156: else {
4157: for (p = 0; p < numPoints; p++) {
4158: PetscInt b = points[2*p];
4159: PetscInt bDof = 0, bSecDof;
4161: PetscSectionGetDof(section,b,&bSecDof);
4162: if (!bSecDof) {
4163: newPointOffsets[0][p + 1] = 0;
4164: pointMatOffsets[0][p + 1] = 0;
4165: continue;
4166: }
4167: if (b >= aStart && b < aEnd) {
4168: PetscSectionGetDof(aSec, b, &bDof);
4169: }
4170: if (bDof) {
4171: PetscInt bOff, q, allDof = 0;
4173: PetscSectionGetOffset(aSec, b, &bOff);
4174: for (q = 0; q < bDof; q++) {
4175: PetscInt a = anchors[bOff + q], aDof;
4177: PetscSectionGetDof(section, a, &aDof);
4178: allDof += aDof;
4179: }
4180: newPointOffsets[0][p+1] = allDof;
4181: pointMatOffsets[0][p+1] = bSecDof * allDof;
4182: }
4183: else {
4184: newPointOffsets[0][p+1] = bSecDof;
4185: pointMatOffsets[0][p+1] = 0;
4186: }
4187: }
4188: newPointOffsets[0][0] = 0;
4189: pointMatOffsets[0][0] = 0;
4190: for (p = 0; p < numPoints; p++) {
4191: newPointOffsets[0][p+1] += newPointOffsets[0][p];
4192: pointMatOffsets[0][p+1] += pointMatOffsets[0][p];
4193: }
4194: DMGetWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4195: }
4197: /* output arrays */
4198: DMGetWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4200: /* get the point-to-point matrices; construct newPoints */
4201: PetscSectionGetMaxDof(aSec, &maxAnchor);
4202: PetscSectionGetMaxDof(section, &maxDof);
4203: DMGetWorkArray(dm,maxDof,PETSC_INT,&indices);
4204: DMGetWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4205: if (numFields) {
4206: for (p = 0, newP = 0; p < numPoints; p++) {
4207: PetscInt b = points[2*p];
4208: PetscInt o = points[2*p+1];
4209: PetscInt bDof = 0, bSecDof;
4211: PetscSectionGetDof(section, b, &bSecDof);
4212: if (!bSecDof) {
4213: continue;
4214: }
4215: if (b >= aStart && b < aEnd) {
4216: PetscSectionGetDof(aSec, b, &bDof);
4217: }
4218: if (bDof) {
4219: PetscInt fStart[32], fEnd[32], fAnchorStart[32], fAnchorEnd[32], bOff, q;
4221: fStart[0] = 0;
4222: fEnd[0] = 0;
4223: for (f = 0; f < numFields; f++) {
4224: PetscInt fDof;
4226: PetscSectionGetFieldDof(cSec, b, f, &fDof);
4227: fStart[f+1] = fStart[f] + fDof;
4228: fEnd[f+1] = fStart[f+1];
4229: }
4230: PetscSectionGetOffset(cSec, b, &bOff);
4231: indicesPointFields_private(cSec, b, bOff, fEnd, PETSC_TRUE, o, indices);
4233: fAnchorStart[0] = 0;
4234: fAnchorEnd[0] = 0;
4235: for (f = 0; f < numFields; f++) {
4236: PetscInt fDof = newPointOffsets[f][p + 1] - newPointOffsets[f][p];
4238: fAnchorStart[f+1] = fAnchorStart[f] + fDof;
4239: fAnchorEnd[f+1] = fAnchorStart[f + 1];
4240: }
4241: PetscSectionGetOffset(aSec, b, &bOff);
4242: for (q = 0; q < bDof; q++) {
4243: PetscInt a = anchors[bOff + q], aOff;
4245: /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4246: newPoints[2*(newP + q)] = a;
4247: newPoints[2*(newP + q) + 1] = 0;
4248: PetscSectionGetOffset(section, a, &aOff);
4249: indicesPointFields_private(section, a, aOff, fAnchorEnd, PETSC_TRUE, 0, newIndices);
4250: }
4251: newP += bDof;
4253: if (outValues) {
4254: /* get the point-to-point submatrix */
4255: for (f = 0; f < numFields; f++) {
4256: MatGetValues(cMat,fEnd[f]-fStart[f],indices + fStart[f],fAnchorEnd[f] - fAnchorStart[f],newIndices + fAnchorStart[f],pointMat[f] + pointMatOffsets[f][p]);
4257: }
4258: }
4259: }
4260: else {
4261: newPoints[2 * newP] = b;
4262: newPoints[2 * newP + 1] = o;
4263: newP++;
4264: }
4265: }
4266: } else {
4267: for (p = 0; p < numPoints; p++) {
4268: PetscInt b = points[2*p];
4269: PetscInt o = points[2*p+1];
4270: PetscInt bDof = 0, bSecDof;
4272: PetscSectionGetDof(section, b, &bSecDof);
4273: if (!bSecDof) {
4274: continue;
4275: }
4276: if (b >= aStart && b < aEnd) {
4277: PetscSectionGetDof(aSec, b, &bDof);
4278: }
4279: if (bDof) {
4280: PetscInt bEnd = 0, bAnchorEnd = 0, bOff;
4282: PetscSectionGetOffset(cSec, b, &bOff);
4283: indicesPoint_private(cSec, b, bOff, &bEnd, PETSC_TRUE, o, indices);
4285: PetscSectionGetOffset (aSec, b, &bOff);
4286: for (q = 0; q < bDof; q++) {
4287: PetscInt a = anchors[bOff + q], aOff;
4289: /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4291: newPoints[2*(newP + q)] = a;
4292: newPoints[2*(newP + q) + 1] = 0;
4293: PetscSectionGetOffset(section, a, &aOff);
4294: indicesPoint_private(section, a, aOff, &bAnchorEnd, PETSC_TRUE, 0, newIndices);
4295: }
4296: newP += bDof;
4298: /* get the point-to-point submatrix */
4299: if (outValues) {
4300: MatGetValues(cMat,bEnd,indices,bAnchorEnd,newIndices,pointMat[0] + pointMatOffsets[0][p]);
4301: }
4302: }
4303: else {
4304: newPoints[2 * newP] = b;
4305: newPoints[2 * newP + 1] = o;
4306: newP++;
4307: }
4308: }
4309: }
4311: if (outValues) {
4312: DMGetWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4313: PetscMemzero(tmpValues,newNumIndices*numIndices*sizeof(*tmpValues));
4314: /* multiply constraints on the right */
4315: if (numFields) {
4316: for (f = 0; f < numFields; f++) {
4317: PetscInt oldOff = offsets[f];
4319: for (p = 0; p < numPoints; p++) {
4320: PetscInt cStart = newPointOffsets[f][p];
4321: PetscInt b = points[2 * p];
4322: PetscInt c, r, k;
4323: PetscInt dof;
4325: PetscSectionGetFieldDof(section,b,f,&dof);
4326: if (!dof) {
4327: continue;
4328: }
4329: if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4330: PetscInt nCols = newPointOffsets[f][p+1]-cStart;
4331: const PetscScalar *mat = pointMat[f] + pointMatOffsets[f][p];
4333: for (r = 0; r < numIndices; r++) {
4334: for (c = 0; c < nCols; c++) {
4335: for (k = 0; k < dof; k++) {
4336: tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4337: }
4338: }
4339: }
4340: }
4341: else {
4342: /* copy this column as is */
4343: for (r = 0; r < numIndices; r++) {
4344: for (c = 0; c < dof; c++) {
4345: tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4346: }
4347: }
4348: }
4349: oldOff += dof;
4350: }
4351: }
4352: }
4353: else {
4354: PetscInt oldOff = 0;
4355: for (p = 0; p < numPoints; p++) {
4356: PetscInt cStart = newPointOffsets[0][p];
4357: PetscInt b = points[2 * p];
4358: PetscInt c, r, k;
4359: PetscInt dof;
4361: PetscSectionGetDof(section,b,&dof);
4362: if (!dof) {
4363: continue;
4364: }
4365: if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4366: PetscInt nCols = newPointOffsets[0][p+1]-cStart;
4367: const PetscScalar *mat = pointMat[0] + pointMatOffsets[0][p];
4369: for (r = 0; r < numIndices; r++) {
4370: for (c = 0; c < nCols; c++) {
4371: for (k = 0; k < dof; k++) {
4372: tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4373: }
4374: }
4375: }
4376: }
4377: else {
4378: /* copy this column as is */
4379: for (r = 0; r < numIndices; r++) {
4380: for (c = 0; c < dof; c++) {
4381: tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4382: }
4383: }
4384: }
4385: oldOff += dof;
4386: }
4387: }
4389: if (multiplyLeft) {
4390: DMGetWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4391: PetscMemzero(newValues,newNumIndices*newNumIndices*sizeof(*newValues));
4392: /* multiply constraints transpose on the left */
4393: if (numFields) {
4394: for (f = 0; f < numFields; f++) {
4395: PetscInt oldOff = offsets[f];
4397: for (p = 0; p < numPoints; p++) {
4398: PetscInt rStart = 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 nRows = newPointOffsets[f][p+1]-rStart;
4406: const PetscScalar *PETSC_RESTRICT mat = pointMat[f] + pointMatOffsets[f][p];
4408: for (r = 0; r < nRows; r++) {
4409: for (c = 0; c < newNumIndices; c++) {
4410: for (k = 0; k < dof; k++) {
4411: newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4412: }
4413: }
4414: }
4415: }
4416: else {
4417: /* copy this row as is */
4418: for (r = 0; r < dof; r++) {
4419: for (c = 0; c < newNumIndices; c++) {
4420: newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4421: }
4422: }
4423: }
4424: oldOff += dof;
4425: }
4426: }
4427: }
4428: else {
4429: PetscInt oldOff = 0;
4431: for (p = 0; p < numPoints; p++) {
4432: PetscInt rStart = newPointOffsets[0][p];
4433: PetscInt b = points[2 * p];
4434: PetscInt c, r, k;
4435: PetscInt dof;
4437: PetscSectionGetDof(section,b,&dof);
4438: if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4439: PetscInt nRows = newPointOffsets[0][p+1]-rStart;
4440: const PetscScalar *PETSC_RESTRICT mat = pointMat[0] + pointMatOffsets[0][p];
4442: for (r = 0; r < nRows; r++) {
4443: for (c = 0; c < newNumIndices; c++) {
4444: for (k = 0; k < dof; k++) {
4445: newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4446: }
4447: }
4448: }
4449: }
4450: else {
4451: /* copy this row as is */
4452: for (r = 0; r < dof; r++) {
4453: for (c = 0; c < newNumIndices; c++) {
4454: newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4455: }
4456: }
4457: }
4458: oldOff += dof;
4459: }
4460: }
4462: DMRestoreWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4463: }
4464: else {
4465: newValues = tmpValues;
4466: }
4467: }
4469: /* clean up */
4470: DMRestoreWorkArray(dm,maxDof,PETSC_INT,&indices);
4471: DMRestoreWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4473: if (numFields) {
4474: for (f = 0; f < numFields; f++) {
4475: DMRestoreWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4476: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4477: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4478: }
4479: }
4480: else {
4481: DMRestoreWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4482: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4483: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[0]);
4484: }
4485: ISRestoreIndices(aIS,&anchors);
4487: /* output */
4488: if (outPoints) {
4489: *outPoints = newPoints;
4490: }
4491: else {
4492: DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4493: }
4494: if (outValues) {
4495: *outValues = newValues;
4496: }
4497: else {
4498: DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4499: }
4500: for (f = 0; f <= numFields; f++) {
4501: offsets[f] = newOffsets[f];
4502: }
4503: return(0);
4504: }
4508: PetscErrorCode DMPlexGetClosureIndices(DM dm, PetscSection section, PetscSection globalSection, PetscInt point, PetscInt *numIndices, PetscInt **indices, PetscInt *outOffsets)
4509: {
4510: PetscSection clSection;
4511: IS clPoints;
4512: const PetscInt *clp;
4513: PetscInt *points = NULL, *pointsNew;
4514: PetscInt numPoints, numPointsNew;
4515: PetscInt offsets[32];
4516: PetscInt Nf, Nind, NindNew, off, globalOff, f, p;
4517: PetscErrorCode ierr;
4525: PetscSectionGetNumFields(section, &Nf);
4526: if (Nf > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", Nf);
4527: PetscMemzero(offsets, 32 * sizeof(PetscInt));
4528: /* Get points in closure */
4529: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4530: if (!clPoints) {
4531: PetscInt pStart, pEnd, q;
4533: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4534: /* Compress out points not in the section */
4535: PetscSectionGetChart(section, &pStart, &pEnd);
4536: for (p = 0, q = 0; p < numPoints*2; p += 2) {
4537: if ((points[p] >= pStart) && (points[p] < pEnd)) {
4538: points[q*2] = points[p];
4539: points[q*2+1] = points[p+1];
4540: ++q;
4541: }
4542: }
4543: numPoints = q;
4544: } else {
4545: PetscInt dof, off;
4547: PetscSectionGetDof(clSection, point, &dof);
4548: numPoints = dof/2;
4549: PetscSectionGetOffset(clSection, point, &off);
4550: ISGetIndices(clPoints, &clp);
4551: points = (PetscInt *) &clp[off];
4552: }
4553: /* Get number of indices and indices per field */
4554: for (p = 0, Nind = 0; p < numPoints*2; p += 2) {
4555: PetscInt dof, fdof;
4557: PetscSectionGetDof(section, points[p], &dof);
4558: for (f = 0; f < Nf; ++f) {
4559: PetscSectionGetFieldDof(section, points[p], f, &fdof);
4560: offsets[f+1] += fdof;
4561: }
4562: Nind += dof;
4563: }
4564: for (f = 1; f < Nf; ++f) offsets[f+1] += offsets[f];
4565: if (Nf && offsets[Nf] != Nind) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[Nf], Nind);
4566: /* Correct for hanging node constraints */
4567: {
4568: DMPlexAnchorsModifyMat(dm, section, numPoints, Nind, points, NULL, &numPointsNew, &NindNew, &pointsNew, NULL, offsets, PETSC_TRUE);
4569: if (numPointsNew) {
4570: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4571: else {ISRestoreIndices(clPoints, &clp);}
4572: numPoints = numPointsNew;
4573: Nind = NindNew;
4574: points = pointsNew;
4575: }
4576: }
4577: /* Calculate indices */
4578: DMGetWorkArray(dm, Nind, PETSC_INT, indices);
4579: if (Nf) {
4580: if (outOffsets) {
4581: PetscInt f;
4583: for (f = 0; f <= Nf; f++) {
4584: outOffsets[f] = offsets[f];
4585: }
4586: }
4587: for (p = 0; p < numPoints*2; p += 2) {
4588: PetscInt o = points[p+1];
4589: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4590: indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, *indices);
4591: }
4592: } else {
4593: for (p = 0, off = 0; p < numPoints*2; p += 2) {
4594: PetscInt o = points[p+1];
4595: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4596: indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, *indices);
4597: }
4598: }
4599: /* Cleanup points */
4600: if (numPointsNew) {
4601: DMRestoreWorkArray(dm, 2*numPointsNew, PETSC_INT, &pointsNew);
4602: } else {
4603: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4604: else {ISRestoreIndices(clPoints, &clp);}
4605: }
4606: if (numIndices) *numIndices = Nind;
4607: return(0);
4608: }
4612: PetscErrorCode DMPlexRestoreClosureIndices(DM dm, PetscSection section, PetscSection globalSection, PetscInt point, PetscInt *numIndices, PetscInt **indices,PetscInt *outOffsets)
4613: {
4619: DMRestoreWorkArray(dm, 0, PETSC_INT, indices);
4620: return(0);
4621: }
4625: /*@C
4626: DMPlexMatSetClosure - Set an array of the values on the closure of 'point'
4628: Not collective
4630: Input Parameters:
4631: + dm - The DM
4632: . section - The section describing the layout in v, or NULL to use the default section
4633: . globalSection - The section describing the layout in v, or NULL to use the default global section
4634: . A - The matrix
4635: . point - The sieve point in the DM
4636: . values - The array of values
4637: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
4639: Fortran Notes:
4640: This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
4642: Level: intermediate
4644: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
4645: @*/
4646: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4647: {
4648: DM_Plex *mesh = (DM_Plex*) dm->data;
4649: PetscSection clSection;
4650: IS clPoints;
4651: PetscInt *points = NULL, *newPoints;
4652: const PetscInt *clp;
4653: PetscInt *indices;
4654: PetscInt offsets[32];
4655: PetscInt numFields, numPoints, newNumPoints, numIndices, newNumIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
4656: PetscScalar *newValues;
4657: PetscErrorCode ierr;
4661: if (!section) {DMGetDefaultSection(dm, §ion);}
4663: if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
4666: PetscSectionGetNumFields(section, &numFields);
4667: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4668: PetscMemzero(offsets, 32 * sizeof(PetscInt));
4669: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4670: if (!clPoints) {
4671: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4672: /* Compress out points not in the section */
4673: PetscSectionGetChart(section, &pStart, &pEnd);
4674: for (p = 0, q = 0; p < numPoints*2; p += 2) {
4675: if ((points[p] >= pStart) && (points[p] < pEnd)) {
4676: points[q*2] = points[p];
4677: points[q*2+1] = points[p+1];
4678: ++q;
4679: }
4680: }
4681: numPoints = q;
4682: } else {
4683: PetscInt dof, off;
4685: PetscSectionGetDof(clSection, point, &dof);
4686: numPoints = dof/2;
4687: PetscSectionGetOffset(clSection, point, &off);
4688: ISGetIndices(clPoints, &clp);
4689: points = (PetscInt *) &clp[off];
4690: }
4691: for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
4692: PetscInt fdof;
4694: PetscSectionGetDof(section, points[p], &dof);
4695: for (f = 0; f < numFields; ++f) {
4696: PetscSectionGetFieldDof(section, points[p], f, &fdof);
4697: offsets[f+1] += fdof;
4698: }
4699: numIndices += dof;
4700: }
4701: for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];
4703: if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
4704: DMPlexAnchorsModifyMat(dm,section,numPoints,numIndices,points,values,&newNumPoints,&newNumIndices,&newPoints,&newValues,offsets,PETSC_TRUE);
4705: if (newNumPoints) {
4706: if (!clPoints) {
4707: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4708: } else {
4709: ISRestoreIndices(clPoints, &clp);
4710: }
4711: numPoints = newNumPoints;
4712: numIndices = newNumIndices;
4713: points = newPoints;
4714: values = newValues;
4715: }
4716: DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
4717: if (numFields) {
4718: for (p = 0; p < numPoints*2; p += 2) {
4719: PetscInt o = points[p+1];
4720: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4721: indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
4722: }
4723: } else {
4724: for (p = 0, off = 0; p < numPoints*2; p += 2) {
4725: PetscInt o = points[p+1];
4726: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4727: indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
4728: }
4729: }
4730: if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
4731: MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
4732: if (mesh->printFEM > 1) {
4733: PetscInt i;
4734: PetscPrintf(PETSC_COMM_SELF, " Indices:");
4735: for (i = 0; i < numIndices; ++i) {PetscPrintf(PETSC_COMM_SELF, " %d", indices[i]);}
4736: PetscPrintf(PETSC_COMM_SELF, "\n");
4737: }
4738: if (ierr) {
4739: PetscMPIInt rank;
4740: PetscErrorCode ierr2;
4742: ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4743: ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4744: ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
4745: ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
4746:
4747: }
4748: if (newNumPoints) {
4749: DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4750: DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4751: }
4752: else {
4753: if (!clPoints) {
4754: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4755: } else {
4756: ISRestoreIndices(clPoints, &clp);
4757: }
4758: }
4759: DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
4760: return(0);
4761: }
4765: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4766: {
4767: DM_Plex *mesh = (DM_Plex*) dmf->data;
4768: PetscInt *fpoints = NULL, *ftotpoints = NULL;
4769: PetscInt *cpoints = NULL;
4770: PetscInt *findices, *cindices;
4771: PetscInt foffsets[32], coffsets[32];
4772: CellRefiner cellRefiner;
4773: PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4774: PetscErrorCode ierr;
4779: if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4781: if (!csection) {DMGetDefaultSection(dmc, &csection);}
4783: if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4785: if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4788: PetscSectionGetNumFields(fsection, &numFields);
4789: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4790: PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4791: PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4792: /* Column indices */
4793: DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4794: maxFPoints = numCPoints;
4795: /* Compress out points not in the section */
4796: /* TODO: Squeeze out points with 0 dof as well */
4797: PetscSectionGetChart(csection, &pStart, &pEnd);
4798: for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4799: if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4800: cpoints[q*2] = cpoints[p];
4801: cpoints[q*2+1] = cpoints[p+1];
4802: ++q;
4803: }
4804: }
4805: numCPoints = q;
4806: for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4807: PetscInt fdof;
4809: PetscSectionGetDof(csection, cpoints[p], &dof);
4810: if (!dof) continue;
4811: for (f = 0; f < numFields; ++f) {
4812: PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4813: coffsets[f+1] += fdof;
4814: }
4815: numCIndices += dof;
4816: }
4817: for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4818: /* Row indices */
4819: DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4820: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4821: DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4822: for (r = 0, q = 0; r < numSubcells; ++r) {
4823: /* TODO Map from coarse to fine cells */
4824: DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4825: /* Compress out points not in the section */
4826: PetscSectionGetChart(fsection, &pStart, &pEnd);
4827: for (p = 0; p < numFPoints*2; p += 2) {
4828: if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4829: PetscSectionGetDof(fsection, fpoints[p], &dof);
4830: if (!dof) continue;
4831: for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4832: if (s < q) continue;
4833: ftotpoints[q*2] = fpoints[p];
4834: ftotpoints[q*2+1] = fpoints[p+1];
4835: ++q;
4836: }
4837: }
4838: DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4839: }
4840: numFPoints = q;
4841: for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4842: PetscInt fdof;
4844: PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4845: if (!dof) continue;
4846: for (f = 0; f < numFields; ++f) {
4847: PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4848: foffsets[f+1] += fdof;
4849: }
4850: numFIndices += dof;
4851: }
4852: for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];
4854: if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4855: if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4856: DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4857: DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4858: if (numFields) {
4859: for (p = 0; p < numFPoints*2; p += 2) {
4860: PetscInt o = ftotpoints[p+1];
4861: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4862: indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4863: }
4864: for (p = 0; p < numCPoints*2; p += 2) {
4865: PetscInt o = cpoints[p+1];
4866: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4867: indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4868: }
4869: } else {
4870: for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4871: PetscInt o = ftotpoints[p+1];
4872: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4873: indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4874: }
4875: for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4876: PetscInt o = cpoints[p+1];
4877: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4878: indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4879: }
4880: }
4881: if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
4882: MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
4883: if (ierr) {
4884: PetscMPIInt rank;
4885: PetscErrorCode ierr2;
4887: ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4888: ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4889: ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
4890: ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
4891: ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
4892:
4893: }
4894: DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4895: DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4896: DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4897: DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4898: return(0);
4899: }
4903: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
4904: {
4905: PetscInt *fpoints = NULL, *ftotpoints = NULL;
4906: PetscInt *cpoints = NULL;
4907: PetscInt foffsets[32], coffsets[32];
4908: CellRefiner cellRefiner;
4909: PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4915: if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4917: if (!csection) {DMGetDefaultSection(dmc, &csection);}
4919: if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4921: if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4923: PetscSectionGetNumFields(fsection, &numFields);
4924: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4925: PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4926: PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4927: /* Column indices */
4928: DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4929: maxFPoints = numCPoints;
4930: /* Compress out points not in the section */
4931: /* TODO: Squeeze out points with 0 dof as well */
4932: PetscSectionGetChart(csection, &pStart, &pEnd);
4933: for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4934: if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4935: cpoints[q*2] = cpoints[p];
4936: cpoints[q*2+1] = cpoints[p+1];
4937: ++q;
4938: }
4939: }
4940: numCPoints = q;
4941: for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4942: PetscInt fdof;
4944: PetscSectionGetDof(csection, cpoints[p], &dof);
4945: if (!dof) continue;
4946: for (f = 0; f < numFields; ++f) {
4947: PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4948: coffsets[f+1] += fdof;
4949: }
4950: numCIndices += dof;
4951: }
4952: for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4953: /* Row indices */
4954: DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4955: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4956: DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4957: for (r = 0, q = 0; r < numSubcells; ++r) {
4958: /* TODO Map from coarse to fine cells */
4959: DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4960: /* Compress out points not in the section */
4961: PetscSectionGetChart(fsection, &pStart, &pEnd);
4962: for (p = 0; p < numFPoints*2; p += 2) {
4963: if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4964: PetscSectionGetDof(fsection, fpoints[p], &dof);
4965: if (!dof) continue;
4966: for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4967: if (s < q) continue;
4968: ftotpoints[q*2] = fpoints[p];
4969: ftotpoints[q*2+1] = fpoints[p+1];
4970: ++q;
4971: }
4972: }
4973: DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4974: }
4975: numFPoints = q;
4976: for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4977: PetscInt fdof;
4979: PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4980: if (!dof) continue;
4981: for (f = 0; f < numFields; ++f) {
4982: PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4983: foffsets[f+1] += fdof;
4984: }
4985: numFIndices += dof;
4986: }
4987: for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];
4989: if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4990: if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4991: if (numFields) {
4992: for (p = 0; p < numFPoints*2; p += 2) {
4993: PetscInt o = ftotpoints[p+1];
4994: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4995: indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4996: }
4997: for (p = 0; p < numCPoints*2; p += 2) {
4998: PetscInt o = cpoints[p+1];
4999: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5000: indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5001: }
5002: } else {
5003: for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5004: PetscInt o = ftotpoints[p+1];
5005: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5006: indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5007: }
5008: for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5009: PetscInt o = cpoints[p+1];
5010: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5011: indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5012: }
5013: }
5014: DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5015: DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5016: return(0);
5017: }
5021: /*@
5022: DMPlexGetHybridBounds - Get the first mesh point of each dimension which is a hybrid
5024: Input Parameter:
5025: . dm - The DMPlex object
5027: Output Parameters:
5028: + cMax - The first hybrid cell
5029: . fMax - The first hybrid face
5030: . eMax - The first hybrid edge
5031: - vMax - The first hybrid vertex
5033: Level: developer
5035: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
5036: @*/
5037: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
5038: {
5039: DM_Plex *mesh = (DM_Plex*) dm->data;
5040: PetscInt dim;
5045: DMGetDimension(dm, &dim);
5046: if (cMax) *cMax = mesh->hybridPointMax[dim];
5047: if (fMax) *fMax = mesh->hybridPointMax[dim-1];
5048: if (eMax) *eMax = mesh->hybridPointMax[1];
5049: if (vMax) *vMax = mesh->hybridPointMax[0];
5050: return(0);
5051: }
5055: /*@
5056: DMPlexSetHybridBounds - Set the first mesh point of each dimension which is a hybrid
5058: Input Parameters:
5059: . dm - The DMPlex object
5060: . cMax - The first hybrid cell
5061: . fMax - The first hybrid face
5062: . eMax - The first hybrid edge
5063: - vMax - The first hybrid vertex
5065: Level: developer
5067: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5068: @*/
5069: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5070: {
5071: DM_Plex *mesh = (DM_Plex*) dm->data;
5072: PetscInt dim;
5077: DMGetDimension(dm, &dim);
5078: if (cMax >= 0) mesh->hybridPointMax[dim] = cMax;
5079: if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5080: if (eMax >= 0) mesh->hybridPointMax[1] = eMax;
5081: if (vMax >= 0) mesh->hybridPointMax[0] = vMax;
5082: return(0);
5083: }
5087: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5088: {
5089: DM_Plex *mesh = (DM_Plex*) dm->data;
5094: *cellHeight = mesh->vtkCellHeight;
5095: return(0);
5096: }
5100: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5101: {
5102: DM_Plex *mesh = (DM_Plex*) dm->data;
5106: mesh->vtkCellHeight = cellHeight;
5107: return(0);
5108: }
5112: /* We can easily have a form that takes an IS instead */
5113: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5114: {
5115: PetscSection section, globalSection;
5116: PetscInt *numbers, p;
5120: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
5121: PetscSectionSetChart(section, pStart, pEnd);
5122: for (p = pStart; p < pEnd; ++p) {
5123: PetscSectionSetDof(section, p, 1);
5124: }
5125: PetscSectionSetUp(section);
5126: PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
5127: PetscMalloc1(pEnd - pStart, &numbers);
5128: for (p = pStart; p < pEnd; ++p) {
5129: PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5130: if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5131: else numbers[p-pStart] += shift;
5132: }
5133: ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5134: if (globalSize) {
5135: PetscLayout layout;
5136: PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5137: PetscLayoutGetSize(layout, globalSize);
5138: PetscLayoutDestroy(&layout);
5139: }
5140: PetscSectionDestroy(§ion);
5141: PetscSectionDestroy(&globalSection);
5142: return(0);
5143: }
5147: PetscErrorCode DMPlexCreateCellNumbering_Internal(DM dm, PetscBool includeHybrid, IS *globalCellNumbers)
5148: {
5149: PetscInt cellHeight, cStart, cEnd, cMax;
5153: DMPlexGetVTKCellHeight(dm, &cellHeight);
5154: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5155: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5156: if (cMax >= 0 && !includeHybrid) cEnd = PetscMin(cEnd, cMax);
5157: DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, globalCellNumbers);
5158: return(0);
5159: }
5163: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5164: {
5165: DM_Plex *mesh = (DM_Plex*) dm->data;
5170: if (!mesh->globalCellNumbers) {DMPlexCreateCellNumbering_Internal(dm, PETSC_FALSE, &mesh->globalCellNumbers);}
5171: *globalCellNumbers = mesh->globalCellNumbers;
5172: return(0);
5173: }
5177: PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM dm, PetscBool includeHybrid, IS *globalVertexNumbers)
5178: {
5179: PetscInt vStart, vEnd, vMax;
5184: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5185: DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5186: if (vMax >= 0 && !includeHybrid) vEnd = PetscMin(vEnd, vMax);
5187: DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, globalVertexNumbers);
5188: return(0);
5189: }
5193: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5194: {
5195: DM_Plex *mesh = (DM_Plex*) dm->data;
5200: if (!mesh->globalVertexNumbers) {DMPlexCreateVertexNumbering_Internal(dm, PETSC_FALSE, &mesh->globalVertexNumbers);}
5201: *globalVertexNumbers = mesh->globalVertexNumbers;
5202: return(0);
5203: }
5207: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5208: {
5209: IS nums[4];
5210: PetscInt depths[4];
5211: PetscInt depth, d, shift = 0;
5216: DMPlexGetDepth(dm, &depth);
5217: /* For unstratified meshes use dim instead of depth */
5218: if (depth < 0) {DMGetDimension(dm, &depth);}
5219: depths[0] = depth; depths[1] = 0;
5220: for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5221: for (d = 0; d <= depth; ++d) {
5222: PetscInt pStart, pEnd, gsize;
5224: DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5225: DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5226: shift += gsize;
5227: }
5228: ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5229: for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5230: return(0);
5231: }
5235: /*@
5236: DMPlexCheckSymmetry - Check that the adjacency information in the mesh is symmetric.
5238: Input Parameters:
5239: + dm - The DMPlex object
5241: Note: This is a useful diagnostic when creating meshes programmatically.
5243: Level: developer
5245: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5246: @*/
5247: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5248: {
5249: PetscSection coneSection, supportSection;
5250: const PetscInt *cone, *support;
5251: PetscInt coneSize, c, supportSize, s;
5252: PetscInt pStart, pEnd, p, csize, ssize;
5253: PetscErrorCode ierr;
5257: DMPlexGetConeSection(dm, &coneSection);
5258: DMPlexGetSupportSection(dm, &supportSection);
5259: /* Check that point p is found in the support of its cone points, and vice versa */
5260: DMPlexGetChart(dm, &pStart, &pEnd);
5261: for (p = pStart; p < pEnd; ++p) {
5262: DMPlexGetConeSize(dm, p, &coneSize);
5263: DMPlexGetCone(dm, p, &cone);
5264: for (c = 0; c < coneSize; ++c) {
5265: PetscBool dup = PETSC_FALSE;
5266: PetscInt d;
5267: for (d = c-1; d >= 0; --d) {
5268: if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5269: }
5270: DMPlexGetSupportSize(dm, cone[c], &supportSize);
5271: DMPlexGetSupport(dm, cone[c], &support);
5272: for (s = 0; s < supportSize; ++s) {
5273: if (support[s] == p) break;
5274: }
5275: if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5276: PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5277: for (s = 0; s < coneSize; ++s) {
5278: PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
5279: }
5280: PetscPrintf(PETSC_COMM_SELF, "\n");
5281: PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
5282: for (s = 0; s < supportSize; ++s) {
5283: PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
5284: }
5285: PetscPrintf(PETSC_COMM_SELF, "\n");
5286: if (dup) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
5287: else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
5288: }
5289: }
5290: DMPlexGetSupportSize(dm, p, &supportSize);
5291: DMPlexGetSupport(dm, p, &support);
5292: for (s = 0; s < supportSize; ++s) {
5293: DMPlexGetConeSize(dm, support[s], &coneSize);
5294: DMPlexGetCone(dm, support[s], &cone);
5295: for (c = 0; c < coneSize; ++c) {
5296: if (cone[c] == p) break;
5297: }
5298: if (c >= coneSize) {
5299: PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
5300: for (c = 0; c < supportSize; ++c) {
5301: PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
5302: }
5303: PetscPrintf(PETSC_COMM_SELF, "\n");
5304: PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
5305: for (c = 0; c < coneSize; ++c) {
5306: PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
5307: }
5308: PetscPrintf(PETSC_COMM_SELF, "\n");
5309: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
5310: }
5311: }
5312: }
5313: PetscSectionGetStorageSize(coneSection, &csize);
5314: PetscSectionGetStorageSize(supportSection, &ssize);
5315: if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
5316: return(0);
5317: }
5321: /*@
5322: DMPlexCheckSkeleton - Check that each cell has the correct number of vertices
5324: Input Parameters:
5325: + dm - The DMPlex object
5326: . isSimplex - Are the cells simplices or tensor products
5327: - cellHeight - Normally 0
5329: Note: This is a useful diagnostic when creating meshes programmatically.
5331: Level: developer
5333: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
5334: @*/
5335: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5336: {
5337: PetscInt dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;
5342: DMGetDimension(dm, &dim);
5343: switch (dim) {
5344: case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
5345: case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
5346: case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
5347: default:
5348: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
5349: }
5350: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5351: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5352: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5353: cMax = cMax >= 0 ? cMax : cEnd;
5354: for (c = cStart; c < cMax; ++c) {
5355: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
5357: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5358: for (cl = 0; cl < closureSize*2; cl += 2) {
5359: const PetscInt p = closure[cl];
5360: if ((p >= vStart) && (p < vEnd)) ++coneSize;
5361: }
5362: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5363: if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d vertices != %d", c, coneSize, numCorners);
5364: }
5365: for (c = cMax; c < cEnd; ++c) {
5366: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
5368: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5369: for (cl = 0; cl < closureSize*2; cl += 2) {
5370: const PetscInt p = closure[cl];
5371: if ((p >= vStart) && (p < vEnd)) ++coneSize;
5372: }
5373: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5374: if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has %d vertices > %d", c, coneSize, numHybridCorners);
5375: }
5376: return(0);
5377: }
5381: /*@
5382: DMPlexCheckFaces - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type
5384: Input Parameters:
5385: + dm - The DMPlex object
5386: . isSimplex - Are the cells simplices or tensor products
5387: - cellHeight - Normally 0
5389: Note: This is a useful diagnostic when creating meshes programmatically.
5391: Level: developer
5393: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
5394: @*/
5395: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5396: {
5397: PetscInt pMax[4];
5398: PetscInt dim, vStart, vEnd, cStart, cEnd, c, h;
5403: DMGetDimension(dm, &dim);
5404: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5405: DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
5406: for (h = cellHeight; h < dim; ++h) {
5407: DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
5408: for (c = cStart; c < cEnd; ++c) {
5409: const PetscInt *cone, *ornt, *faces;
5410: PetscInt numFaces, faceSize, coneSize,f;
5411: PetscInt *closure = NULL, closureSize, cl, numCorners = 0;
5413: if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
5414: DMPlexGetConeSize(dm, c, &coneSize);
5415: DMPlexGetCone(dm, c, &cone);
5416: DMPlexGetConeOrientation(dm, c, &ornt);
5417: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5418: for (cl = 0; cl < closureSize*2; cl += 2) {
5419: const PetscInt p = closure[cl];
5420: if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
5421: }
5422: DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
5423: if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
5424: for (f = 0; f < numFaces; ++f) {
5425: PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;
5427: DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
5428: for (cl = 0; cl < fclosureSize*2; cl += 2) {
5429: const PetscInt p = fclosure[cl];
5430: if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
5431: }
5432: 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);
5433: for (v = 0; v < fnumCorners; ++v) {
5434: 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]);
5435: }
5436: DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
5437: }
5438: DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
5439: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5440: }
5441: }
5442: return(0);
5443: }
5447: /* Pointwise interpolation
5448: Just code FEM for now
5449: u^f = I u^c
5450: sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
5451: u^f_i = sum_j psi^f_i I phi^c_j u^c_j
5452: I_{ij} = psi^f_i phi^c_j
5453: */
5454: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
5455: {
5456: PetscSection gsc, gsf;
5457: PetscInt m, n;
5458: void *ctx;
5459: DM cdm;
5460: PetscBool regular;
5464: DMGetDefaultGlobalSection(dmFine, &gsf);
5465: PetscSectionGetConstrainedStorageSize(gsf, &m);
5466: DMGetDefaultGlobalSection(dmCoarse, &gsc);
5467: PetscSectionGetConstrainedStorageSize(gsc, &n);
5469: MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
5470: MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
5471: MatSetType(*interpolation, dmCoarse->mattype);
5472: DMGetApplicationContext(dmFine, &ctx);
5474: DMGetCoarseDM(dmFine, &cdm);
5475: DMPlexGetRegularRefinement(dmFine, ®ular);
5476: if (regular && cdm == dmCoarse) {DMPlexComputeInterpolatorNested(dmCoarse, dmFine, *interpolation, ctx);}
5477: else {DMPlexComputeInterpolatorGeneral(dmCoarse, dmFine, *interpolation, ctx);}
5478: MatViewFromOptions(*interpolation, NULL, "-interp_mat_view");
5479: /* Use naive scaling */
5480: DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
5481: return(0);
5482: }
5486: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat)
5487: {
5489: VecScatter ctx;
5492: DMPlexComputeInjectorFEM(dmCoarse, dmFine, &ctx, NULL);
5493: MatCreateScatter(PetscObjectComm((PetscObject)ctx), ctx, mat);
5494: VecScatterDestroy(&ctx);
5495: return(0);
5496: }
5500: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
5501: {
5502: PetscSection section;
5503: IS *bcPoints, *bcComps;
5504: PetscBool *isFE;
5505: PetscInt *bcFields, *numComp, *numDof;
5506: PetscInt depth, dim, numBd, numBC = 0, numFields, bd, bc = 0, f;
5507: PetscInt cStart, cEnd, cEndInterior;
5511: DMGetNumFields(dm, &numFields);
5512: if (!numFields) return(0);
5513: /* FE and FV boundary conditions are handled slightly differently */
5514: PetscMalloc1(numFields, &isFE);
5515: for (f = 0; f < numFields; ++f) {
5516: PetscObject obj;
5517: PetscClassId id;
5519: DMGetField(dm, f, &obj);
5520: PetscObjectGetClassId(obj, &id);
5521: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
5522: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
5523: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
5524: }
5525: /* Allocate boundary point storage for FEM boundaries */
5526: DMPlexGetDepth(dm, &depth);
5527: DMGetDimension(dm, &dim);
5528: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
5529: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
5530: DMGetNumBoundary(dm, &numBd);
5531: for (bd = 0; bd < numBd; ++bd) {
5532: PetscInt field;
5533: PetscBool isEssential;
5535: DMGetBoundary(dm, bd, &isEssential, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL, NULL);
5536: if (isFE[field] && isEssential) ++numBC;
5537: }
5538: /* Add ghost cell boundaries for FVM */
5539: for (f = 0; f < numFields; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
5540: PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);
5541: /* Constrain ghost cells for FV */
5542: for (f = 0; f < numFields; ++f) {
5543: PetscInt *newidx, c;
5545: if (isFE[f] || cEndInterior < 0) continue;
5546: PetscMalloc1(cEnd-cEndInterior,&newidx);
5547: for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
5548: bcFields[bc] = f;
5549: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5550: }
5551: /* Handle FEM Dirichlet boundaries */
5552: for (bd = 0; bd < numBd; ++bd) {
5553: const char *bdLabel;
5554: DMLabel label;
5555: const PetscInt *comps;
5556: const PetscInt *values;
5557: PetscInt bd2, field, numComps, numValues;
5558: PetscBool isEssential, duplicate = PETSC_FALSE;
5560: DMGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
5561: if (!isFE[field]) continue;
5562: DMGetLabel(dm, bdLabel, &label);
5563: /* Only want to modify label once */
5564: for (bd2 = 0; bd2 < bd; ++bd2) {
5565: const char *bdname;
5566: DMGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
5567: PetscStrcmp(bdname, bdLabel, &duplicate);
5568: if (duplicate) break;
5569: }
5570: if (!duplicate && (isFE[field])) {
5571: /* don't complete cells, which are just present to give orientation to the boundary */
5572: DMPlexLabelComplete(dm, label);
5573: }
5574: /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
5575: if (isEssential) {
5576: PetscInt *newidx;
5577: PetscInt n, newn = 0, p, v;
5579: bcFields[bc] = field;
5580: if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
5581: for (v = 0; v < numValues; ++v) {
5582: IS tmp;
5583: const PetscInt *idx;
5585: DMGetStratumIS(dm, bdLabel, values[v], &tmp);
5586: if (!tmp) continue;
5587: ISGetLocalSize(tmp, &n);
5588: ISGetIndices(tmp, &idx);
5589: if (isFE[field]) {
5590: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
5591: } else {
5592: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
5593: }
5594: ISRestoreIndices(tmp, &idx);
5595: ISDestroy(&tmp);
5596: }
5597: PetscMalloc1(newn,&newidx);
5598: newn = 0;
5599: for (v = 0; v < numValues; ++v) {
5600: IS tmp;
5601: const PetscInt *idx;
5603: DMGetStratumIS(dm, bdLabel, values[v], &tmp);
5604: if (!tmp) continue;
5605: ISGetLocalSize(tmp, &n);
5606: ISGetIndices(tmp, &idx);
5607: if (isFE[field]) {
5608: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
5609: } else {
5610: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
5611: }
5612: ISRestoreIndices(tmp, &idx);
5613: ISDestroy(&tmp);
5614: }
5615: ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5616: }
5617: }
5618: /* Handle discretization */
5619: PetscCalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
5620: for (f = 0; f < numFields; ++f) {
5621: PetscObject obj;
5623: DMGetField(dm, f, &obj);
5624: if (isFE[f]) {
5625: PetscFE fe = (PetscFE) obj;
5626: const PetscInt *numFieldDof;
5627: PetscInt d;
5629: PetscFEGetNumComponents(fe, &numComp[f]);
5630: PetscFEGetNumDof(fe, &numFieldDof);
5631: for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
5632: } else {
5633: PetscFV fv = (PetscFV) obj;
5635: PetscFVGetNumComponents(fv, &numComp[f]);
5636: numDof[f*(dim+1)+dim] = numComp[f];
5637: }
5638: }
5639: for (f = 0; f < numFields; ++f) {
5640: PetscInt d;
5641: for (d = 1; d < dim; ++d) {
5642: 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.");
5643: }
5644: }
5645: DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, §ion);
5646: for (f = 0; f < numFields; ++f) {
5647: PetscFE fe;
5648: const char *name;
5650: DMGetField(dm, f, (PetscObject *) &fe);
5651: PetscObjectGetName((PetscObject) fe, &name);
5652: PetscSectionSetFieldName(section, f, name);
5653: }
5654: DMSetDefaultSection(dm, section);
5655: PetscSectionDestroy(§ion);
5656: for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
5657: PetscFree3(bcFields,bcPoints,bcComps);
5658: PetscFree2(numComp,numDof);
5659: PetscFree(isFE);
5660: return(0);
5661: }
5665: /*@
5666: DMPlexGetRegularRefinement - Get the flag indicating that this mesh was obtained by regular refinement from its coarse mesh
5668: Input Parameter:
5669: . dm - The DMPlex object
5671: Output Parameter:
5672: . regular - The flag
5674: Level: intermediate
5676: .seealso: DMPlexSetRegularRefinement()
5677: @*/
5678: PetscErrorCode DMPlexGetRegularRefinement(DM dm, PetscBool *regular)
5679: {
5683: *regular = ((DM_Plex *) dm->data)->regularRefinement;
5684: return(0);
5685: }
5689: /*@
5690: DMPlexSetRegularRefinement - Set the flag indicating that this mesh was obtained by regular refinement from its coarse mesh
5692: Input Parameters:
5693: + dm - The DMPlex object
5694: - regular - The flag
5696: Level: intermediate
5698: .seealso: DMPlexGetRegularRefinement()
5699: @*/
5700: PetscErrorCode DMPlexSetRegularRefinement(DM dm, PetscBool regular)
5701: {
5704: ((DM_Plex *) dm->data)->regularRefinement = regular;
5705: return(0);
5706: }
5708: /* anchors */
5711: /*@
5712: DMPlexGetAnchors - Get the layout of the anchor (point-to-point) constraints. Typically, the user will not have to
5713: call DMPlexGetAnchors() directly: if there are anchors, then DMPlexGetAnchors() is called during DMGetConstraints().
5715: not collective
5717: Input Parameters:
5718: . dm - The DMPlex object
5720: Output Parameters:
5721: + anchorSection - If not NULL, set to the section describing which points anchor the constrained points.
5722: - anchorIS - If not NULL, set to the list of anchors indexed by anchorSection
5725: Level: intermediate
5727: .seealso: DMPlexSetAnchors(), DMGetConstraints(), DMSetConstraints()
5728: @*/
5729: PetscErrorCode DMPlexGetAnchors(DM dm, PetscSection *anchorSection, IS *anchorIS)
5730: {
5731: DM_Plex *plex = (DM_Plex *)dm->data;
5736: if (!plex->anchorSection && !plex->anchorIS && plex->createanchors) {(*plex->createanchors)(dm);}
5737: if (anchorSection) *anchorSection = plex->anchorSection;
5738: if (anchorIS) *anchorIS = plex->anchorIS;
5739: return(0);
5740: }
5744: /*@
5745: DMPlexSetAnchors - Set the layout of the local anchor (point-to-point) constraints. Unlike boundary conditions,
5746: when a point's degrees of freedom in a section are constrained to an outside value, the anchor constraints set a
5747: point's degrees of freedom to be a linear combination of other points' degrees of freedom.
5749: After specifying the layout of constraints with DMPlexSetAnchors(), one specifies the constraints by calling
5750: DMGetConstraints() and filling in the entries in the constraint matrix.
5752: collective on dm
5754: Input Parameters:
5755: + dm - The DMPlex object
5756: . 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).
5757: - anchorIS - The list of all anchor points. Must have a local communicator (PETSC_COMM_SELF or derivative).
5759: The reference counts of anchorSection and anchorIS are incremented.
5761: Level: intermediate
5763: .seealso: DMPlexGetAnchors(), DMGetConstraints(), DMSetConstraints()
5764: @*/
5765: PetscErrorCode DMPlexSetAnchors(DM dm, PetscSection anchorSection, IS anchorIS)
5766: {
5767: DM_Plex *plex = (DM_Plex *)dm->data;
5768: PetscMPIInt result;
5773: if (anchorSection) {
5775: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorSection),&result);
5776: if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor section must have local communicator");
5777: }
5778: if (anchorIS) {
5780: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorIS),&result);
5781: if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor IS must have local communicator");
5782: }
5784: PetscObjectReference((PetscObject)anchorSection);
5785: PetscSectionDestroy(&plex->anchorSection);
5786: plex->anchorSection = anchorSection;
5788: PetscObjectReference((PetscObject)anchorIS);
5789: ISDestroy(&plex->anchorIS);
5790: plex->anchorIS = anchorIS;
5792: #if defined(PETSC_USE_DEBUG)
5793: if (anchorIS && anchorSection) {
5794: PetscInt size, a, pStart, pEnd;
5795: const PetscInt *anchors;
5797: PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5798: ISGetLocalSize(anchorIS,&size);
5799: ISGetIndices(anchorIS,&anchors);
5800: for (a = 0; a < size; a++) {
5801: PetscInt p;
5803: p = anchors[a];
5804: if (p >= pStart && p < pEnd) {
5805: PetscInt dof;
5807: PetscSectionGetDof(anchorSection,p,&dof);
5808: if (dof) {
5809: PetscErrorCode ierr2;
5811: ierr2 = ISRestoreIndices(anchorIS,&anchors);CHKERRQ(ierr2);
5812: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point %d cannot be constrained and an anchor",p);
5813: }
5814: }
5815: }
5816: ISRestoreIndices(anchorIS,&anchors);
5817: }
5818: #endif
5819: /* reset the generic constraints */
5820: DMSetDefaultConstraints(dm,NULL,NULL);
5821: return(0);
5822: }
5826: static PetscErrorCode DMPlexCreateConstraintSection_Anchors(DM dm, PetscSection section, PetscSection *cSec)
5827: {
5828: PetscSection anchorSection;
5829: PetscInt pStart, pEnd, sStart, sEnd, p, dof, numFields, f;
5834: DMPlexGetAnchors(dm,&anchorSection,NULL);
5835: PetscSectionCreate(PETSC_COMM_SELF,cSec);
5836: PetscSectionGetNumFields(section,&numFields);
5837: if (numFields) {
5838: PetscInt f;
5839: PetscSectionSetNumFields(*cSec,numFields);
5841: for (f = 0; f < numFields; f++) {
5842: PetscInt numComp;
5844: PetscSectionGetFieldComponents(section,f,&numComp);
5845: PetscSectionSetFieldComponents(*cSec,f,numComp);
5846: }
5847: }
5848: PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5849: PetscSectionGetChart(section,&sStart,&sEnd);
5850: pStart = PetscMax(pStart,sStart);
5851: pEnd = PetscMin(pEnd,sEnd);
5852: pEnd = PetscMax(pStart,pEnd);
5853: PetscSectionSetChart(*cSec,pStart,pEnd);
5854: for (p = pStart; p < pEnd; p++) {
5855: PetscSectionGetDof(anchorSection,p,&dof);
5856: if (dof) {
5857: PetscSectionGetDof(section,p,&dof);
5858: PetscSectionSetDof(*cSec,p,dof);
5859: for (f = 0; f < numFields; f++) {
5860: PetscSectionGetFieldDof(section,p,f,&dof);
5861: PetscSectionSetFieldDof(*cSec,p,f,dof);
5862: }
5863: }
5864: }
5865: PetscSectionSetUp(*cSec);
5866: return(0);
5867: }
5871: static PetscErrorCode DMPlexCreateConstraintMatrix_Anchors(DM dm, PetscSection section, PetscSection cSec, Mat *cMat)
5872: {
5873: PetscSection aSec;
5874: PetscInt pStart, pEnd, p, dof, aDof, aOff, off, nnz, annz, m, n, q, a, offset, *i, *j;
5875: const PetscInt *anchors;
5876: PetscInt numFields, f;
5877: IS aIS;
5882: PetscSectionGetStorageSize(cSec, &m);
5883: PetscSectionGetStorageSize(section, &n);
5884: MatCreate(PETSC_COMM_SELF,cMat);
5885: MatSetSizes(*cMat,m,n,m,n);
5886: MatSetType(*cMat,MATSEQAIJ);
5887: DMPlexGetAnchors(dm,&aSec,&aIS);
5888: ISGetIndices(aIS,&anchors);
5889: /* cSec will be a subset of aSec and section */
5890: PetscSectionGetChart(cSec,&pStart,&pEnd);
5891: PetscMalloc1(m+1,&i);
5892: i[0] = 0;
5893: PetscSectionGetNumFields(section,&numFields);
5894: for (p = pStart; p < pEnd; p++) {
5895: PetscInt rDof, rOff, r;
5897: PetscSectionGetDof(aSec,p,&rDof);
5898: if (!rDof) continue;
5899: PetscSectionGetOffset(aSec,p,&rOff);
5900: if (numFields) {
5901: for (f = 0; f < numFields; f++) {
5902: annz = 0;
5903: for (r = 0; r < rDof; r++) {
5904: a = anchors[rOff + r];
5905: PetscSectionGetFieldDof(section,a,f,&aDof);
5906: annz += aDof;
5907: }
5908: PetscSectionGetFieldDof(cSec,p,f,&dof);
5909: PetscSectionGetFieldOffset(cSec,p,f,&off);
5910: for (q = 0; q < dof; q++) {
5911: i[off + q + 1] = i[off + q] + annz;
5912: }
5913: }
5914: }
5915: else {
5916: annz = 0;
5917: for (q = 0; q < dof; q++) {
5918: a = anchors[off + q];
5919: PetscSectionGetDof(section,a,&aDof);
5920: annz += aDof;
5921: }
5922: PetscSectionGetDof(cSec,p,&dof);
5923: PetscSectionGetOffset(cSec,p,&off);
5924: for (q = 0; q < dof; q++) {
5925: i[off + q + 1] = i[off + q] + annz;
5926: }
5927: }
5928: }
5929: nnz = i[m];
5930: PetscMalloc1(nnz,&j);
5931: offset = 0;
5932: for (p = pStart; p < pEnd; p++) {
5933: if (numFields) {
5934: for (f = 0; f < numFields; f++) {
5935: PetscSectionGetFieldDof(cSec,p,f,&dof);
5936: for (q = 0; q < dof; q++) {
5937: PetscInt rDof, rOff, r;
5938: PetscSectionGetDof(aSec,p,&rDof);
5939: PetscSectionGetOffset(aSec,p,&rOff);
5940: for (r = 0; r < rDof; r++) {
5941: PetscInt s;
5943: a = anchors[rOff + r];
5944: PetscSectionGetFieldDof(section,a,f,&aDof);
5945: PetscSectionGetFieldOffset(section,a,f,&aOff);
5946: for (s = 0; s < aDof; s++) {
5947: j[offset++] = aOff + s;
5948: }
5949: }
5950: }
5951: }
5952: }
5953: else {
5954: PetscSectionGetDof(cSec,p,&dof);
5955: for (q = 0; q < dof; q++) {
5956: PetscInt rDof, rOff, r;
5957: PetscSectionGetDof(aSec,p,&rDof);
5958: PetscSectionGetOffset(aSec,p,&rOff);
5959: for (r = 0; r < rDof; r++) {
5960: PetscInt s;
5962: a = anchors[rOff + r];
5963: PetscSectionGetDof(section,a,&aDof);
5964: PetscSectionGetOffset(section,a,&aOff);
5965: for (s = 0; s < aDof; s++) {
5966: j[offset++] = aOff + s;
5967: }
5968: }
5969: }
5970: }
5971: }
5972: MatSeqAIJSetPreallocationCSR(*cMat,i,j,NULL);
5973: PetscFree(i);
5974: PetscFree(j);
5975: ISRestoreIndices(aIS,&anchors);
5976: return(0);
5977: }
5981: PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm)
5982: {
5983: DM_Plex *plex = (DM_Plex *)dm->data;
5984: PetscSection anchorSection, section, cSec;
5985: Mat cMat;
5990: DMPlexGetAnchors(dm,&anchorSection,NULL);
5991: if (anchorSection) {
5992: PetscDS ds;
5993: PetscInt nf;
5995: DMGetDefaultSection(dm,§ion);
5996: DMPlexCreateConstraintSection_Anchors(dm,section,&cSec);
5997: DMPlexCreateConstraintMatrix_Anchors(dm,section,cSec,&cMat);
5998: DMGetDS(dm,&ds);
5999: PetscDSGetNumFields(ds,&nf);
6000: if (nf && plex->computeanchormatrix) {(*plex->computeanchormatrix)(dm,section,cSec,cMat);}
6001: DMSetDefaultConstraints(dm,cSec,cMat);
6002: PetscSectionDestroy(&cSec);
6003: MatDestroy(&cMat);
6004: }
6005: return(0);
6006: }