Actual source code: plex.c
petsc-3.7.3 2016-08-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) dm); /* viewer drops reference */
84: PetscObjectReference((PetscObject) v); /* viewer drops reference */
85: PetscViewerVTKAddField(viewer, (PetscObject) dm, DMPlexVTKWriteAll, ft, (PetscObject) v);
86: } else if (ishdf5) {
87: #if defined(PETSC_HAVE_HDF5)
88: VecView_Plex_Local_HDF5(v, viewer);
89: #else
90: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
91: #endif
92: } else {
93: if (isseq) {VecView_Seq(v, viewer);}
94: else {VecView_MPI(v, viewer);}
95: }
96: return(0);
97: }
101: PetscErrorCode VecView_Plex(Vec v, PetscViewer viewer)
102: {
103: DM dm;
104: PetscBool isvtk, ishdf5, isseq;
108: VecGetDM(v, &dm);
109: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
110: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
111: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
112: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
113: if (isvtk) {
114: Vec locv;
115: const char *name;
117: DMGetLocalVector(dm, &locv);
118: PetscObjectGetName((PetscObject) v, &name);
119: PetscObjectSetName((PetscObject) locv, name);
120: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
121: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
122: VecView_Plex_Local(locv, viewer);
123: DMRestoreLocalVector(dm, &locv);
124: } else if (ishdf5) {
125: #if defined(PETSC_HAVE_HDF5)
126: VecView_Plex_HDF5(v, viewer);
127: #else
128: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
129: #endif
130: } else {
131: if (isseq) {VecView_Seq(v, viewer);}
132: else {VecView_MPI(v, viewer);}
133: }
134: return(0);
135: }
139: PetscErrorCode VecView_Plex_Native(Vec originalv, PetscViewer viewer)
140: {
141: DM dm;
142: MPI_Comm comm;
143: PetscViewerFormat format;
144: Vec v;
145: PetscBool isvtk, ishdf5;
146: PetscErrorCode ierr;
149: VecGetDM(originalv, &dm);
150: PetscObjectGetComm((PetscObject) originalv, &comm);
151: if (!dm) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
152: PetscViewerGetFormat(viewer, &format);
153: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
154: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
155: if (format == PETSC_VIEWER_NATIVE) {
156: const char *vecname;
157: PetscInt n, nroots;
159: if (dm->sfNatural) {
160: VecGetLocalSize(originalv, &n);
161: PetscSFGetGraph(dm->sfNatural, &nroots, NULL, NULL, NULL);
162: if (n == nroots) {
163: DMGetGlobalVector(dm, &v);
164: DMPlexGlobalToNaturalBegin(dm, originalv, v);
165: DMPlexGlobalToNaturalEnd(dm, originalv, v);
166: PetscObjectGetName((PetscObject) originalv, &vecname);
167: PetscObjectSetName((PetscObject) v, vecname);
168: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "DM global to natural SF only handles global vectors");
169: } else SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created");
170: } else {
171: /* we are viewing a natural DMPlex vec. */
172: v = originalv;
173: }
174: if (ishdf5) {
175: #if defined(PETSC_HAVE_HDF5)
176: VecView_Plex_HDF5_Native(v, viewer);
177: #else
178: SETERRQ(comm, PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
179: #endif
180: } else if (isvtk) {
181: SETERRQ(comm, PETSC_ERR_SUP, "VTK format does not support viewing in natural order. Please switch to HDF5.");
182: } else {
183: PetscBool isseq;
185: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
186: if (isseq) {VecView_Seq(v, viewer);}
187: else {VecView_MPI(v, viewer);}
188: }
189: if (format == PETSC_VIEWER_NATIVE) {DMRestoreGlobalVector(dm, &v);}
190: return(0);
191: }
195: PetscErrorCode VecLoad_Plex_Local(Vec v, PetscViewer viewer)
196: {
197: DM dm;
198: PetscBool ishdf5;
202: VecGetDM(v, &dm);
203: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
204: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
205: if (ishdf5) {
206: DM dmBC;
207: Vec gv;
208: const char *name;
210: DMGetOutputDM(dm, &dmBC);
211: DMGetGlobalVector(dmBC, &gv);
212: PetscObjectGetName((PetscObject) v, &name);
213: PetscObjectSetName((PetscObject) gv, name);
214: VecLoad_Default(gv, viewer);
215: DMGlobalToLocalBegin(dmBC, gv, INSERT_VALUES, v);
216: DMGlobalToLocalEnd(dmBC, gv, INSERT_VALUES, v);
217: DMRestoreGlobalVector(dmBC, &gv);
218: } else {
219: VecLoad_Default(v, viewer);
220: }
221: return(0);
222: }
226: PetscErrorCode VecLoad_Plex(Vec v, PetscViewer viewer)
227: {
228: DM dm;
229: PetscBool ishdf5;
233: VecGetDM(v, &dm);
234: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
235: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
236: if (ishdf5) {
237: #if defined(PETSC_HAVE_HDF5)
238: VecLoad_Plex_HDF5(v, viewer);
239: #else
240: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
241: #endif
242: } else {
243: VecLoad_Default(v, viewer);
244: }
245: return(0);
246: }
250: PetscErrorCode VecLoad_Plex_Native(Vec originalv, PetscViewer viewer)
251: {
252: DM dm;
253: PetscViewerFormat format;
254: PetscBool ishdf5;
255: PetscErrorCode ierr;
258: VecGetDM(originalv, &dm);
259: if (!dm) SETERRQ(PetscObjectComm((PetscObject) originalv), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
260: PetscViewerGetFormat(viewer, &format);
261: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
262: if (format == PETSC_VIEWER_NATIVE) {
263: if (dm->sfNatural) {
264: if (ishdf5) {
265: #if defined(PETSC_HAVE_HDF5)
266: Vec v;
267: const char *vecname;
269: DMGetGlobalVector(dm, &v);
270: PetscObjectGetName((PetscObject) originalv, &vecname);
271: PetscObjectSetName((PetscObject) v, vecname);
272: VecLoad_Plex_HDF5_Native(v, viewer);
273: DMPlexNaturalToGlobalBegin(dm, v, originalv);
274: DMPlexNaturalToGlobalEnd(dm, v, originalv);
275: DMRestoreGlobalVector(dm, &v);
276: #else
277: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
278: #endif
279: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Reading in natural order is not supported for anything but HDF5.");
280: }
281: }
282: return(0);
283: }
287: PetscErrorCode DMPlexView_Ascii_Geometry(DM dm, PetscViewer viewer)
288: {
289: PetscSection coordSection;
290: Vec coordinates;
291: DMLabel depthLabel;
292: const char *name[4];
293: const PetscScalar *a;
294: PetscInt dim, pStart, pEnd, cStart, cEnd, c;
295: PetscErrorCode ierr;
298: DMGetDimension(dm, &dim);
299: DMGetCoordinatesLocal(dm, &coordinates);
300: DMGetCoordinateSection(dm, &coordSection);
301: DMPlexGetDepthLabel(dm, &depthLabel);
302: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
303: PetscSectionGetChart(coordSection, &pStart, &pEnd);
304: VecGetArrayRead(coordinates, &a);
305: name[0] = "vertex";
306: name[1] = "edge";
307: name[dim-1] = "face";
308: name[dim] = "cell";
309: for (c = cStart; c < cEnd; ++c) {
310: PetscInt *closure = NULL;
311: PetscInt closureSize, cl;
313: PetscViewerASCIIPrintf(viewer, "Geometry for cell %D:\n", c);
314: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
315: PetscViewerASCIIPushTab(viewer);
316: for (cl = 0; cl < closureSize*2; cl += 2) {
317: PetscInt point = closure[cl], depth, dof, off, d, p;
319: if ((point < pStart) || (point >= pEnd)) continue;
320: PetscSectionGetDof(coordSection, point, &dof);
321: if (!dof) continue;
322: DMLabelGetValue(depthLabel, point, &depth);
323: PetscSectionGetOffset(coordSection, point, &off);
324: PetscViewerASCIIPrintf(viewer, "%s %D coords:", name[depth], point);
325: for (p = 0; p < dof/dim; ++p) {
326: PetscViewerASCIIPrintf(viewer, " (");
327: for (d = 0; d < dim; ++d) {
328: if (d > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
329: PetscViewerASCIIPrintf(viewer, "%g", PetscRealPart(a[off+p*dim+d]));
330: }
331: PetscViewerASCIIPrintf(viewer, ")");
332: }
333: PetscViewerASCIIPrintf(viewer, "\n");
334: }
335: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
336: PetscViewerASCIIPopTab(viewer);
337: }
338: VecRestoreArrayRead(coordinates, &a);
339: return(0);
340: }
344: PetscErrorCode DMPlexView_Ascii(DM dm, PetscViewer viewer)
345: {
346: DM_Plex *mesh = (DM_Plex*) dm->data;
347: DM cdm;
348: DMLabel markers;
349: PetscSection coordSection;
350: Vec coordinates;
351: PetscViewerFormat format;
352: PetscErrorCode ierr;
355: DMGetCoordinateDM(dm, &cdm);
356: DMGetDefaultSection(cdm, &coordSection);
357: DMGetCoordinatesLocal(dm, &coordinates);
358: PetscViewerGetFormat(viewer, &format);
359: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
360: const char *name;
361: PetscInt maxConeSize, maxSupportSize;
362: PetscInt pStart, pEnd, p;
363: PetscMPIInt rank, size;
365: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
366: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
367: PetscObjectGetName((PetscObject) dm, &name);
368: DMPlexGetChart(dm, &pStart, &pEnd);
369: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
370: PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
371: PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
372: PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
373: PetscViewerASCIIPushSynchronized(viewer);
374: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Max sizes cone: %D support: %D\n", rank,maxConeSize, maxSupportSize);
375: for (p = pStart; p < pEnd; ++p) {
376: PetscInt dof, off, s;
378: PetscSectionGetDof(mesh->supportSection, p, &dof);
379: PetscSectionGetOffset(mesh->supportSection, p, &off);
380: for (s = off; s < off+dof; ++s) {
381: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D ----> %D\n", rank, p, mesh->supports[s]);
382: }
383: }
384: PetscViewerFlush(viewer);
385: PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
386: for (p = pStart; p < pEnd; ++p) {
387: PetscInt dof, off, c;
389: PetscSectionGetDof(mesh->coneSection, p, &dof);
390: PetscSectionGetOffset(mesh->coneSection, p, &off);
391: for (c = off; c < off+dof; ++c) {
392: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D <---- %D (%D)\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]);
393: }
394: }
395: PetscViewerFlush(viewer);
396: PetscViewerASCIIPopSynchronized(viewer);
397: PetscSectionGetChart(coordSection, &pStart, NULL);
398: if (pStart >= 0) {PetscSectionVecView(coordSection, coordinates, viewer);}
399: DMGetLabel(dm, "marker", &markers);
400: DMLabelView(markers,viewer);
401: if (size > 1) {
402: PetscSF sf;
404: DMGetPointSF(dm, &sf);
405: PetscSFView(sf, viewer);
406: }
407: PetscViewerFlush(viewer);
408: } else if (format == PETSC_VIEWER_ASCII_LATEX) {
409: const char *name, *color;
410: const char *defcolors[3] = {"gray", "orange", "green"};
411: const char *deflcolors[4] = {"blue", "cyan", "red", "magenta"};
412: PetscReal scale = 2.0;
413: PetscBool useNumbers = PETSC_TRUE, useLabels, useColors;
414: double tcoords[3];
415: PetscScalar *coords;
416: PetscInt numLabels, l, numColors, numLColors, dim, depth, cStart, cEnd, c, vStart, vEnd, v, eStart = 0, eEnd = 0, e, p;
417: PetscMPIInt rank, size;
418: char **names, **colors, **lcolors;
420: DMGetDimension(dm, &dim);
421: DMPlexGetDepth(dm, &depth);
422: DMGetNumLabels(dm, &numLabels);
423: numLabels = PetscMax(numLabels, 10);
424: numColors = 10;
425: numLColors = 10;
426: PetscCalloc3(numLabels, &names, numColors, &colors, numLColors, &lcolors);
427: PetscOptionsGetReal(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_scale", &scale, NULL);
428: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_numbers", &useNumbers, NULL);
429: PetscOptionsGetStringArray(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_labels", names, &numLabels, &useLabels);
430: if (!useLabels) numLabels = 0;
431: PetscOptionsGetStringArray(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_colors", colors, &numColors, &useColors);
432: if (!useColors) {
433: numColors = 3;
434: for (c = 0; c < numColors; ++c) {PetscStrallocpy(defcolors[c], &colors[c]);}
435: }
436: PetscOptionsGetStringArray(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_view_lcolors", lcolors, &numLColors, &useColors);
437: if (!useColors) {
438: numLColors = 4;
439: for (c = 0; c < numLColors; ++c) {PetscStrallocpy(deflcolors[c], &lcolors[c]);}
440: }
441: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
442: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
443: PetscObjectGetName((PetscObject) dm, &name);
444: PetscViewerASCIIPrintf(viewer, "\
445: \\documentclass[tikz]{standalone}\n\n\
446: \\usepackage{pgflibraryshapes}\n\
447: \\usetikzlibrary{backgrounds}\n\
448: \\usetikzlibrary{arrows}\n\
449: \\begin{document}\n");
450: if (size > 1) {
451: PetscViewerASCIIPrintf(viewer, "%s for process ", name);
452: for (p = 0; p < size; ++p) {
453: if (p > 0 && p == size-1) {
454: PetscViewerASCIIPrintf(viewer, ", and ", colors[p%numColors], p);
455: } else if (p > 0) {
456: PetscViewerASCIIPrintf(viewer, ", ", colors[p%numColors], p);
457: }
458: PetscViewerASCIIPrintf(viewer, "{\\textcolor{%s}%D}", colors[p%numColors], p);
459: }
460: PetscViewerASCIIPrintf(viewer, ".\n\n\n");
461: }
462: PetscViewerASCIIPrintf(viewer, "\\begin{tikzpicture}[scale = %g,font=\\fontsize{8}{8}\\selectfont]\n", 1.0);
463: /* Plot vertices */
464: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
465: VecGetArray(coordinates, &coords);
466: PetscViewerASCIIPushSynchronized(viewer);
467: for (v = vStart; v < vEnd; ++v) {
468: PetscInt off, dof, d;
469: PetscBool isLabeled = PETSC_FALSE;
471: PetscSectionGetDof(coordSection, v, &dof);
472: PetscSectionGetOffset(coordSection, v, &off);
473: PetscViewerASCIISynchronizedPrintf(viewer, "\\path (");
474: if (PetscUnlikely(dof > 3)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"coordSection vertex %D has dof %D > 3",v,dof);
475: for (d = 0; d < dof; ++d) {
476: tcoords[d] = (double) (scale*PetscRealPart(coords[off+d]));
477: tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
478: }
479: /* Rotate coordinates since PGF makes z point out of the page instead of up */
480: if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
481: for (d = 0; d < dof; ++d) {
482: if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
483: PetscViewerASCIISynchronizedPrintf(viewer, "%g", tcoords[d]);
484: }
485: color = colors[rank%numColors];
486: for (l = 0; l < numLabels; ++l) {
487: PetscInt val;
488: DMGetLabelValue(dm, names[l], v, &val);
489: if (val >= 0) {color = lcolors[l%numLColors]; isLabeled = PETSC_TRUE; break;}
490: }
491: if (useNumbers) {
492: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [draw,shape=circle,color=%s] {%D};\n", v, rank, color, v);
493: } else {
494: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [fill,inner sep=%dpt,shape=circle,color=%s] {};\n", v, rank, !isLabeled ? 1 : 2, color);
495: }
496: }
497: VecRestoreArray(coordinates, &coords);
498: PetscViewerFlush(viewer);
499: /* Plot edges */
500: if (depth > 1) {DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);}
501: if (dim < 3 && useNumbers) {
502: VecGetArray(coordinates, &coords);
503: PetscViewerASCIIPrintf(viewer, "\\path\n");
504: for (e = eStart; e < eEnd; ++e) {
505: const PetscInt *cone;
506: PetscInt coneSize, offA, offB, dof, d;
508: DMPlexGetConeSize(dm, e, &coneSize);
509: if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %D cone should have two vertices, not %D", e, coneSize);
510: DMPlexGetCone(dm, e, &cone);
511: PetscSectionGetDof(coordSection, cone[0], &dof);
512: PetscSectionGetOffset(coordSection, cone[0], &offA);
513: PetscSectionGetOffset(coordSection, cone[1], &offB);
514: PetscViewerASCIISynchronizedPrintf(viewer, "(");
515: for (d = 0; d < dof; ++d) {
516: tcoords[d] = (double) (scale*PetscRealPart(coords[offA+d]+coords[offB+d]));
517: tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
518: }
519: /* Rotate coordinates since PGF makes z point out of the page instead of up */
520: if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
521: for (d = 0; d < dof; ++d) {
522: if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
523: PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)tcoords[d]);
524: }
525: color = colors[rank%numColors];
526: for (l = 0; l < numLabels; ++l) {
527: PetscInt val;
528: DMGetLabelValue(dm, names[l], v, &val);
529: if (val >= 0) {color = lcolors[l%numLColors]; break;}
530: }
531: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [draw,shape=circle,color=%s] {%D} --\n", e, rank, color, e);
532: }
533: VecRestoreArray(coordinates, &coords);
534: PetscViewerFlush(viewer);
535: PetscViewerASCIIPrintf(viewer, "(0,0);\n");
536: }
537: /* Plot cells */
538: if (dim == 3 || !useNumbers) {
539: for (e = eStart; e < eEnd; ++e) {
540: const PetscInt *cone;
542: color = colors[rank%numColors];
543: for (l = 0; l < numLabels; ++l) {
544: PetscInt val;
545: DMGetLabelValue(dm, names[l], e, &val);
546: if (val >= 0) {color = lcolors[l%numLColors]; break;}
547: }
548: DMPlexGetCone(dm, e, &cone);
549: PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] (%D_%d) -- (%D_%d);\n", color, cone[0], rank, cone[1], rank);
550: }
551: } else {
552: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
553: for (c = cStart; c < cEnd; ++c) {
554: PetscInt *closure = NULL;
555: PetscInt closureSize, firstPoint = -1;
557: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
558: PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] ", colors[rank%numColors]);
559: for (p = 0; p < closureSize*2; p += 2) {
560: const PetscInt point = closure[p];
562: if ((point < vStart) || (point >= vEnd)) continue;
563: if (firstPoint >= 0) {PetscViewerASCIISynchronizedPrintf(viewer, " -- ");}
564: PetscViewerASCIISynchronizedPrintf(viewer, "(%D_%d)", point, rank);
565: if (firstPoint < 0) firstPoint = point;
566: }
567: /* Why doesn't this work? PetscViewerASCIISynchronizedPrintf(viewer, " -- cycle;\n"); */
568: PetscViewerASCIISynchronizedPrintf(viewer, " -- (%D_%d);\n", firstPoint, rank);
569: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
570: }
571: }
572: PetscViewerFlush(viewer);
573: PetscViewerASCIIPopSynchronized(viewer);
574: PetscViewerASCIIPrintf(viewer, "\\end{tikzpicture}\n");
575: PetscViewerASCIIPrintf(viewer, "\\end{document}\n", name);
576: for (l = 0; l < numLabels; ++l) {PetscFree(names[l]);}
577: for (c = 0; c < numColors; ++c) {PetscFree(colors[c]);}
578: for (c = 0; c < numLColors; ++c) {PetscFree(lcolors[c]);}
579: PetscFree3(names, colors, lcolors);
580: } else {
581: MPI_Comm comm;
582: PetscInt *sizes, *hybsizes;
583: PetscInt locDepth, depth, dim, d, pMax[4];
584: PetscInt pStart, pEnd, p;
585: PetscInt numLabels, l;
586: const char *name;
587: PetscMPIInt size;
589: PetscObjectGetComm((PetscObject)dm,&comm);
590: MPI_Comm_size(comm, &size);
591: DMGetDimension(dm, &dim);
592: PetscObjectGetName((PetscObject) dm, &name);
593: if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dim);}
594: else {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dim);}
595: DMPlexGetDepth(dm, &locDepth);
596: MPIU_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
597: DMPlexGetHybridBounds(dm, &pMax[depth], depth > 0 ? &pMax[depth-1] : NULL, &pMax[1], &pMax[0]);
598: PetscMalloc2(size,&sizes,size,&hybsizes);
599: if (depth == 1) {
600: DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd);
601: pEnd = pEnd - pStart;
602: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
603: PetscViewerASCIIPrintf(viewer, " %d-cells:", 0);
604: for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
605: PetscViewerASCIIPrintf(viewer, "\n");
606: DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd);
607: pEnd = pEnd - pStart;
608: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
609: PetscViewerASCIIPrintf(viewer, " %D-cells:", dim);
610: for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
611: PetscViewerASCIIPrintf(viewer, "\n");
612: } else {
613: for (d = 0; d <= dim; d++) {
614: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
615: pEnd -= pStart;
616: pMax[d] -= pStart;
617: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
618: MPI_Gather(&pMax[d], 1, MPIU_INT, hybsizes, 1, MPIU_INT, 0, comm);
619: PetscViewerASCIIPrintf(viewer, " %D-cells:", d);
620: for (p = 0; p < size; ++p) {
621: if (hybsizes[p] >= 0) {PetscViewerASCIIPrintf(viewer, " %D (%D)", sizes[p], sizes[p] - hybsizes[p]);}
622: else {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
623: }
624: PetscViewerASCIIPrintf(viewer, "\n");
625: }
626: }
627: PetscFree2(sizes,hybsizes);
628: DMGetNumLabels(dm, &numLabels);
629: if (numLabels) {PetscViewerASCIIPrintf(viewer, "Labels:\n");}
630: for (l = 0; l < numLabels; ++l) {
631: DMLabel label;
632: const char *name;
633: IS valueIS;
634: const PetscInt *values;
635: PetscInt numValues, v;
637: DMGetLabelName(dm, l, &name);
638: DMGetLabel(dm, name, &label);
639: DMLabelGetNumValues(label, &numValues);
640: PetscViewerASCIIPrintf(viewer, " %s: %D strata of sizes (", name, numValues);
641: DMLabelGetValueIS(label, &valueIS);
642: ISGetIndices(valueIS, &values);
643: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
644: for (v = 0; v < numValues; ++v) {
645: PetscInt size;
647: DMLabelGetStratumSize(label, values[v], &size);
648: if (v > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
649: PetscViewerASCIIPrintf(viewer, "%D", size);
650: }
651: PetscViewerASCIIPrintf(viewer, ")\n");
652: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
653: ISRestoreIndices(valueIS, &values);
654: ISDestroy(&valueIS);
655: }
656: DMGetCoarseDM(dm, &cdm);
657: if (cdm) {
658: PetscViewerASCIIPushTab(viewer);
659: DMPlexView_Ascii(cdm, viewer);
660: PetscViewerASCIIPopTab(viewer);
661: }
662: }
663: return(0);
664: }
668: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
669: {
670: PetscBool iascii, ishdf5, isvtk;
676: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
677: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
678: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
679: if (iascii) {
680: DMPlexView_Ascii(dm, viewer);
681: } else if (ishdf5) {
682: #if defined(PETSC_HAVE_HDF5)
683: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
684: DMPlexView_HDF5(dm, viewer);
685: PetscViewerPopFormat(viewer);
686: #else
687: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
688: #endif
689: }
690: else if (isvtk) {
691: DMPlexVTKWriteAll((PetscObject) dm,viewer);
692: }
693: return(0);
694: }
698: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
699: {
700: PetscBool isbinary, ishdf5;
706: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
707: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
708: if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");}
709: else if (ishdf5) {
710: #if defined(PETSC_HAVE_HDF5)
711: DMPlexLoad_HDF5(dm, viewer);
712: #else
713: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
714: #endif
715: }
716: return(0);
717: }
722: PetscErrorCode DMDestroy_Plex(DM dm)
723: {
724: DM_Plex *mesh = (DM_Plex*) dm->data;
728: if (--mesh->refct > 0) return(0);
729: PetscSectionDestroy(&mesh->coneSection);
730: PetscFree(mesh->cones);
731: PetscFree(mesh->coneOrientations);
732: PetscSectionDestroy(&mesh->supportSection);
733: PetscFree(mesh->supports);
734: PetscFree(mesh->facesTmp);
735: PetscFree(mesh->tetgenOpts);
736: PetscFree(mesh->triangleOpts);
737: PetscPartitionerDestroy(&mesh->partitioner);
738: DMLabelDestroy(&mesh->subpointMap);
739: ISDestroy(&mesh->globalVertexNumbers);
740: ISDestroy(&mesh->globalCellNumbers);
741: PetscSectionDestroy(&mesh->anchorSection);
742: ISDestroy(&mesh->anchorIS);
743: PetscSectionDestroy(&mesh->parentSection);
744: PetscFree(mesh->parents);
745: PetscFree(mesh->childIDs);
746: PetscSectionDestroy(&mesh->childSection);
747: PetscFree(mesh->children);
748: DMDestroy(&mesh->referenceTree);
749: PetscGridHashDestroy(&mesh->lbox);
750: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
751: PetscFree(mesh);
752: return(0);
753: }
757: PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
758: {
759: PetscSection sectionGlobal;
760: PetscInt bs = -1;
761: PetscInt localSize;
762: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock;
764: MatType mtype;
765: ISLocalToGlobalMapping ltog;
768: MatInitializePackage();
769: mtype = dm->mattype;
770: DMGetDefaultGlobalSection(dm, §ionGlobal);
771: /* PetscSectionGetStorageSize(sectionGlobal, &localSize); */
772: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
773: MatCreate(PetscObjectComm((PetscObject)dm), J);
774: MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
775: MatSetType(*J, mtype);
776: MatSetFromOptions(*J);
777: PetscStrcmp(mtype, MATSHELL, &isShell);
778: PetscStrcmp(mtype, MATBAIJ, &isBlock);
779: PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
780: PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
781: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
782: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
783: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
784: if (!isShell) {
785: PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
786: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin;
788: if (bs < 0) {
789: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
790: PetscInt pStart, pEnd, p, dof, cdof;
792: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
793: for (p = pStart; p < pEnd; ++p) {
794: PetscInt bdof;
796: PetscSectionGetDof(sectionGlobal, p, &dof);
797: PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
798: bdof = PetscAbs(dof) - cdof;
799: if (bdof) {
800: if (bs < 0) {
801: bs = bdof;
802: } else if (bs != bdof) {
803: /* Layout does not admit a pointwise block size */
804: bs = 1;
805: break;
806: }
807: }
808: }
809: /* Must have same blocksize on all procs (some might have no points) */
810: bsLocal = bs;
811: MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
812: bsLocal = bs < 0 ? bsMax : bs;
813: MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
814: if (bsMin != bsMax) {
815: bs = 1;
816: } else {
817: bs = bsMax;
818: }
819: } else {
820: bs = 1;
821: }
822: }
823: PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
824: DMPlexPreallocateOperator(dm, bs, dnz, onz, dnzu, onzu, *J, fillMatrix);
825: PetscFree4(dnz, onz, dnzu, onzu);
827: /* Set localtoglobalmapping on the matrix for MatSetValuesLocal() to work */
828: DMGetLocalToGlobalMapping(dm,<og);
829: MatSetLocalToGlobalMapping(*J,ltog,ltog);
830: }
831: return(0);
832: }
836: /*@
837: DMPlexGetChart - Return the interval for all mesh points [pStart, pEnd)
839: Not collective
841: Input Parameter:
842: . mesh - The DMPlex
844: Output Parameters:
845: + pStart - The first mesh point
846: - pEnd - The upper bound for mesh points
848: Level: beginner
850: .seealso: DMPlexCreate(), DMPlexSetChart()
851: @*/
852: PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
853: {
854: DM_Plex *mesh = (DM_Plex*) dm->data;
859: PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
860: return(0);
861: }
865: /*@
866: DMPlexSetChart - Set the interval for all mesh points [pStart, pEnd)
868: Not collective
870: Input Parameters:
871: + mesh - The DMPlex
872: . pStart - The first mesh point
873: - pEnd - The upper bound for mesh points
875: Output Parameters:
877: Level: beginner
879: .seealso: DMPlexCreate(), DMPlexGetChart()
880: @*/
881: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
882: {
883: DM_Plex *mesh = (DM_Plex*) dm->data;
888: PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
889: PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
890: return(0);
891: }
895: /*@
896: DMPlexGetConeSize - Return the number of in-edges for this point in the Sieve DAG
898: Not collective
900: Input Parameters:
901: + mesh - The DMPlex
902: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
904: Output Parameter:
905: . size - The cone size for point p
907: Level: beginner
909: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
910: @*/
911: PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size)
912: {
913: DM_Plex *mesh = (DM_Plex*) dm->data;
919: PetscSectionGetDof(mesh->coneSection, p, size);
920: return(0);
921: }
925: /*@
926: DMPlexSetConeSize - Set the number of in-edges for this point in the Sieve DAG
928: Not collective
930: Input Parameters:
931: + mesh - The DMPlex
932: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
933: - size - The cone size for point p
935: Output Parameter:
937: Note:
938: This should be called after DMPlexSetChart().
940: Level: beginner
942: .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart()
943: @*/
944: PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size)
945: {
946: DM_Plex *mesh = (DM_Plex*) dm->data;
951: PetscSectionSetDof(mesh->coneSection, p, size);
953: mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
954: return(0);
955: }
959: /*@
960: DMPlexAddConeSize - Add the given number of in-edges to this point in the Sieve DAG
962: Not collective
964: Input Parameters:
965: + mesh - The DMPlex
966: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
967: - size - The additional cone size for point p
969: Output Parameter:
971: Note:
972: This should be called after DMPlexSetChart().
974: Level: beginner
976: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart()
977: @*/
978: PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size)
979: {
980: DM_Plex *mesh = (DM_Plex*) dm->data;
981: PetscInt csize;
986: PetscSectionAddDof(mesh->coneSection, p, size);
987: PetscSectionGetDof(mesh->coneSection, p, &csize);
989: mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
990: return(0);
991: }
995: /*@C
996: DMPlexGetCone - Return the points on the in-edges for this point in the Sieve DAG
998: Not collective
1000: Input Parameters:
1001: + mesh - The DMPlex
1002: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1004: Output Parameter:
1005: . cone - An array of points which are on the in-edges for point p
1007: Level: beginner
1009: Fortran Notes:
1010: Since it returns an array, this routine is only available in Fortran 90, and you must
1011: include petsc.h90 in your code.
1013: You must also call DMPlexRestoreCone() after you finish using the returned array.
1015: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart()
1016: @*/
1017: PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
1018: {
1019: DM_Plex *mesh = (DM_Plex*) dm->data;
1020: PetscInt off;
1026: PetscSectionGetOffset(mesh->coneSection, p, &off);
1027: *cone = &mesh->cones[off];
1028: return(0);
1029: }
1033: /*@
1034: DMPlexSetCone - Set the points on the in-edges for this point in the Sieve DAG
1036: Not collective
1038: Input Parameters:
1039: + mesh - The DMPlex
1040: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1041: - cone - An array of points which are on the in-edges for point p
1043: Output Parameter:
1045: Note:
1046: This should be called after all calls to DMPlexSetConeSize() and DMSetUp().
1048: Level: beginner
1050: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1051: @*/
1052: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
1053: {
1054: DM_Plex *mesh = (DM_Plex*) dm->data;
1055: PetscInt pStart, pEnd;
1056: PetscInt dof, off, c;
1061: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1062: PetscSectionGetDof(mesh->coneSection, p, &dof);
1064: PetscSectionGetOffset(mesh->coneSection, p, &off);
1065: 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);
1066: for (c = 0; c < dof; ++c) {
1067: 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);
1068: mesh->cones[off+c] = cone[c];
1069: }
1070: return(0);
1071: }
1075: /*@C
1076: DMPlexGetConeOrientation - Return the orientations on the in-edges for this point in the Sieve DAG
1078: Not collective
1080: Input Parameters:
1081: + mesh - The DMPlex
1082: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1084: Output Parameter:
1085: . coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1086: integer giving the prescription for cone traversal. If it is negative, the cone is
1087: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1088: the index of the cone point on which to start.
1090: Level: beginner
1092: Fortran Notes:
1093: Since it returns an array, this routine is only available in Fortran 90, and you must
1094: include petsc.h90 in your code.
1096: You must also call DMPlexRestoreConeOrientation() after you finish using the returned array.
1098: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
1099: @*/
1100: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
1101: {
1102: DM_Plex *mesh = (DM_Plex*) dm->data;
1103: PetscInt off;
1108: #if defined(PETSC_USE_DEBUG)
1109: {
1110: PetscInt dof;
1111: PetscSectionGetDof(mesh->coneSection, p, &dof);
1113: }
1114: #endif
1115: PetscSectionGetOffset(mesh->coneSection, p, &off);
1117: *coneOrientation = &mesh->coneOrientations[off];
1118: return(0);
1119: }
1123: /*@
1124: DMPlexSetConeOrientation - Set the orientations on the in-edges for this point in the Sieve DAG
1126: Not collective
1128: Input Parameters:
1129: + mesh - The DMPlex
1130: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1131: - coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1132: integer giving the prescription for cone traversal. If it is negative, the cone is
1133: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1134: the index of the cone point on which to start.
1136: Output Parameter:
1138: Note:
1139: This should be called after all calls to DMPlexSetConeSize() and DMSetUp().
1141: Level: beginner
1143: .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1144: @*/
1145: PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
1146: {
1147: DM_Plex *mesh = (DM_Plex*) dm->data;
1148: PetscInt pStart, pEnd;
1149: PetscInt dof, off, c;
1154: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1155: PetscSectionGetDof(mesh->coneSection, p, &dof);
1157: PetscSectionGetOffset(mesh->coneSection, p, &off);
1158: 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);
1159: for (c = 0; c < dof; ++c) {
1160: PetscInt cdof, o = coneOrientation[c];
1162: PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1163: 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);
1164: mesh->coneOrientations[off+c] = o;
1165: }
1166: return(0);
1167: }
1171: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
1172: {
1173: DM_Plex *mesh = (DM_Plex*) dm->data;
1174: PetscInt pStart, pEnd;
1175: PetscInt dof, off;
1180: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1181: 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);
1182: 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);
1183: PetscSectionGetDof(mesh->coneSection, p, &dof);
1184: PetscSectionGetOffset(mesh->coneSection, p, &off);
1185: 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);
1186: mesh->cones[off+conePos] = conePoint;
1187: return(0);
1188: }
1192: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1193: {
1194: DM_Plex *mesh = (DM_Plex*) dm->data;
1195: PetscInt pStart, pEnd;
1196: PetscInt dof, off;
1201: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1202: 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);
1203: PetscSectionGetDof(mesh->coneSection, p, &dof);
1204: PetscSectionGetOffset(mesh->coneSection, p, &off);
1205: 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);
1206: mesh->coneOrientations[off+conePos] = coneOrientation;
1207: return(0);
1208: }
1212: /*@
1213: DMPlexGetSupportSize - Return the number of out-edges for this point in the Sieve DAG
1215: Not collective
1217: Input Parameters:
1218: + mesh - The DMPlex
1219: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1221: Output Parameter:
1222: . size - The support size for point p
1224: Level: beginner
1226: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1227: @*/
1228: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1229: {
1230: DM_Plex *mesh = (DM_Plex*) dm->data;
1236: PetscSectionGetDof(mesh->supportSection, p, size);
1237: return(0);
1238: }
1242: /*@
1243: DMPlexSetSupportSize - Set the number of out-edges for this point in the Sieve DAG
1245: Not collective
1247: Input Parameters:
1248: + mesh - The DMPlex
1249: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1250: - size - The support size for point p
1252: Output Parameter:
1254: Note:
1255: This should be called after DMPlexSetChart().
1257: Level: beginner
1259: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1260: @*/
1261: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1262: {
1263: DM_Plex *mesh = (DM_Plex*) dm->data;
1268: PetscSectionSetDof(mesh->supportSection, p, size);
1270: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1271: return(0);
1272: }
1276: /*@C
1277: DMPlexGetSupport - Return the points on the out-edges for this point in the Sieve DAG
1279: Not collective
1281: Input Parameters:
1282: + mesh - The DMPlex
1283: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1285: Output Parameter:
1286: . support - An array of points which are on the out-edges for point p
1288: Level: beginner
1290: Fortran Notes:
1291: Since it returns an array, this routine is only available in Fortran 90, and you must
1292: include petsc.h90 in your code.
1294: You must also call DMPlexRestoreSupport() after you finish using the returned array.
1296: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1297: @*/
1298: PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1299: {
1300: DM_Plex *mesh = (DM_Plex*) dm->data;
1301: PetscInt off;
1307: PetscSectionGetOffset(mesh->supportSection, p, &off);
1308: *support = &mesh->supports[off];
1309: return(0);
1310: }
1314: /*@
1315: DMPlexSetSupport - Set the points on the out-edges for this point in the Sieve DAG
1317: Not collective
1319: Input Parameters:
1320: + mesh - The DMPlex
1321: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1322: - support - An array of points which are on the in-edges for point p
1324: Output Parameter:
1326: Note:
1327: This should be called after all calls to DMPlexSetSupportSize() and DMSetUp().
1329: Level: beginner
1331: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1332: @*/
1333: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1334: {
1335: DM_Plex *mesh = (DM_Plex*) dm->data;
1336: PetscInt pStart, pEnd;
1337: PetscInt dof, off, c;
1342: PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1343: PetscSectionGetDof(mesh->supportSection, p, &dof);
1345: PetscSectionGetOffset(mesh->supportSection, p, &off);
1346: 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);
1347: for (c = 0; c < dof; ++c) {
1348: 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);
1349: mesh->supports[off+c] = support[c];
1350: }
1351: return(0);
1352: }
1356: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1357: {
1358: DM_Plex *mesh = (DM_Plex*) dm->data;
1359: PetscInt pStart, pEnd;
1360: PetscInt dof, off;
1365: PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1366: PetscSectionGetDof(mesh->supportSection, p, &dof);
1367: PetscSectionGetOffset(mesh->supportSection, p, &off);
1368: 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);
1369: 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);
1370: 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);
1371: mesh->supports[off+supportPos] = supportPoint;
1372: return(0);
1373: }
1377: /*@C
1378: DMPlexGetTransitiveClosure - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1380: Not collective
1382: Input Parameters:
1383: + mesh - The DMPlex
1384: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1385: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1386: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
1388: Output Parameters:
1389: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1390: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
1392: Note:
1393: If using internal storage (points is NULL on input), each call overwrites the last output.
1395: Fortran Notes:
1396: Since it returns an array, this routine is only available in Fortran 90, and you must
1397: include petsc.h90 in your code.
1399: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1401: Level: beginner
1403: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1404: @*/
1405: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1406: {
1407: DM_Plex *mesh = (DM_Plex*) dm->data;
1408: PetscInt *closure, *fifo;
1409: const PetscInt *tmp = NULL, *tmpO = NULL;
1410: PetscInt tmpSize, t;
1411: PetscInt depth = 0, maxSize;
1412: PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0;
1413: PetscErrorCode ierr;
1417: DMPlexGetDepth(dm, &depth);
1418: /* This is only 1-level */
1419: if (useCone) {
1420: DMPlexGetConeSize(dm, p, &tmpSize);
1421: DMPlexGetCone(dm, p, &tmp);
1422: DMPlexGetConeOrientation(dm, p, &tmpO);
1423: } else {
1424: DMPlexGetSupportSize(dm, p, &tmpSize);
1425: DMPlexGetSupport(dm, p, &tmp);
1426: }
1427: if (depth == 1) {
1428: if (*points) {
1429: closure = *points;
1430: } else {
1431: maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1432: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1433: }
1434: closure[0] = p; closure[1] = 0;
1435: for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1436: closure[closureSize] = tmp[t];
1437: closure[closureSize+1] = tmpO ? tmpO[t] : 0;
1438: }
1439: if (numPoints) *numPoints = closureSize/2;
1440: if (points) *points = closure;
1441: return(0);
1442: }
1443: {
1444: PetscInt c, coneSeries, s,supportSeries;
1446: c = mesh->maxConeSize;
1447: coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1448: s = mesh->maxSupportSize;
1449: supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1450: maxSize = 2*PetscMax(coneSeries,supportSeries);
1451: }
1452: DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1453: if (*points) {
1454: closure = *points;
1455: } else {
1456: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1457: }
1458: closure[0] = p; closure[1] = 0;
1459: for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1460: const PetscInt cp = tmp[t];
1461: const PetscInt co = tmpO ? tmpO[t] : 0;
1463: closure[closureSize] = cp;
1464: closure[closureSize+1] = co;
1465: fifo[fifoSize] = cp;
1466: fifo[fifoSize+1] = co;
1467: }
1468: /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1469: while (fifoSize - fifoStart) {
1470: const PetscInt q = fifo[fifoStart];
1471: const PetscInt o = fifo[fifoStart+1];
1472: const PetscInt rev = o >= 0 ? 0 : 1;
1473: const PetscInt off = rev ? -(o+1) : o;
1475: if (useCone) {
1476: DMPlexGetConeSize(dm, q, &tmpSize);
1477: DMPlexGetCone(dm, q, &tmp);
1478: DMPlexGetConeOrientation(dm, q, &tmpO);
1479: } else {
1480: DMPlexGetSupportSize(dm, q, &tmpSize);
1481: DMPlexGetSupport(dm, q, &tmp);
1482: tmpO = NULL;
1483: }
1484: for (t = 0; t < tmpSize; ++t) {
1485: const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize;
1486: const PetscInt cp = tmp[i];
1487: /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1488: /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1489: const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1490: PetscInt co = tmpO ? tmpO[i] : 0;
1491: PetscInt c;
1493: if (rev) {
1494: PetscInt childSize, coff;
1495: DMPlexGetConeSize(dm, cp, &childSize);
1496: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1497: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1498: }
1499: /* Check for duplicate */
1500: for (c = 0; c < closureSize; c += 2) {
1501: if (closure[c] == cp) break;
1502: }
1503: if (c == closureSize) {
1504: closure[closureSize] = cp;
1505: closure[closureSize+1] = co;
1506: fifo[fifoSize] = cp;
1507: fifo[fifoSize+1] = co;
1508: closureSize += 2;
1509: fifoSize += 2;
1510: }
1511: }
1512: fifoStart += 2;
1513: }
1514: if (numPoints) *numPoints = closureSize/2;
1515: if (points) *points = closure;
1516: DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1517: return(0);
1518: }
1522: /*@C
1523: 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
1525: Not collective
1527: Input Parameters:
1528: + mesh - The DMPlex
1529: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1530: . orientation - The orientation of the point
1531: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1532: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
1534: Output Parameters:
1535: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1536: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
1538: Note:
1539: If using internal storage (points is NULL on input), each call overwrites the last output.
1541: Fortran Notes:
1542: Since it returns an array, this routine is only available in Fortran 90, and you must
1543: include petsc.h90 in your code.
1545: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1547: Level: beginner
1549: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1550: @*/
1551: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1552: {
1553: DM_Plex *mesh = (DM_Plex*) dm->data;
1554: PetscInt *closure, *fifo;
1555: const PetscInt *tmp = NULL, *tmpO = NULL;
1556: PetscInt tmpSize, t;
1557: PetscInt depth = 0, maxSize;
1558: PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0;
1559: PetscErrorCode ierr;
1563: DMPlexGetDepth(dm, &depth);
1564: /* This is only 1-level */
1565: if (useCone) {
1566: DMPlexGetConeSize(dm, p, &tmpSize);
1567: DMPlexGetCone(dm, p, &tmp);
1568: DMPlexGetConeOrientation(dm, p, &tmpO);
1569: } else {
1570: DMPlexGetSupportSize(dm, p, &tmpSize);
1571: DMPlexGetSupport(dm, p, &tmp);
1572: }
1573: if (depth == 1) {
1574: if (*points) {
1575: closure = *points;
1576: } else {
1577: maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1578: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1579: }
1580: closure[0] = p; closure[1] = ornt;
1581: for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1582: const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1583: closure[closureSize] = tmp[i];
1584: closure[closureSize+1] = tmpO ? tmpO[i] : 0;
1585: }
1586: if (numPoints) *numPoints = closureSize/2;
1587: if (points) *points = closure;
1588: return(0);
1589: }
1590: {
1591: PetscInt c, coneSeries, s,supportSeries;
1593: c = mesh->maxConeSize;
1594: coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1595: s = mesh->maxSupportSize;
1596: supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1597: maxSize = 2*PetscMax(coneSeries,supportSeries);
1598: }
1599: DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1600: if (*points) {
1601: closure = *points;
1602: } else {
1603: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1604: }
1605: closure[0] = p; closure[1] = ornt;
1606: for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1607: const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1608: const PetscInt cp = tmp[i];
1609: PetscInt co = tmpO ? tmpO[i] : 0;
1611: if (ornt < 0) {
1612: PetscInt childSize, coff;
1613: DMPlexGetConeSize(dm, cp, &childSize);
1614: coff = co < 0 ? -(tmpO[i]+1) : tmpO[i];
1615: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1616: }
1617: closure[closureSize] = cp;
1618: closure[closureSize+1] = co;
1619: fifo[fifoSize] = cp;
1620: fifo[fifoSize+1] = co;
1621: }
1622: /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1623: while (fifoSize - fifoStart) {
1624: const PetscInt q = fifo[fifoStart];
1625: const PetscInt o = fifo[fifoStart+1];
1626: const PetscInt rev = o >= 0 ? 0 : 1;
1627: const PetscInt off = rev ? -(o+1) : o;
1629: if (useCone) {
1630: DMPlexGetConeSize(dm, q, &tmpSize);
1631: DMPlexGetCone(dm, q, &tmp);
1632: DMPlexGetConeOrientation(dm, q, &tmpO);
1633: } else {
1634: DMPlexGetSupportSize(dm, q, &tmpSize);
1635: DMPlexGetSupport(dm, q, &tmp);
1636: tmpO = NULL;
1637: }
1638: for (t = 0; t < tmpSize; ++t) {
1639: const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize;
1640: const PetscInt cp = tmp[i];
1641: /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1642: /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1643: const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1644: PetscInt co = tmpO ? tmpO[i] : 0;
1645: PetscInt c;
1647: if (rev) {
1648: PetscInt childSize, coff;
1649: DMPlexGetConeSize(dm, cp, &childSize);
1650: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1651: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1652: }
1653: /* Check for duplicate */
1654: for (c = 0; c < closureSize; c += 2) {
1655: if (closure[c] == cp) break;
1656: }
1657: if (c == closureSize) {
1658: closure[closureSize] = cp;
1659: closure[closureSize+1] = co;
1660: fifo[fifoSize] = cp;
1661: fifo[fifoSize+1] = co;
1662: closureSize += 2;
1663: fifoSize += 2;
1664: }
1665: }
1666: fifoStart += 2;
1667: }
1668: if (numPoints) *numPoints = closureSize/2;
1669: if (points) *points = closure;
1670: DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1671: return(0);
1672: }
1676: /*@C
1677: DMPlexRestoreTransitiveClosure - Restore the array of points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1679: Not collective
1681: Input Parameters:
1682: + mesh - The DMPlex
1683: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1684: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1685: . numPoints - The number of points in the closure, so points[] is of size 2*numPoints, zeroed on exit
1686: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...], zeroed on exit
1688: Note:
1689: If not using internal storage (points is not NULL on input), this call is unnecessary
1691: Fortran Notes:
1692: Since it returns an array, this routine is only available in Fortran 90, and you must
1693: include petsc.h90 in your code.
1695: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1697: Level: beginner
1699: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1700: @*/
1701: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1702: {
1709: DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1710: if (numPoints) *numPoints = 0;
1711: return(0);
1712: }
1716: /*@
1717: DMPlexGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG
1719: Not collective
1721: Input Parameter:
1722: . mesh - The DMPlex
1724: Output Parameters:
1725: + maxConeSize - The maximum number of in-edges
1726: - maxSupportSize - The maximum number of out-edges
1728: Level: beginner
1730: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1731: @*/
1732: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1733: {
1734: DM_Plex *mesh = (DM_Plex*) dm->data;
1738: if (maxConeSize) *maxConeSize = mesh->maxConeSize;
1739: if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1740: return(0);
1741: }
1745: PetscErrorCode DMSetUp_Plex(DM dm)
1746: {
1747: DM_Plex *mesh = (DM_Plex*) dm->data;
1748: PetscInt size;
1753: PetscSectionSetUp(mesh->coneSection);
1754: PetscSectionGetStorageSize(mesh->coneSection, &size);
1755: PetscMalloc1(size, &mesh->cones);
1756: PetscCalloc1(size, &mesh->coneOrientations);
1757: if (mesh->maxSupportSize) {
1758: PetscSectionSetUp(mesh->supportSection);
1759: PetscSectionGetStorageSize(mesh->supportSection, &size);
1760: PetscMalloc1(size, &mesh->supports);
1761: }
1762: return(0);
1763: }
1767: PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1768: {
1772: if (subdm) {DMClone(dm, subdm);}
1773: DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1774: return(0);
1775: }
1779: /*@
1780: DMPlexSymmetrize - Creates support (out-edge) information from cone (in-edge) inoformation
1782: Not collective
1784: Input Parameter:
1785: . mesh - The DMPlex
1787: Output Parameter:
1789: Note:
1790: This should be called after all calls to DMPlexSetCone()
1792: Level: beginner
1794: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1795: @*/
1796: PetscErrorCode DMPlexSymmetrize(DM dm)
1797: {
1798: DM_Plex *mesh = (DM_Plex*) dm->data;
1799: PetscInt *offsets;
1800: PetscInt supportSize;
1801: PetscInt pStart, pEnd, p;
1806: if (mesh->supports) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Supports were already setup in this DMPlex");
1807: /* Calculate support sizes */
1808: DMPlexGetChart(dm, &pStart, &pEnd);
1809: for (p = pStart; p < pEnd; ++p) {
1810: PetscInt dof, off, c;
1812: PetscSectionGetDof(mesh->coneSection, p, &dof);
1813: PetscSectionGetOffset(mesh->coneSection, p, &off);
1814: for (c = off; c < off+dof; ++c) {
1815: PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1816: }
1817: }
1818: for (p = pStart; p < pEnd; ++p) {
1819: PetscInt dof;
1821: PetscSectionGetDof(mesh->supportSection, p, &dof);
1823: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1824: }
1825: PetscSectionSetUp(mesh->supportSection);
1826: /* Calculate supports */
1827: PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1828: PetscMalloc1(supportSize, &mesh->supports);
1829: PetscCalloc1(pEnd - pStart, &offsets);
1830: for (p = pStart; p < pEnd; ++p) {
1831: PetscInt dof, off, c;
1833: PetscSectionGetDof(mesh->coneSection, p, &dof);
1834: PetscSectionGetOffset(mesh->coneSection, p, &off);
1835: for (c = off; c < off+dof; ++c) {
1836: const PetscInt q = mesh->cones[c];
1837: PetscInt offS;
1839: PetscSectionGetOffset(mesh->supportSection, q, &offS);
1841: mesh->supports[offS+offsets[q]] = p;
1842: ++offsets[q];
1843: }
1844: }
1845: PetscFree(offsets);
1846: return(0);
1847: }
1851: /*@
1852: DMPlexStratify - The Sieve DAG for most topologies is a graded poset (http://en.wikipedia.org/wiki/Graded_poset), and
1853: can be illustrated by Hasse Diagram (a http://en.wikipedia.org/wiki/Hasse_diagram). The strata group all points of the
1854: same grade, and this function calculates the strata. This grade can be seen as the height (or depth) of the point in
1855: the DAG.
1857: Collective on dm
1859: Input Parameter:
1860: . mesh - The DMPlex
1862: Output Parameter:
1864: Notes:
1865: Concretely, DMPlexStratify() creates a new label named "depth" containing the dimension of each element: 0 for vertices,
1866: 1 for edges, and so on. The depth label can be accessed through DMPlexGetDepthLabel() or DMPlexGetDepthStratum(), or
1867: manually via DMGetLabel(). The height is defined implicitly by height = maxDimension - depth, and can be accessed
1868: via DMPlexGetHeightStratum(). For example, cells have height 0 and faces have height 1.
1870: DMPlexStratify() should be called after all calls to DMPlexSymmetrize()
1872: Level: beginner
1874: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1875: @*/
1876: PetscErrorCode DMPlexStratify(DM dm)
1877: {
1878: DM_Plex *mesh = (DM_Plex*) dm->data;
1879: DMLabel label;
1880: PetscInt pStart, pEnd, p;
1881: PetscInt numRoots = 0, numLeaves = 0;
1886: PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1887: /* Calculate depth */
1888: DMPlexGetChart(dm, &pStart, &pEnd);
1889: DMCreateLabel(dm, "depth");
1890: DMPlexGetDepthLabel(dm, &label);
1891: /* Initialize roots and count leaves */
1892: for (p = pStart; p < pEnd; ++p) {
1893: PetscInt coneSize, supportSize;
1895: DMPlexGetConeSize(dm, p, &coneSize);
1896: DMPlexGetSupportSize(dm, p, &supportSize);
1897: if (!coneSize && supportSize) {
1898: ++numRoots;
1899: DMLabelSetValue(label, p, 0);
1900: } else if (!supportSize && coneSize) {
1901: ++numLeaves;
1902: } else if (!supportSize && !coneSize) {
1903: /* Isolated points */
1904: DMLabelSetValue(label, p, 0);
1905: }
1906: }
1907: if (numRoots + numLeaves == (pEnd - pStart)) {
1908: for (p = pStart; p < pEnd; ++p) {
1909: PetscInt coneSize, supportSize;
1911: DMPlexGetConeSize(dm, p, &coneSize);
1912: DMPlexGetSupportSize(dm, p, &supportSize);
1913: if (!supportSize && coneSize) {
1914: DMLabelSetValue(label, p, 1);
1915: }
1916: }
1917: } else {
1918: IS pointIS;
1919: PetscInt numPoints = 0, level = 0;
1921: DMLabelGetStratumIS(label, level, &pointIS);
1922: if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1923: while (numPoints) {
1924: const PetscInt *points;
1925: const PetscInt newLevel = level+1;
1927: ISGetIndices(pointIS, &points);
1928: for (p = 0; p < numPoints; ++p) {
1929: const PetscInt point = points[p];
1930: const PetscInt *support;
1931: PetscInt supportSize, s;
1933: DMPlexGetSupportSize(dm, point, &supportSize);
1934: DMPlexGetSupport(dm, point, &support);
1935: for (s = 0; s < supportSize; ++s) {
1936: DMLabelSetValue(label, support[s], newLevel);
1937: }
1938: }
1939: ISRestoreIndices(pointIS, &points);
1940: ++level;
1941: ISDestroy(&pointIS);
1942: DMLabelGetStratumIS(label, level, &pointIS);
1943: if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1944: else {numPoints = 0;}
1945: }
1946: ISDestroy(&pointIS);
1947: }
1948: { /* just in case there is an empty process */
1949: PetscInt numValues, maxValues = 0, v;
1951: DMLabelGetNumValues(label,&numValues);
1952: MPI_Allreduce(&numValues,&maxValues,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
1953: for (v = numValues; v < maxValues; v++) {
1954: DMLabelAddStratum(label,v);
1955: }
1956: }
1958: DMLabelGetState(label, &mesh->depthState);
1959: PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1960: return(0);
1961: }
1965: /*@C
1966: DMPlexGetJoin - Get an array for the join of the set of points
1968: Not Collective
1970: Input Parameters:
1971: + dm - The DMPlex object
1972: . numPoints - The number of input points for the join
1973: - points - The input points
1975: Output Parameters:
1976: + numCoveredPoints - The number of points in the join
1977: - coveredPoints - The points in the join
1979: Level: intermediate
1981: Note: Currently, this is restricted to a single level join
1983: Fortran Notes:
1984: Since it returns an array, this routine is only available in Fortran 90, and you must
1985: include petsc.h90 in your code.
1987: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1989: .keywords: mesh
1990: .seealso: DMPlexRestoreJoin(), DMPlexGetMeet()
1991: @*/
1992: PetscErrorCode DMPlexGetJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1993: {
1994: DM_Plex *mesh = (DM_Plex*) dm->data;
1995: PetscInt *join[2];
1996: PetscInt joinSize, i = 0;
1997: PetscInt dof, off, p, c, m;
2005: DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[0]);
2006: DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1]);
2007: /* Copy in support of first point */
2008: PetscSectionGetDof(mesh->supportSection, points[0], &dof);
2009: PetscSectionGetOffset(mesh->supportSection, points[0], &off);
2010: for (joinSize = 0; joinSize < dof; ++joinSize) {
2011: join[i][joinSize] = mesh->supports[off+joinSize];
2012: }
2013: /* Check each successive support */
2014: for (p = 1; p < numPoints; ++p) {
2015: PetscInt newJoinSize = 0;
2017: PetscSectionGetDof(mesh->supportSection, points[p], &dof);
2018: PetscSectionGetOffset(mesh->supportSection, points[p], &off);
2019: for (c = 0; c < dof; ++c) {
2020: const PetscInt point = mesh->supports[off+c];
2022: for (m = 0; m < joinSize; ++m) {
2023: if (point == join[i][m]) {
2024: join[1-i][newJoinSize++] = point;
2025: break;
2026: }
2027: }
2028: }
2029: joinSize = newJoinSize;
2030: i = 1-i;
2031: }
2032: *numCoveredPoints = joinSize;
2033: *coveredPoints = join[i];
2034: DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2035: return(0);
2036: }
2040: /*@C
2041: DMPlexRestoreJoin - Restore an array for the join of the set of points
2043: Not Collective
2045: Input Parameters:
2046: + dm - The DMPlex object
2047: . numPoints - The number of input points for the join
2048: - points - The input points
2050: Output Parameters:
2051: + numCoveredPoints - The number of points in the join
2052: - coveredPoints - The points in the join
2054: Fortran Notes:
2055: Since it returns an array, this routine is only available in Fortran 90, and you must
2056: include petsc.h90 in your code.
2058: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2060: Level: intermediate
2062: .keywords: mesh
2063: .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet()
2064: @*/
2065: PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2066: {
2074: DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2075: if (numCoveredPoints) *numCoveredPoints = 0;
2076: return(0);
2077: }
2081: /*@C
2082: DMPlexGetFullJoin - Get an array for the join of the set of points
2084: Not Collective
2086: Input Parameters:
2087: + dm - The DMPlex object
2088: . numPoints - The number of input points for the join
2089: - points - The input points
2091: Output Parameters:
2092: + numCoveredPoints - The number of points in the join
2093: - coveredPoints - The points in the join
2095: Fortran Notes:
2096: Since it returns an array, this routine is only available in Fortran 90, and you must
2097: include petsc.h90 in your code.
2099: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2101: Level: intermediate
2103: .keywords: mesh
2104: .seealso: DMPlexGetJoin(), DMPlexRestoreJoin(), DMPlexGetMeet()
2105: @*/
2106: PetscErrorCode DMPlexGetFullJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2107: {
2108: DM_Plex *mesh = (DM_Plex*) dm->data;
2109: PetscInt *offsets, **closures;
2110: PetscInt *join[2];
2111: PetscInt depth = 0, maxSize, joinSize = 0, i = 0;
2112: PetscInt p, d, c, m, ms;
2121: DMPlexGetDepth(dm, &depth);
2122: PetscCalloc1(numPoints, &closures);
2123: DMGetWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2124: ms = mesh->maxSupportSize;
2125: maxSize = (ms > 1) ? ((PetscPowInt(ms,depth+1)-1)/(ms-1)) : depth + 1;
2126: DMGetWorkArray(dm, maxSize, PETSC_INT, &join[0]);
2127: DMGetWorkArray(dm, maxSize, PETSC_INT, &join[1]);
2129: for (p = 0; p < numPoints; ++p) {
2130: PetscInt closureSize;
2132: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closures[p]);
2134: offsets[p*(depth+2)+0] = 0;
2135: for (d = 0; d < depth+1; ++d) {
2136: PetscInt pStart, pEnd, i;
2138: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
2139: for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) {
2140: if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2141: offsets[p*(depth+2)+d+1] = i;
2142: break;
2143: }
2144: }
2145: if (i == closureSize) offsets[p*(depth+2)+d+1] = i;
2146: }
2147: 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);
2148: }
2149: for (d = 0; d < depth+1; ++d) {
2150: PetscInt dof;
2152: /* Copy in support of first point */
2153: dof = offsets[d+1] - offsets[d];
2154: for (joinSize = 0; joinSize < dof; ++joinSize) {
2155: join[i][joinSize] = closures[0][(offsets[d]+joinSize)*2];
2156: }
2157: /* Check each successive cone */
2158: for (p = 1; p < numPoints && joinSize; ++p) {
2159: PetscInt newJoinSize = 0;
2161: dof = offsets[p*(depth+2)+d+1] - offsets[p*(depth+2)+d];
2162: for (c = 0; c < dof; ++c) {
2163: const PetscInt point = closures[p][(offsets[p*(depth+2)+d]+c)*2];
2165: for (m = 0; m < joinSize; ++m) {
2166: if (point == join[i][m]) {
2167: join[1-i][newJoinSize++] = point;
2168: break;
2169: }
2170: }
2171: }
2172: joinSize = newJoinSize;
2173: i = 1-i;
2174: }
2175: if (joinSize) break;
2176: }
2177: *numCoveredPoints = joinSize;
2178: *coveredPoints = join[i];
2179: for (p = 0; p < numPoints; ++p) {
2180: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, NULL, &closures[p]);
2181: }
2182: PetscFree(closures);
2183: DMRestoreWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2184: DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2185: return(0);
2186: }
2190: /*@C
2191: DMPlexGetMeet - Get an array for the meet of the set of points
2193: Not Collective
2195: Input Parameters:
2196: + dm - The DMPlex object
2197: . numPoints - The number of input points for the meet
2198: - points - The input points
2200: Output Parameters:
2201: + numCoveredPoints - The number of points in the meet
2202: - coveredPoints - The points in the meet
2204: Level: intermediate
2206: Note: Currently, this is restricted to a single level meet
2208: Fortran Notes:
2209: Since it returns an array, this routine is only available in Fortran 90, and you must
2210: include petsc.h90 in your code.
2212: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2214: .keywords: mesh
2215: .seealso: DMPlexRestoreMeet(), DMPlexGetJoin()
2216: @*/
2217: PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2218: {
2219: DM_Plex *mesh = (DM_Plex*) dm->data;
2220: PetscInt *meet[2];
2221: PetscInt meetSize, i = 0;
2222: PetscInt dof, off, p, c, m;
2230: DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[0]);
2231: DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1]);
2232: /* Copy in cone of first point */
2233: PetscSectionGetDof(mesh->coneSection, points[0], &dof);
2234: PetscSectionGetOffset(mesh->coneSection, points[0], &off);
2235: for (meetSize = 0; meetSize < dof; ++meetSize) {
2236: meet[i][meetSize] = mesh->cones[off+meetSize];
2237: }
2238: /* Check each successive cone */
2239: for (p = 1; p < numPoints; ++p) {
2240: PetscInt newMeetSize = 0;
2242: PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2243: PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2244: for (c = 0; c < dof; ++c) {
2245: const PetscInt point = mesh->cones[off+c];
2247: for (m = 0; m < meetSize; ++m) {
2248: if (point == meet[i][m]) {
2249: meet[1-i][newMeetSize++] = point;
2250: break;
2251: }
2252: }
2253: }
2254: meetSize = newMeetSize;
2255: i = 1-i;
2256: }
2257: *numCoveringPoints = meetSize;
2258: *coveringPoints = meet[i];
2259: DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2260: return(0);
2261: }
2265: /*@C
2266: DMPlexRestoreMeet - Restore an array for the meet of the set of points
2268: Not Collective
2270: Input Parameters:
2271: + dm - The DMPlex object
2272: . numPoints - The number of input points for the meet
2273: - points - The input points
2275: Output Parameters:
2276: + numCoveredPoints - The number of points in the meet
2277: - coveredPoints - The points in the meet
2279: Level: intermediate
2281: Fortran Notes:
2282: Since it returns an array, this routine is only available in Fortran 90, and you must
2283: include petsc.h90 in your code.
2285: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2287: .keywords: mesh
2288: .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin()
2289: @*/
2290: PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2291: {
2299: DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2300: if (numCoveredPoints) *numCoveredPoints = 0;
2301: return(0);
2302: }
2306: /*@C
2307: DMPlexGetFullMeet - Get an array for the meet of the set of points
2309: Not Collective
2311: Input Parameters:
2312: + dm - The DMPlex object
2313: . numPoints - The number of input points for the meet
2314: - points - The input points
2316: Output Parameters:
2317: + numCoveredPoints - The number of points in the meet
2318: - coveredPoints - The points in the meet
2320: Level: intermediate
2322: Fortran Notes:
2323: Since it returns an array, this routine is only available in Fortran 90, and you must
2324: include petsc.h90 in your code.
2326: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2328: .keywords: mesh
2329: .seealso: DMPlexGetMeet(), DMPlexRestoreMeet(), DMPlexGetJoin()
2330: @*/
2331: PetscErrorCode DMPlexGetFullMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2332: {
2333: DM_Plex *mesh = (DM_Plex*) dm->data;
2334: PetscInt *offsets, **closures;
2335: PetscInt *meet[2];
2336: PetscInt height = 0, maxSize, meetSize = 0, i = 0;
2337: PetscInt p, h, c, m, mc;
2346: DMPlexGetDepth(dm, &height);
2347: PetscMalloc1(numPoints, &closures);
2348: DMGetWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2349: mc = mesh->maxConeSize;
2350: maxSize = (mc > 1) ? ((PetscPowInt(mc,height+1)-1)/(mc-1)) : height + 1;
2351: DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[0]);
2352: DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[1]);
2354: for (p = 0; p < numPoints; ++p) {
2355: PetscInt closureSize;
2357: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closures[p]);
2359: offsets[p*(height+2)+0] = 0;
2360: for (h = 0; h < height+1; ++h) {
2361: PetscInt pStart, pEnd, i;
2363: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
2364: for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) {
2365: if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2366: offsets[p*(height+2)+h+1] = i;
2367: break;
2368: }
2369: }
2370: if (i == closureSize) offsets[p*(height+2)+h+1] = i;
2371: }
2372: 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);
2373: }
2374: for (h = 0; h < height+1; ++h) {
2375: PetscInt dof;
2377: /* Copy in cone of first point */
2378: dof = offsets[h+1] - offsets[h];
2379: for (meetSize = 0; meetSize < dof; ++meetSize) {
2380: meet[i][meetSize] = closures[0][(offsets[h]+meetSize)*2];
2381: }
2382: /* Check each successive cone */
2383: for (p = 1; p < numPoints && meetSize; ++p) {
2384: PetscInt newMeetSize = 0;
2386: dof = offsets[p*(height+2)+h+1] - offsets[p*(height+2)+h];
2387: for (c = 0; c < dof; ++c) {
2388: const PetscInt point = closures[p][(offsets[p*(height+2)+h]+c)*2];
2390: for (m = 0; m < meetSize; ++m) {
2391: if (point == meet[i][m]) {
2392: meet[1-i][newMeetSize++] = point;
2393: break;
2394: }
2395: }
2396: }
2397: meetSize = newMeetSize;
2398: i = 1-i;
2399: }
2400: if (meetSize) break;
2401: }
2402: *numCoveredPoints = meetSize;
2403: *coveredPoints = meet[i];
2404: for (p = 0; p < numPoints; ++p) {
2405: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, NULL, &closures[p]);
2406: }
2407: PetscFree(closures);
2408: DMRestoreWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2409: DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2410: return(0);
2411: }
2415: /*@C
2416: DMPlexEqual - Determine if two DMs have the same topology
2418: Not Collective
2420: Input Parameters:
2421: + dmA - A DMPlex object
2422: - dmB - A DMPlex object
2424: Output Parameters:
2425: . equal - PETSC_TRUE if the topologies are identical
2427: Level: intermediate
2429: Notes:
2430: We are not solving graph isomorphism, so we do not permutation.
2432: .keywords: mesh
2433: .seealso: DMPlexGetCone()
2434: @*/
2435: PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal)
2436: {
2437: PetscInt depth, depthB, pStart, pEnd, pStartB, pEndB, p;
2445: *equal = PETSC_FALSE;
2446: DMPlexGetDepth(dmA, &depth);
2447: DMPlexGetDepth(dmB, &depthB);
2448: if (depth != depthB) return(0);
2449: DMPlexGetChart(dmA, &pStart, &pEnd);
2450: DMPlexGetChart(dmB, &pStartB, &pEndB);
2451: if ((pStart != pStartB) || (pEnd != pEndB)) return(0);
2452: for (p = pStart; p < pEnd; ++p) {
2453: const PetscInt *cone, *coneB, *ornt, *orntB, *support, *supportB;
2454: PetscInt coneSize, coneSizeB, c, supportSize, supportSizeB, s;
2456: DMPlexGetConeSize(dmA, p, &coneSize);
2457: DMPlexGetCone(dmA, p, &cone);
2458: DMPlexGetConeOrientation(dmA, p, &ornt);
2459: DMPlexGetConeSize(dmB, p, &coneSizeB);
2460: DMPlexGetCone(dmB, p, &coneB);
2461: DMPlexGetConeOrientation(dmB, p, &orntB);
2462: if (coneSize != coneSizeB) return(0);
2463: for (c = 0; c < coneSize; ++c) {
2464: if (cone[c] != coneB[c]) return(0);
2465: if (ornt[c] != orntB[c]) return(0);
2466: }
2467: DMPlexGetSupportSize(dmA, p, &supportSize);
2468: DMPlexGetSupport(dmA, p, &support);
2469: DMPlexGetSupportSize(dmB, p, &supportSizeB);
2470: DMPlexGetSupport(dmB, p, &supportB);
2471: if (supportSize != supportSizeB) return(0);
2472: for (s = 0; s < supportSize; ++s) {
2473: if (support[s] != supportB[s]) return(0);
2474: }
2475: }
2476: *equal = PETSC_TRUE;
2477: return(0);
2478: }
2482: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2483: {
2484: MPI_Comm comm;
2488: PetscObjectGetComm((PetscObject)dm,&comm);
2490: switch (cellDim) {
2491: case 0:
2492: *numFaceVertices = 0;
2493: break;
2494: case 1:
2495: *numFaceVertices = 1;
2496: break;
2497: case 2:
2498: switch (numCorners) {
2499: case 3: /* triangle */
2500: *numFaceVertices = 2; /* Edge has 2 vertices */
2501: break;
2502: case 4: /* quadrilateral */
2503: *numFaceVertices = 2; /* Edge has 2 vertices */
2504: break;
2505: case 6: /* quadratic triangle, tri and quad cohesive Lagrange cells */
2506: *numFaceVertices = 3; /* Edge has 3 vertices */
2507: break;
2508: case 9: /* quadratic quadrilateral, quadratic quad cohesive Lagrange cells */
2509: *numFaceVertices = 3; /* Edge has 3 vertices */
2510: break;
2511: default:
2512: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %D for dimension %D", numCorners, cellDim);
2513: }
2514: break;
2515: case 3:
2516: switch (numCorners) {
2517: case 4: /* tetradehdron */
2518: *numFaceVertices = 3; /* Face has 3 vertices */
2519: break;
2520: case 6: /* tet cohesive cells */
2521: *numFaceVertices = 4; /* Face has 4 vertices */
2522: break;
2523: case 8: /* hexahedron */
2524: *numFaceVertices = 4; /* Face has 4 vertices */
2525: break;
2526: case 9: /* tet cohesive Lagrange cells */
2527: *numFaceVertices = 6; /* Face has 6 vertices */
2528: break;
2529: case 10: /* quadratic tetrahedron */
2530: *numFaceVertices = 6; /* Face has 6 vertices */
2531: break;
2532: case 12: /* hex cohesive Lagrange cells */
2533: *numFaceVertices = 6; /* Face has 6 vertices */
2534: break;
2535: case 18: /* quadratic tet cohesive Lagrange cells */
2536: *numFaceVertices = 6; /* Face has 6 vertices */
2537: break;
2538: case 27: /* quadratic hexahedron, quadratic hex cohesive Lagrange cells */
2539: *numFaceVertices = 9; /* Face has 9 vertices */
2540: break;
2541: default:
2542: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %D for dimension %D", numCorners, cellDim);
2543: }
2544: break;
2545: default:
2546: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid cell dimension %D", cellDim);
2547: }
2548: return(0);
2549: }
2553: /*@
2554: DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point
2556: Not Collective
2558: Input Parameter:
2559: . dm - The DMPlex object
2561: Output Parameter:
2562: . depthLabel - The DMLabel recording point depth
2564: Level: developer
2566: .keywords: mesh, points
2567: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2568: @*/
2569: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
2570: {
2576: if (!dm->depthLabel) {DMGetLabel(dm, "depth", &dm->depthLabel);}
2577: *depthLabel = dm->depthLabel;
2578: return(0);
2579: }
2583: /*@
2584: DMPlexGetDepth - Get the depth of the DAG representing this mesh
2586: Not Collective
2588: Input Parameter:
2589: . dm - The DMPlex object
2591: Output Parameter:
2592: . depth - The number of strata (breadth first levels) in the DAG
2594: Level: developer
2596: .keywords: mesh, points
2597: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2598: @*/
2599: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
2600: {
2601: DMLabel label;
2602: PetscInt d = 0;
2608: DMPlexGetDepthLabel(dm, &label);
2609: if (label) {DMLabelGetNumValues(label, &d);}
2610: *depth = d-1;
2611: return(0);
2612: }
2616: /*@
2617: DMPlexGetDepthStratum - Get the bounds [start, end) for all points at a certain depth.
2619: Not Collective
2621: Input Parameters:
2622: + dm - The DMPlex object
2623: - stratumValue - The requested depth
2625: Output Parameters:
2626: + start - The first point at this depth
2627: - end - One beyond the last point at this depth
2629: Level: developer
2631: .keywords: mesh, points
2632: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
2633: @*/
2634: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2635: {
2636: DMLabel label;
2637: PetscInt pStart, pEnd;
2644: DMPlexGetChart(dm, &pStart, &pEnd);
2645: if (pStart == pEnd) return(0);
2646: if (stratumValue < 0) {
2647: if (start) *start = pStart;
2648: if (end) *end = pEnd;
2649: return(0);
2650: }
2651: DMPlexGetDepthLabel(dm, &label);
2652: if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2653: DMLabelGetStratumBounds(label, stratumValue, start, end);
2654: return(0);
2655: }
2659: /*@
2660: DMPlexGetHeightStratum - Get the bounds [start, end) for all points at a certain height.
2662: Not Collective
2664: Input Parameters:
2665: + dm - The DMPlex object
2666: - stratumValue - The requested height
2668: Output Parameters:
2669: + start - The first point at this height
2670: - end - One beyond the last point at this height
2672: Level: developer
2674: .keywords: mesh, points
2675: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
2676: @*/
2677: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2678: {
2679: DMLabel label;
2680: PetscInt depth, pStart, pEnd;
2687: DMPlexGetChart(dm, &pStart, &pEnd);
2688: if (pStart == pEnd) return(0);
2689: if (stratumValue < 0) {
2690: if (start) *start = pStart;
2691: if (end) *end = pEnd;
2692: return(0);
2693: }
2694: DMPlexGetDepthLabel(dm, &label);
2695: if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2696: DMLabelGetNumValues(label, &depth);
2697: DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
2698: return(0);
2699: }
2703: /* Set the number of dof on each point and separate by fields */
2704: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
2705: {
2706: PetscInt *pMax;
2707: PetscInt depth, pStart = 0, pEnd = 0;
2708: PetscInt Nf, p, d, dep, f;
2709: PetscBool *isFE;
2713: PetscMalloc1(numFields, &isFE);
2714: DMGetNumFields(dm, &Nf);
2715: for (f = 0; f < numFields; ++f) {
2716: PetscObject obj;
2717: PetscClassId id;
2719: isFE[f] = PETSC_FALSE;
2720: if (f >= Nf) continue;
2721: DMGetField(dm, f, &obj);
2722: PetscObjectGetClassId(obj, &id);
2723: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
2724: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
2725: }
2726: PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
2727: if (numFields > 0) {
2728: PetscSectionSetNumFields(*section, numFields);
2729: if (numComp) {
2730: for (f = 0; f < numFields; ++f) {
2731: PetscSectionSetFieldComponents(*section, f, numComp[f]);
2732: }
2733: }
2734: }
2735: DMPlexGetChart(dm, &pStart, &pEnd);
2736: PetscSectionSetChart(*section, pStart, pEnd);
2737: DMPlexGetDepth(dm, &depth);
2738: PetscMalloc1(depth+1,&pMax);
2739: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2740: for (dep = 0; dep <= depth; ++dep) {
2741: d = dim == depth ? dep : (!dep ? 0 : dim);
2742: DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
2743: pMax[dep] = pMax[dep] < 0 ? pEnd : pMax[dep];
2744: for (p = pStart; p < pEnd; ++p) {
2745: PetscInt tot = 0;
2747: for (f = 0; f < numFields; ++f) {
2748: if (isFE[f] && p >= pMax[dep]) continue;
2749: PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
2750: tot += numDof[f*(dim+1)+d];
2751: }
2752: PetscSectionSetDof(*section, p, tot);
2753: }
2754: }
2755: PetscFree(pMax);
2756: PetscFree(isFE);
2757: return(0);
2758: }
2762: /* Set the number of dof on each point and separate by fields
2763: If bcComps is NULL or the IS is NULL, constrain every dof on the point
2764: */
2765: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2766: {
2767: PetscInt numFields;
2768: PetscInt bc;
2769: PetscSection aSec;
2773: PetscSectionGetNumFields(section, &numFields);
2774: for (bc = 0; bc < numBC; ++bc) {
2775: PetscInt field = 0;
2776: const PetscInt *comp;
2777: const PetscInt *idx;
2778: PetscInt Nc = -1, n, i;
2780: if (numFields) field = bcField[bc];
2781: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2782: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2783: ISGetLocalSize(bcPoints[bc], &n);
2784: ISGetIndices(bcPoints[bc], &idx);
2785: for (i = 0; i < n; ++i) {
2786: const PetscInt p = idx[i];
2787: PetscInt numConst;
2789: if (numFields) {
2790: PetscSectionGetFieldDof(section, p, field, &numConst);
2791: } else {
2792: PetscSectionGetDof(section, p, &numConst);
2793: }
2794: /* If Nc < 0, constrain every dof on the point */
2795: if (Nc > 0) numConst = PetscMin(numConst, Nc);
2796: if (numFields) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
2797: PetscSectionAddConstraintDof(section, p, numConst);
2798: }
2799: ISRestoreIndices(bcPoints[bc], &idx);
2800: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2801: }
2802: DMPlexGetAnchors(dm, &aSec, NULL);
2803: if (aSec) {
2804: PetscInt aStart, aEnd, a;
2806: PetscSectionGetChart(aSec, &aStart, &aEnd);
2807: for (a = aStart; a < aEnd; a++) {
2808: PetscInt dof, f;
2810: PetscSectionGetDof(aSec, a, &dof);
2811: if (dof) {
2812: /* if there are point-to-point constraints, then all dofs are constrained */
2813: PetscSectionGetDof(section, a, &dof);
2814: PetscSectionSetConstraintDof(section, a, dof);
2815: for (f = 0; f < numFields; f++) {
2816: PetscSectionGetFieldDof(section, a, f, &dof);
2817: PetscSectionSetFieldConstraintDof(section, a, f, dof);
2818: }
2819: }
2820: }
2821: }
2822: return(0);
2823: }
2827: /* Set the constrained field indices on each point
2828: If bcComps is NULL or the IS is NULL, constrain every dof on the point
2829: */
2830: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2831: {
2832: PetscSection aSec;
2833: PetscInt *indices;
2834: PetscInt numFields, maxDof, pStart, pEnd, p, bc, f, d;
2838: PetscSectionGetNumFields(section, &numFields);
2839: if (!numFields) return(0);
2840: /* Initialize all field indices to -1 */
2841: PetscSectionGetChart(section, &pStart, &pEnd);
2842: PetscSectionGetMaxDof(section, &maxDof);
2843: PetscMalloc1(maxDof, &indices);
2844: for (d = 0; d < maxDof; ++d) indices[d] = -1;
2845: for (p = pStart; p < pEnd; ++p) for (f = 0; f < numFields; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
2846: /* Handle BC constraints */
2847: for (bc = 0; bc < numBC; ++bc) {
2848: const PetscInt field = bcField[bc];
2849: const PetscInt *comp, *idx;
2850: PetscInt Nc = -1, n, i;
2852: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2853: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2854: ISGetLocalSize(bcPoints[bc], &n);
2855: ISGetIndices(bcPoints[bc], &idx);
2856: for (i = 0; i < n; ++i) {
2857: const PetscInt p = idx[i];
2858: const PetscInt *find;
2859: PetscInt fcdof, c;
2861: PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
2862: if (Nc < 0) {
2863: for (d = 0; d < fcdof; ++d) indices[d] = d;
2864: } else {
2865: PetscSectionGetFieldConstraintIndices(section, p, field, &find);
2866: for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
2867: for (c = 0; c < Nc; ++c) indices[d+c] = comp[c];
2868: PetscSortInt(d+Nc, indices);
2869: for (c = d+Nc; c < fcdof; ++c) indices[c] = -1;
2870: }
2871: PetscSectionSetFieldConstraintIndices(section, p, field, indices);
2872: }
2873: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2874: ISRestoreIndices(bcPoints[bc], &idx);
2875: }
2876: /* Handle anchors */
2877: DMPlexGetAnchors(dm, &aSec, NULL);
2878: if (aSec) {
2879: PetscInt aStart, aEnd, a;
2881: for (d = 0; d < maxDof; ++d) indices[d] = d;
2882: PetscSectionGetChart(aSec, &aStart, &aEnd);
2883: for (a = aStart; a < aEnd; a++) {
2884: PetscInt dof, fdof, f;
2886: PetscSectionGetDof(aSec, a, &dof);
2887: if (dof) {
2888: /* if there are point-to-point constraints, then all dofs are constrained */
2889: for (f = 0; f < numFields; f++) {
2890: PetscSectionGetFieldDof(section, a, f, &fdof);
2891: PetscSectionSetFieldConstraintIndices(section, a, f, indices);
2892: }
2893: }
2894: }
2895: }
2896: PetscFree(indices);
2897: return(0);
2898: }
2902: /* Set the constrained indices on each point */
2903: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
2904: {
2905: PetscInt *indices;
2906: PetscInt numFields, maxDof, pStart, pEnd, p, f, d;
2910: PetscSectionGetNumFields(section, &numFields);
2911: PetscSectionGetMaxDof(section, &maxDof);
2912: PetscSectionGetChart(section, &pStart, &pEnd);
2913: PetscMalloc1(maxDof, &indices);
2914: for (d = 0; d < maxDof; ++d) indices[d] = -1;
2915: for (p = pStart; p < pEnd; ++p) {
2916: PetscInt cdof, d;
2918: PetscSectionGetConstraintDof(section, p, &cdof);
2919: if (cdof) {
2920: if (numFields) {
2921: PetscInt numConst = 0, foff = 0;
2923: for (f = 0; f < numFields; ++f) {
2924: const PetscInt *find;
2925: PetscInt fcdof, fdof;
2927: PetscSectionGetFieldDof(section, p, f, &fdof);
2928: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
2929: /* Change constraint numbering from field component to local dof number */
2930: PetscSectionGetFieldConstraintIndices(section, p, f, &find);
2931: for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
2932: numConst += fcdof;
2933: foff += fdof;
2934: }
2935: if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
2936: } else {
2937: for (d = 0; d < cdof; ++d) indices[d] = d;
2938: }
2939: PetscSectionSetConstraintIndices(section, p, indices);
2940: }
2941: }
2942: PetscFree(indices);
2943: return(0);
2944: }
2948: /*@C
2949: DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided.
2951: Not Collective
2953: Input Parameters:
2954: + dm - The DMPlex object
2955: . dim - The spatial dimension of the problem
2956: . numFields - The number of fields in the problem
2957: . numComp - An array of size numFields that holds the number of components for each field
2958: . numDof - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
2959: . numBC - The number of boundary conditions
2960: . bcField - An array of size numBC giving the field number for each boundry condition
2961: . bcComps - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
2962: . bcPoints - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
2963: - perm - Optional permutation of the chart, or NULL
2965: Output Parameter:
2966: . section - The PetscSection object
2968: 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
2969: number of dof for field 0 on each edge.
2971: The chart permutation is the same one set using PetscSectionSetPermutation()
2973: Level: developer
2975: Fortran Notes:
2976: A Fortran 90 version is available as DMPlexCreateSectionF90()
2978: .keywords: mesh, elements
2979: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
2980: @*/
2981: 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)
2982: {
2983: PetscSection aSec;
2987: DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
2988: DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
2989: if (perm) {PetscSectionSetPermutation(*section, perm);}
2990: PetscSectionSetUp(*section);
2991: DMPlexGetAnchors(dm,&aSec,NULL);
2992: if (numBC || aSec) {
2993: DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
2994: DMPlexCreateSectionBCIndices(dm, *section);
2995: }
2996: PetscSectionViewFromOptions(*section,NULL,"-section_view");
2997: return(0);
2998: }
3002: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
3003: {
3004: PetscSection section, s;
3005: Mat m;
3009: DMClone(dm, cdm);
3010: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
3011: DMSetDefaultSection(*cdm, section);
3012: PetscSectionDestroy(§ion);
3013: PetscSectionCreate(PETSC_COMM_SELF, &s);
3014: MatCreate(PETSC_COMM_SELF, &m);
3015: DMSetDefaultConstraints(*cdm, s, m);
3016: PetscSectionDestroy(&s);
3017: MatDestroy(&m);
3018: return(0);
3019: }
3023: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
3024: {
3025: DM_Plex *mesh = (DM_Plex*) dm->data;
3029: if (section) *section = mesh->coneSection;
3030: return(0);
3031: }
3035: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
3036: {
3037: DM_Plex *mesh = (DM_Plex*) dm->data;
3041: if (section) *section = mesh->supportSection;
3042: return(0);
3043: }
3047: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
3048: {
3049: DM_Plex *mesh = (DM_Plex*) dm->data;
3053: if (cones) *cones = mesh->cones;
3054: return(0);
3055: }
3059: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
3060: {
3061: DM_Plex *mesh = (DM_Plex*) dm->data;
3065: if (coneOrientations) *coneOrientations = mesh->coneOrientations;
3066: return(0);
3067: }
3069: /******************************** FEM Support **********************************/
3073: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3074: {
3075: PetscScalar *array, *vArray;
3076: const PetscInt *cone, *coneO;
3077: PetscInt pStart, pEnd, p, numPoints, size = 0, offset = 0;
3078: PetscErrorCode ierr;
3081: PetscSectionGetChart(section, &pStart, &pEnd);
3082: DMPlexGetConeSize(dm, point, &numPoints);
3083: DMPlexGetCone(dm, point, &cone);
3084: DMPlexGetConeOrientation(dm, point, &coneO);
3085: if (!values || !*values) {
3086: if ((point >= pStart) && (point < pEnd)) {
3087: PetscInt dof;
3089: PetscSectionGetDof(section, point, &dof);
3090: size += dof;
3091: }
3092: for (p = 0; p < numPoints; ++p) {
3093: const PetscInt cp = cone[p];
3094: PetscInt dof;
3096: if ((cp < pStart) || (cp >= pEnd)) continue;
3097: PetscSectionGetDof(section, cp, &dof);
3098: size += dof;
3099: }
3100: if (!values) {
3101: if (csize) *csize = size;
3102: return(0);
3103: }
3104: DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
3105: } else {
3106: array = *values;
3107: }
3108: size = 0;
3109: VecGetArray(v, &vArray);
3110: if ((point >= pStart) && (point < pEnd)) {
3111: PetscInt dof, off, d;
3112: PetscScalar *varr;
3114: PetscSectionGetDof(section, point, &dof);
3115: PetscSectionGetOffset(section, point, &off);
3116: varr = &vArray[off];
3117: for (d = 0; d < dof; ++d, ++offset) {
3118: array[offset] = varr[d];
3119: }
3120: size += dof;
3121: }
3122: for (p = 0; p < numPoints; ++p) {
3123: const PetscInt cp = cone[p];
3124: PetscInt o = coneO[p];
3125: PetscInt dof, off, d;
3126: PetscScalar *varr;
3128: if ((cp < pStart) || (cp >= pEnd)) continue;
3129: PetscSectionGetDof(section, cp, &dof);
3130: PetscSectionGetOffset(section, cp, &off);
3131: varr = &vArray[off];
3132: if (o >= 0) {
3133: for (d = 0; d < dof; ++d, ++offset) {
3134: array[offset] = varr[d];
3135: }
3136: } else {
3137: for (d = dof-1; d >= 0; --d, ++offset) {
3138: array[offset] = varr[d];
3139: }
3140: }
3141: size += dof;
3142: }
3143: VecRestoreArray(v, &vArray);
3144: if (!*values) {
3145: if (csize) *csize = size;
3146: *values = array;
3147: } else {
3148: if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3149: *csize = size;
3150: }
3151: return(0);
3152: }
3156: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3157: {
3158: PetscInt offset = 0, p;
3162: *size = 0;
3163: for (p = 0; p < numPoints*2; p += 2) {
3164: const PetscInt point = points[p];
3165: const PetscInt o = points[p+1];
3166: PetscInt dof, off, d;
3167: const PetscScalar *varr;
3169: PetscSectionGetDof(section, point, &dof);
3170: PetscSectionGetOffset(section, point, &off);
3171: varr = &vArray[off];
3172: if (o >= 0) {
3173: for (d = 0; d < dof; ++d, ++offset) array[offset] = varr[d];
3174: } else {
3175: for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
3176: }
3177: }
3178: *size = offset;
3179: return(0);
3180: }
3184: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Fields_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], PetscInt numFields, const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3185: {
3186: PetscInt offset = 0, f;
3190: *size = 0;
3191: for (f = 0; f < numFields; ++f) {
3192: PetscInt fcomp, p;
3194: PetscSectionGetFieldComponents(section, f, &fcomp);
3195: for (p = 0; p < numPoints*2; p += 2) {
3196: const PetscInt point = points[p];
3197: const PetscInt o = points[p+1];
3198: PetscInt fdof, foff, d, c;
3199: const PetscScalar *varr;
3201: PetscSectionGetFieldDof(section, point, f, &fdof);
3202: PetscSectionGetFieldOffset(section, point, f, &foff);
3203: varr = &vArray[foff];
3204: if (o >= 0) {
3205: for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
3206: } else {
3207: for (d = fdof/fcomp-1; d >= 0; --d) {
3208: for (c = 0; c < fcomp; ++c, ++offset) {
3209: array[offset] = varr[d*fcomp+c];
3210: }
3211: }
3212: }
3213: }
3214: }
3215: *size = offset;
3216: return(0);
3217: }
3221: /*@C
3222: DMPlexVecGetClosure - Get an array of the values on the closure of 'point'
3224: Not collective
3226: Input Parameters:
3227: + dm - The DM
3228: . section - The section describing the layout in v, or NULL to use the default section
3229: . v - The local vector
3230: - point - The sieve point in the DM
3232: Output Parameters:
3233: + csize - The number of values in the closure, or NULL
3234: - values - The array of values, which is a borrowed array and should not be freed
3236: Fortran Notes:
3237: Since it returns an array, this routine is only available in Fortran 90, and you must
3238: include petsc.h90 in your code.
3240: The csize argument is not present in the Fortran 90 binding since it is internal to the array.
3242: Level: intermediate
3244: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3245: @*/
3246: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3247: {
3248: PetscSection clSection;
3249: IS clPoints;
3250: PetscScalar *array, *vArray;
3251: PetscInt *points = NULL;
3252: const PetscInt *clp;
3253: PetscInt depth, numFields, numPoints, size;
3254: PetscErrorCode ierr;
3258: if (!section) {DMGetDefaultSection(dm, §ion);}
3261: DMPlexGetDepth(dm, &depth);
3262: PetscSectionGetNumFields(section, &numFields);
3263: if (depth == 1 && numFields < 2) {
3264: DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
3265: return(0);
3266: }
3267: /* Get points */
3268: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3269: if (!clPoints) {
3270: PetscInt pStart, pEnd, p, q;
3272: PetscSectionGetChart(section, &pStart, &pEnd);
3273: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3274: /* Compress out points not in the section */
3275: for (p = 0, q = 0; p < numPoints*2; p += 2) {
3276: if ((points[p] >= pStart) && (points[p] < pEnd)) {
3277: points[q*2] = points[p];
3278: points[q*2+1] = points[p+1];
3279: ++q;
3280: }
3281: }
3282: numPoints = q;
3283: } else {
3284: PetscInt dof, off;
3286: PetscSectionGetDof(clSection, point, &dof);
3287: PetscSectionGetOffset(clSection, point, &off);
3288: ISGetIndices(clPoints, &clp);
3289: numPoints = dof/2;
3290: points = (PetscInt *) &clp[off];
3291: }
3292: /* Get array */
3293: if (!values || !*values) {
3294: PetscInt asize = 0, dof, p;
3296: for (p = 0; p < numPoints*2; p += 2) {
3297: PetscSectionGetDof(section, points[p], &dof);
3298: asize += dof;
3299: }
3300: if (!values) {
3301: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3302: else {ISRestoreIndices(clPoints, &clp);}
3303: if (csize) *csize = asize;
3304: return(0);
3305: }
3306: DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
3307: } else {
3308: array = *values;
3309: }
3310: VecGetArray(v, &vArray);
3311: /* Get values */
3312: if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);}
3313: else {DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);}
3314: /* Cleanup points */
3315: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3316: else {ISRestoreIndices(clPoints, &clp);}
3317: /* Cleanup array */
3318: VecRestoreArray(v, &vArray);
3319: if (!*values) {
3320: if (csize) *csize = size;
3321: *values = array;
3322: } else {
3323: if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %D < actual size %D", *csize, size);
3324: *csize = size;
3325: }
3326: return(0);
3327: }
3331: /*@C
3332: DMPlexVecRestoreClosure - Restore the array of the values on the closure of 'point'
3334: Not collective
3336: Input Parameters:
3337: + dm - The DM
3338: . section - The section describing the layout in v, or NULL to use the default section
3339: . v - The local vector
3340: . point - The sieve point in the DM
3341: . csize - The number of values in the closure, or NULL
3342: - values - The array of values, which is a borrowed array and should not be freed
3344: Fortran Notes:
3345: Since it returns an array, this routine is only available in Fortran 90, and you must
3346: include petsc.h90 in your code.
3348: The csize argument is not present in the Fortran 90 binding since it is internal to the array.
3350: Level: intermediate
3352: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3353: @*/
3354: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3355: {
3356: PetscInt size = 0;
3360: /* Should work without recalculating size */
3361: DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
3362: return(0);
3363: }
3365: PETSC_STATIC_INLINE void add (PetscScalar *x, PetscScalar y) {*x += y;}
3366: PETSC_STATIC_INLINE void insert(PetscScalar *x, PetscScalar y) {*x = y;}
3370: 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[])
3371: {
3372: PetscInt cdof; /* The number of constraints on this point */
3373: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3374: PetscScalar *a;
3375: PetscInt off, cind = 0, k;
3376: PetscErrorCode ierr;
3379: PetscSectionGetConstraintDof(section, point, &cdof);
3380: PetscSectionGetOffset(section, point, &off);
3381: a = &array[off];
3382: if (!cdof || setBC) {
3383: if (orientation >= 0) {
3384: for (k = 0; k < dof; ++k) {
3385: fuse(&a[k], values[k]);
3386: }
3387: } else {
3388: for (k = 0; k < dof; ++k) {
3389: fuse(&a[k], values[dof-k-1]);
3390: }
3391: }
3392: } else {
3393: PetscSectionGetConstraintIndices(section, point, &cdofs);
3394: if (orientation >= 0) {
3395: for (k = 0; k < dof; ++k) {
3396: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3397: fuse(&a[k], values[k]);
3398: }
3399: } else {
3400: for (k = 0; k < dof; ++k) {
3401: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3402: fuse(&a[k], values[dof-k-1]);
3403: }
3404: }
3405: }
3406: return(0);
3407: }
3411: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3412: {
3413: PetscInt cdof; /* The number of constraints on this point */
3414: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3415: PetscScalar *a;
3416: PetscInt off, cind = 0, k;
3417: PetscErrorCode ierr;
3420: PetscSectionGetConstraintDof(section, point, &cdof);
3421: PetscSectionGetOffset(section, point, &off);
3422: a = &array[off];
3423: if (cdof) {
3424: PetscSectionGetConstraintIndices(section, point, &cdofs);
3425: if (orientation >= 0) {
3426: for (k = 0; k < dof; ++k) {
3427: if ((cind < cdof) && (k == cdofs[cind])) {
3428: fuse(&a[k], values[k]);
3429: ++cind;
3430: }
3431: }
3432: } else {
3433: for (k = 0; k < dof; ++k) {
3434: if ((cind < cdof) && (k == cdofs[cind])) {
3435: fuse(&a[k], values[dof-k-1]);
3436: ++cind;
3437: }
3438: }
3439: }
3440: }
3441: return(0);
3442: }
3446: 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[])
3447: {
3448: PetscScalar *a;
3449: PetscInt fdof, foff, fcdof, foffset = *offset;
3450: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3451: PetscInt cind = 0, k, c;
3452: PetscErrorCode ierr;
3455: PetscSectionGetFieldDof(section, point, f, &fdof);
3456: PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3457: PetscSectionGetFieldOffset(section, point, f, &foff);
3458: a = &array[foff];
3459: if (!fcdof || setBC) {
3460: if (o >= 0) {
3461: for (k = 0; k < fdof; ++k) fuse(&a[k], values[foffset+k]);
3462: } else {
3463: for (k = fdof/fcomp-1; k >= 0; --k) {
3464: for (c = 0; c < fcomp; ++c) {
3465: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3466: }
3467: }
3468: }
3469: } else {
3470: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3471: if (o >= 0) {
3472: for (k = 0; k < fdof; ++k) {
3473: if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
3474: fuse(&a[k], values[foffset+k]);
3475: }
3476: } else {
3477: for (k = fdof/fcomp-1; k >= 0; --k) {
3478: for (c = 0; c < fcomp; ++c) {
3479: if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
3480: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3481: }
3482: }
3483: }
3484: }
3485: *offset += fdof;
3486: return(0);
3487: }
3491: 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[])
3492: {
3493: PetscScalar *a;
3494: PetscInt fdof, foff, fcdof, foffset = *offset;
3495: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3496: PetscInt cind = 0, k, c;
3497: PetscErrorCode ierr;
3500: PetscSectionGetFieldDof(section, point, f, &fdof);
3501: PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3502: PetscSectionGetFieldOffset(section, point, f, &foff);
3503: a = &array[foff];
3504: if (fcdof) {
3505: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3506: if (o >= 0) {
3507: for (k = 0; k < fdof; ++k) {
3508: if ((cind < fcdof) && (k == fcdofs[cind])) {
3509: fuse(&a[k], values[foffset+k]);
3510: ++cind;
3511: }
3512: }
3513: } else {
3514: for (k = fdof/fcomp-1; k >= 0; --k) {
3515: for (c = 0; c < fcomp; ++c) {
3516: if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {
3517: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3518: ++cind;
3519: }
3520: }
3521: }
3522: }
3523: }
3524: *offset += fdof;
3525: return(0);
3526: }
3530: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3531: {
3532: PetscScalar *array;
3533: const PetscInt *cone, *coneO;
3534: PetscInt pStart, pEnd, p, numPoints, off, dof;
3535: PetscErrorCode ierr;
3538: PetscSectionGetChart(section, &pStart, &pEnd);
3539: DMPlexGetConeSize(dm, point, &numPoints);
3540: DMPlexGetCone(dm, point, &cone);
3541: DMPlexGetConeOrientation(dm, point, &coneO);
3542: VecGetArray(v, &array);
3543: for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
3544: const PetscInt cp = !p ? point : cone[p-1];
3545: const PetscInt o = !p ? 0 : coneO[p-1];
3547: if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
3548: PetscSectionGetDof(section, cp, &dof);
3549: /* ADD_VALUES */
3550: {
3551: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3552: PetscScalar *a;
3553: PetscInt cdof, coff, cind = 0, k;
3555: PetscSectionGetConstraintDof(section, cp, &cdof);
3556: PetscSectionGetOffset(section, cp, &coff);
3557: a = &array[coff];
3558: if (!cdof) {
3559: if (o >= 0) {
3560: for (k = 0; k < dof; ++k) {
3561: a[k] += values[off+k];
3562: }
3563: } else {
3564: for (k = 0; k < dof; ++k) {
3565: a[k] += values[off+dof-k-1];
3566: }
3567: }
3568: } else {
3569: PetscSectionGetConstraintIndices(section, cp, &cdofs);
3570: if (o >= 0) {
3571: for (k = 0; k < dof; ++k) {
3572: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3573: a[k] += values[off+k];
3574: }
3575: } else {
3576: for (k = 0; k < dof; ++k) {
3577: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3578: a[k] += values[off+dof-k-1];
3579: }
3580: }
3581: }
3582: }
3583: }
3584: VecRestoreArray(v, &array);
3585: return(0);
3586: }
3590: /*@C
3591: DMPlexVecSetClosure - Set an array of the values on the closure of 'point'
3593: Not collective
3595: Input Parameters:
3596: + dm - The DM
3597: . section - The section describing the layout in v, or NULL to use the default section
3598: . v - The local vector
3599: . point - The sieve point in the DM
3600: . values - The array of values
3601: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
3603: Fortran Notes:
3604: This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
3606: Level: intermediate
3608: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
3609: @*/
3610: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3611: {
3612: PetscSection clSection;
3613: IS clPoints;
3614: PetscScalar *array;
3615: PetscInt *points = NULL;
3616: const PetscInt *clp;
3617: PetscInt depth, numFields, numPoints, p;
3618: PetscErrorCode ierr;
3622: if (!section) {DMGetDefaultSection(dm, §ion);}
3625: DMPlexGetDepth(dm, &depth);
3626: PetscSectionGetNumFields(section, &numFields);
3627: if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
3628: DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
3629: return(0);
3630: }
3631: /* Get points */
3632: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3633: if (!clPoints) {
3634: PetscInt pStart, pEnd, q;
3636: PetscSectionGetChart(section, &pStart, &pEnd);
3637: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3638: /* Compress out points not in the section */
3639: for (p = 0, q = 0; p < numPoints*2; p += 2) {
3640: if ((points[p] >= pStart) && (points[p] < pEnd)) {
3641: points[q*2] = points[p];
3642: points[q*2+1] = points[p+1];
3643: ++q;
3644: }
3645: }
3646: numPoints = q;
3647: } else {
3648: PetscInt dof, off;
3650: PetscSectionGetDof(clSection, point, &dof);
3651: PetscSectionGetOffset(clSection, point, &off);
3652: ISGetIndices(clPoints, &clp);
3653: numPoints = dof/2;
3654: points = (PetscInt *) &clp[off];
3655: }
3656: /* Get array */
3657: VecGetArray(v, &array);
3658: /* Get values */
3659: if (numFields > 0) {
3660: PetscInt offset = 0, fcomp, f;
3661: for (f = 0; f < numFields; ++f) {
3662: PetscSectionGetFieldComponents(section, f, &fcomp);
3663: switch (mode) {
3664: case INSERT_VALUES:
3665: for (p = 0; p < numPoints*2; p += 2) {
3666: const PetscInt point = points[p];
3667: const PetscInt o = points[p+1];
3668: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3669: } break;
3670: case INSERT_ALL_VALUES:
3671: for (p = 0; p < numPoints*2; p += 2) {
3672: const PetscInt point = points[p];
3673: const PetscInt o = points[p+1];
3674: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3675: } break;
3676: case INSERT_BC_VALUES:
3677: for (p = 0; p < numPoints*2; p += 2) {
3678: const PetscInt point = points[p];
3679: const PetscInt o = points[p+1];
3680: updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3681: } break;
3682: case ADD_VALUES:
3683: for (p = 0; p < numPoints*2; p += 2) {
3684: const PetscInt point = points[p];
3685: const PetscInt o = points[p+1];
3686: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3687: } break;
3688: case ADD_ALL_VALUES:
3689: for (p = 0; p < numPoints*2; p += 2) {
3690: const PetscInt point = points[p];
3691: const PetscInt o = points[p+1];
3692: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3693: } break;
3694: default:
3695: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3696: }
3697: }
3698: } else {
3699: PetscInt dof, off;
3701: switch (mode) {
3702: case INSERT_VALUES:
3703: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3704: PetscInt o = points[p+1];
3705: PetscSectionGetDof(section, points[p], &dof);
3706: updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
3707: } break;
3708: case INSERT_ALL_VALUES:
3709: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3710: PetscInt o = points[p+1];
3711: PetscSectionGetDof(section, points[p], &dof);
3712: updatePoint_private(section, points[p], dof, insert, PETSC_TRUE, o, &values[off], array);
3713: } break;
3714: case INSERT_BC_VALUES:
3715: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3716: PetscInt o = points[p+1];
3717: PetscSectionGetDof(section, points[p], &dof);
3718: updatePointBC_private(section, points[p], dof, insert, o, &values[off], array);
3719: } break;
3720: case ADD_VALUES:
3721: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3722: PetscInt o = points[p+1];
3723: PetscSectionGetDof(section, points[p], &dof);
3724: updatePoint_private(section, points[p], dof, add, PETSC_FALSE, o, &values[off], array);
3725: } break;
3726: case ADD_ALL_VALUES:
3727: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3728: PetscInt o = points[p+1];
3729: PetscSectionGetDof(section, points[p], &dof);
3730: updatePoint_private(section, points[p], dof, add, PETSC_TRUE, o, &values[off], array);
3731: } break;
3732: default:
3733: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3734: }
3735: }
3736: /* Cleanup points */
3737: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3738: else {ISRestoreIndices(clPoints, &clp);}
3739: /* Cleanup array */
3740: VecRestoreArray(v, &array);
3741: return(0);
3742: }
3746: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
3747: {
3748: PetscSection clSection;
3749: IS clPoints;
3750: PetscScalar *array;
3751: PetscInt *points = NULL;
3752: const PetscInt *clp;
3753: PetscInt numFields, numPoints, p;
3754: PetscInt offset = 0, fcomp, f;
3755: PetscErrorCode ierr;
3759: if (!section) {DMGetDefaultSection(dm, §ion);}
3762: PetscSectionGetNumFields(section, &numFields);
3763: /* Get points */
3764: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3765: if (!clPoints) {
3766: PetscInt pStart, pEnd, q;
3768: PetscSectionGetChart(section, &pStart, &pEnd);
3769: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3770: /* Compress out points not in the section */
3771: for (p = 0, q = 0; p < numPoints*2; p += 2) {
3772: if ((points[p] >= pStart) && (points[p] < pEnd)) {
3773: points[q*2] = points[p];
3774: points[q*2+1] = points[p+1];
3775: ++q;
3776: }
3777: }
3778: numPoints = q;
3779: } else {
3780: PetscInt dof, off;
3782: PetscSectionGetDof(clSection, point, &dof);
3783: PetscSectionGetOffset(clSection, point, &off);
3784: ISGetIndices(clPoints, &clp);
3785: numPoints = dof/2;
3786: points = (PetscInt *) &clp[off];
3787: }
3788: /* Get array */
3789: VecGetArray(v, &array);
3790: /* Get values */
3791: for (f = 0; f < numFields; ++f) {
3792: PetscSectionGetFieldComponents(section, f, &fcomp);
3793: if (!fieldActive[f]) {
3794: for (p = 0; p < numPoints*2; p += 2) {
3795: PetscInt fdof;
3796: PetscSectionGetFieldDof(section, points[p], f, &fdof);
3797: offset += fdof;
3798: }
3799: continue;
3800: }
3801: switch (mode) {
3802: case INSERT_VALUES:
3803: for (p = 0; p < numPoints*2; p += 2) {
3804: const PetscInt point = points[p];
3805: const PetscInt o = points[p+1];
3806: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3807: } break;
3808: case INSERT_ALL_VALUES:
3809: for (p = 0; p < numPoints*2; p += 2) {
3810: const PetscInt point = points[p];
3811: const PetscInt o = points[p+1];
3812: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3813: } break;
3814: case INSERT_BC_VALUES:
3815: for (p = 0; p < numPoints*2; p += 2) {
3816: const PetscInt point = points[p];
3817: const PetscInt o = points[p+1];
3818: updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3819: } break;
3820: case ADD_VALUES:
3821: for (p = 0; p < numPoints*2; p += 2) {
3822: const PetscInt point = points[p];
3823: const PetscInt o = points[p+1];
3824: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3825: } break;
3826: case ADD_ALL_VALUES:
3827: for (p = 0; p < numPoints*2; p += 2) {
3828: const PetscInt point = points[p];
3829: const PetscInt o = points[p+1];
3830: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3831: } break;
3832: default:
3833: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %d", mode);
3834: }
3835: }
3836: /* Cleanup points */
3837: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3838: else {ISRestoreIndices(clPoints, &clp);}
3839: /* Cleanup array */
3840: VecRestoreArray(v, &array);
3841: return(0);
3842: }
3846: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
3847: {
3848: PetscMPIInt rank;
3849: PetscInt i, j;
3853: MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
3854: PetscViewerASCIIPrintf(viewer, "[%d]mat for sieve point %D\n", rank, point);
3855: for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
3856: for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
3857: numCIndices = numCIndices ? numCIndices : numRIndices;
3858: for (i = 0; i < numRIndices; i++) {
3859: PetscViewerASCIIPrintf(viewer, "[%d]", rank);
3860: for (j = 0; j < numCIndices; j++) {
3861: #if defined(PETSC_USE_COMPLEX)
3862: PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
3863: #else
3864: PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
3865: #endif
3866: }
3867: PetscViewerASCIIPrintf(viewer, "\n");
3868: }
3869: return(0);
3870: }
3874: /* . off - The global offset of this point */
3875: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
3876: {
3877: PetscInt dof; /* The number of unknowns on this point */
3878: PetscInt cdof; /* The number of constraints on this point */
3879: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3880: PetscInt cind = 0, k;
3881: PetscErrorCode ierr;
3884: PetscSectionGetDof(section, point, &dof);
3885: PetscSectionGetConstraintDof(section, point, &cdof);
3886: if (!cdof || setBC) {
3887: if (orientation >= 0) {
3888: for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
3889: } else {
3890: for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
3891: }
3892: } else {
3893: PetscSectionGetConstraintIndices(section, point, &cdofs);
3894: if (orientation >= 0) {
3895: for (k = 0; k < dof; ++k) {
3896: if ((cind < cdof) && (k == cdofs[cind])) {
3897: /* Insert check for returning constrained indices */
3898: indices[*loff+k] = -(off+k+1);
3899: ++cind;
3900: } else {
3901: indices[*loff+k] = off+k-cind;
3902: }
3903: }
3904: } else {
3905: for (k = 0; k < dof; ++k) {
3906: if ((cind < cdof) && (k == cdofs[cind])) {
3907: /* Insert check for returning constrained indices */
3908: indices[*loff+dof-k-1] = -(off+k+1);
3909: ++cind;
3910: } else {
3911: indices[*loff+dof-k-1] = off+k-cind;
3912: }
3913: }
3914: }
3915: }
3916: *loff += dof;
3917: return(0);
3918: }
3922: /* . off - The global offset of this point */
3923: PetscErrorCode indicesPointFields_private(PetscSection section, PetscInt point, PetscInt off, PetscInt foffs[], PetscBool setBC, PetscInt orientation, PetscInt indices[])
3924: {
3925: PetscInt numFields, foff, f;
3929: PetscSectionGetNumFields(section, &numFields);
3930: for (f = 0, foff = 0; f < numFields; ++f) {
3931: PetscInt fdof, fcomp, cfdof;
3932: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3933: PetscInt cind = 0, k, c;
3935: PetscSectionGetFieldComponents(section, f, &fcomp);
3936: PetscSectionGetFieldDof(section, point, f, &fdof);
3937: PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
3938: if (!cfdof || setBC) {
3939: if (orientation >= 0) {
3940: for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
3941: } else {
3942: for (k = fdof/fcomp-1; k >= 0; --k) {
3943: for (c = 0; c < fcomp; ++c) {
3944: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
3945: }
3946: }
3947: }
3948: } else {
3949: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3950: if (orientation >= 0) {
3951: for (k = 0; k < fdof; ++k) {
3952: if ((cind < cfdof) && (k == fcdofs[cind])) {
3953: indices[foffs[f]+k] = -(off+foff+k+1);
3954: ++cind;
3955: } else {
3956: indices[foffs[f]+k] = off+foff+k-cind;
3957: }
3958: }
3959: } else {
3960: for (k = fdof/fcomp-1; k >= 0; --k) {
3961: for (c = 0; c < fcomp; ++c) {
3962: if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
3963: indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
3964: ++cind;
3965: } else {
3966: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
3967: }
3968: }
3969: }
3970: }
3971: }
3972: foff += (setBC ? fdof : (fdof - cfdof));
3973: foffs[f] += fdof;
3974: }
3975: return(0);
3976: }
3980: 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)
3981: {
3982: Mat cMat;
3983: PetscSection aSec, cSec;
3984: IS aIS;
3985: PetscInt aStart = -1, aEnd = -1;
3986: const PetscInt *anchors;
3987: PetscInt numFields, f, p, q, newP = 0;
3988: PetscInt newNumPoints = 0, newNumIndices = 0;
3989: PetscInt *newPoints, *indices, *newIndices;
3990: PetscInt maxAnchor, maxDof;
3991: PetscInt newOffsets[32];
3992: PetscInt *pointMatOffsets[32];
3993: PetscInt *newPointOffsets[32];
3994: PetscScalar *pointMat[32];
3995: PetscScalar *newValues=NULL,*tmpValues;
3996: PetscBool anyConstrained = PETSC_FALSE;
3997: PetscErrorCode ierr;
4002: PetscSectionGetNumFields(section, &numFields);
4004: DMPlexGetAnchors(dm,&aSec,&aIS);
4005: /* if there are point-to-point constraints */
4006: if (aSec) {
4007: PetscMemzero(newOffsets, 32 * sizeof(PetscInt));
4008: ISGetIndices(aIS,&anchors);
4009: PetscSectionGetChart(aSec,&aStart,&aEnd);
4010: /* figure out how many points are going to be in the new element matrix
4011: * (we allow double counting, because it's all just going to be summed
4012: * into the global matrix anyway) */
4013: for (p = 0; p < 2*numPoints; p+=2) {
4014: PetscInt b = points[p];
4015: PetscInt bDof = 0, bSecDof;
4017: PetscSectionGetDof(section,b,&bSecDof);
4018: if (!bSecDof) {
4019: continue;
4020: }
4021: if (b >= aStart && b < aEnd) {
4022: PetscSectionGetDof(aSec,b,&bDof);
4023: }
4024: if (bDof) {
4025: /* this point is constrained */
4026: /* it is going to be replaced by its anchors */
4027: PetscInt bOff, q;
4029: anyConstrained = PETSC_TRUE;
4030: newNumPoints += bDof;
4031: PetscSectionGetOffset(aSec,b,&bOff);
4032: for (q = 0; q < bDof; q++) {
4033: PetscInt a = anchors[bOff + q];
4034: PetscInt aDof;
4036: PetscSectionGetDof(section,a,&aDof);
4037: newNumIndices += aDof;
4038: for (f = 0; f < numFields; ++f) {
4039: PetscInt fDof;
4041: PetscSectionGetFieldDof(section, a, f, &fDof);
4042: newOffsets[f+1] += fDof;
4043: }
4044: }
4045: }
4046: else {
4047: /* this point is not constrained */
4048: newNumPoints++;
4049: newNumIndices += bSecDof;
4050: for (f = 0; f < numFields; ++f) {
4051: PetscInt fDof;
4053: PetscSectionGetFieldDof(section, b, f, &fDof);
4054: newOffsets[f+1] += fDof;
4055: }
4056: }
4057: }
4058: }
4059: if (!anyConstrained) {
4060: if (outNumPoints) *outNumPoints = 0;
4061: if (outNumIndices) *outNumIndices = 0;
4062: if (outPoints) *outPoints = NULL;
4063: if (outValues) *outValues = NULL;
4064: if (aSec) {ISRestoreIndices(aIS,&anchors);}
4065: return(0);
4066: }
4068: if (outNumPoints) *outNumPoints = newNumPoints;
4069: if (outNumIndices) *outNumIndices = newNumIndices;
4071: for (f = 1; f < numFields; ++f) newOffsets[f+1] += newOffsets[f];
4073: if (!outPoints && !outValues) {
4074: if (offsets) {
4075: for (f = 0; f <= numFields; f++) {
4076: offsets[f] = newOffsets[f];
4077: }
4078: }
4079: if (aSec) {ISRestoreIndices(aIS,&anchors);}
4080: return(0);
4081: }
4083: if (numFields && newOffsets[numFields] != newNumIndices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid size for closure %D should be %D", newOffsets[numFields], newNumIndices);
4085: DMGetDefaultConstraints(dm, &cSec, &cMat);
4087: /* workspaces */
4088: if (numFields) {
4089: for (f = 0; f < numFields; f++) {
4090: DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4091: DMGetWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4092: }
4093: }
4094: else {
4095: DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4096: DMGetWorkArray(dm,numPoints,PETSC_INT,&newPointOffsets[0]);
4097: }
4099: /* get workspaces for the point-to-point matrices */
4100: if (numFields) {
4101: PetscInt totalOffset, totalMatOffset;
4103: for (p = 0; p < numPoints; p++) {
4104: PetscInt b = points[2*p];
4105: PetscInt bDof = 0, bSecDof;
4107: PetscSectionGetDof(section,b,&bSecDof);
4108: if (!bSecDof) {
4109: for (f = 0; f < numFields; f++) {
4110: newPointOffsets[f][p + 1] = 0;
4111: pointMatOffsets[f][p + 1] = 0;
4112: }
4113: continue;
4114: }
4115: if (b >= aStart && b < aEnd) {
4116: PetscSectionGetDof(aSec, b, &bDof);
4117: }
4118: if (bDof) {
4119: for (f = 0; f < numFields; f++) {
4120: PetscInt fDof, q, bOff, allFDof = 0;
4122: PetscSectionGetFieldDof(section, b, f, &fDof);
4123: PetscSectionGetOffset(aSec, b, &bOff);
4124: for (q = 0; q < bDof; q++) {
4125: PetscInt a = anchors[bOff + q];
4126: PetscInt aFDof;
4128: PetscSectionGetFieldDof(section, a, f, &aFDof);
4129: allFDof += aFDof;
4130: }
4131: newPointOffsets[f][p+1] = allFDof;
4132: pointMatOffsets[f][p+1] = fDof * allFDof;
4133: }
4134: }
4135: else {
4136: for (f = 0; f < numFields; f++) {
4137: PetscInt fDof;
4139: PetscSectionGetFieldDof(section, b, f, &fDof);
4140: newPointOffsets[f][p+1] = fDof;
4141: pointMatOffsets[f][p+1] = 0;
4142: }
4143: }
4144: }
4145: for (f = 0, totalOffset = 0, totalMatOffset = 0; f < numFields; f++) {
4146: newPointOffsets[f][0] = totalOffset;
4147: pointMatOffsets[f][0] = totalMatOffset;
4148: for (p = 0; p < numPoints; p++) {
4149: newPointOffsets[f][p+1] += newPointOffsets[f][p];
4150: pointMatOffsets[f][p+1] += pointMatOffsets[f][p];
4151: }
4152: totalOffset = newPointOffsets[f][numPoints];
4153: totalMatOffset = pointMatOffsets[f][numPoints];
4154: DMGetWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4155: }
4156: }
4157: else {
4158: for (p = 0; p < numPoints; p++) {
4159: PetscInt b = points[2*p];
4160: PetscInt bDof = 0, bSecDof;
4162: PetscSectionGetDof(section,b,&bSecDof);
4163: if (!bSecDof) {
4164: newPointOffsets[0][p + 1] = 0;
4165: pointMatOffsets[0][p + 1] = 0;
4166: continue;
4167: }
4168: if (b >= aStart && b < aEnd) {
4169: PetscSectionGetDof(aSec, b, &bDof);
4170: }
4171: if (bDof) {
4172: PetscInt bOff, q, allDof = 0;
4174: PetscSectionGetOffset(aSec, b, &bOff);
4175: for (q = 0; q < bDof; q++) {
4176: PetscInt a = anchors[bOff + q], aDof;
4178: PetscSectionGetDof(section, a, &aDof);
4179: allDof += aDof;
4180: }
4181: newPointOffsets[0][p+1] = allDof;
4182: pointMatOffsets[0][p+1] = bSecDof * allDof;
4183: }
4184: else {
4185: newPointOffsets[0][p+1] = bSecDof;
4186: pointMatOffsets[0][p+1] = 0;
4187: }
4188: }
4189: newPointOffsets[0][0] = 0;
4190: pointMatOffsets[0][0] = 0;
4191: for (p = 0; p < numPoints; p++) {
4192: newPointOffsets[0][p+1] += newPointOffsets[0][p];
4193: pointMatOffsets[0][p+1] += pointMatOffsets[0][p];
4194: }
4195: DMGetWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4196: }
4198: /* output arrays */
4199: DMGetWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4201: /* get the point-to-point matrices; construct newPoints */
4202: PetscSectionGetMaxDof(aSec, &maxAnchor);
4203: PetscSectionGetMaxDof(section, &maxDof);
4204: DMGetWorkArray(dm,maxDof,PETSC_INT,&indices);
4205: DMGetWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4206: if (numFields) {
4207: for (p = 0, newP = 0; p < numPoints; p++) {
4208: PetscInt b = points[2*p];
4209: PetscInt o = points[2*p+1];
4210: PetscInt bDof = 0, bSecDof;
4212: PetscSectionGetDof(section, b, &bSecDof);
4213: if (!bSecDof) {
4214: continue;
4215: }
4216: if (b >= aStart && b < aEnd) {
4217: PetscSectionGetDof(aSec, b, &bDof);
4218: }
4219: if (bDof) {
4220: PetscInt fStart[32], fEnd[32], fAnchorStart[32], fAnchorEnd[32], bOff, q;
4222: fStart[0] = 0;
4223: fEnd[0] = 0;
4224: for (f = 0; f < numFields; f++) {
4225: PetscInt fDof;
4227: PetscSectionGetFieldDof(cSec, b, f, &fDof);
4228: fStart[f+1] = fStart[f] + fDof;
4229: fEnd[f+1] = fStart[f+1];
4230: }
4231: PetscSectionGetOffset(cSec, b, &bOff);
4232: indicesPointFields_private(cSec, b, bOff, fEnd, PETSC_TRUE, o, indices);
4234: fAnchorStart[0] = 0;
4235: fAnchorEnd[0] = 0;
4236: for (f = 0; f < numFields; f++) {
4237: PetscInt fDof = newPointOffsets[f][p + 1] - newPointOffsets[f][p];
4239: fAnchorStart[f+1] = fAnchorStart[f] + fDof;
4240: fAnchorEnd[f+1] = fAnchorStart[f + 1];
4241: }
4242: PetscSectionGetOffset(aSec, b, &bOff);
4243: for (q = 0; q < bDof; q++) {
4244: PetscInt a = anchors[bOff + q], aOff;
4246: /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4247: newPoints[2*(newP + q)] = a;
4248: newPoints[2*(newP + q) + 1] = 0;
4249: PetscSectionGetOffset(section, a, &aOff);
4250: indicesPointFields_private(section, a, aOff, fAnchorEnd, PETSC_TRUE, 0, newIndices);
4251: }
4252: newP += bDof;
4254: if (outValues) {
4255: /* get the point-to-point submatrix */
4256: for (f = 0; f < numFields; f++) {
4257: MatGetValues(cMat,fEnd[f]-fStart[f],indices + fStart[f],fAnchorEnd[f] - fAnchorStart[f],newIndices + fAnchorStart[f],pointMat[f] + pointMatOffsets[f][p]);
4258: }
4259: }
4260: }
4261: else {
4262: newPoints[2 * newP] = b;
4263: newPoints[2 * newP + 1] = o;
4264: newP++;
4265: }
4266: }
4267: } else {
4268: for (p = 0; p < numPoints; p++) {
4269: PetscInt b = points[2*p];
4270: PetscInt o = points[2*p+1];
4271: PetscInt bDof = 0, bSecDof;
4273: PetscSectionGetDof(section, b, &bSecDof);
4274: if (!bSecDof) {
4275: continue;
4276: }
4277: if (b >= aStart && b < aEnd) {
4278: PetscSectionGetDof(aSec, b, &bDof);
4279: }
4280: if (bDof) {
4281: PetscInt bEnd = 0, bAnchorEnd = 0, bOff;
4283: PetscSectionGetOffset(cSec, b, &bOff);
4284: indicesPoint_private(cSec, b, bOff, &bEnd, PETSC_TRUE, o, indices);
4286: PetscSectionGetOffset (aSec, b, &bOff);
4287: for (q = 0; q < bDof; q++) {
4288: PetscInt a = anchors[bOff + q], aOff;
4290: /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4292: newPoints[2*(newP + q)] = a;
4293: newPoints[2*(newP + q) + 1] = 0;
4294: PetscSectionGetOffset(section, a, &aOff);
4295: indicesPoint_private(section, a, aOff, &bAnchorEnd, PETSC_TRUE, 0, newIndices);
4296: }
4297: newP += bDof;
4299: /* get the point-to-point submatrix */
4300: if (outValues) {
4301: MatGetValues(cMat,bEnd,indices,bAnchorEnd,newIndices,pointMat[0] + pointMatOffsets[0][p]);
4302: }
4303: }
4304: else {
4305: newPoints[2 * newP] = b;
4306: newPoints[2 * newP + 1] = o;
4307: newP++;
4308: }
4309: }
4310: }
4312: if (outValues) {
4313: DMGetWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4314: PetscMemzero(tmpValues,newNumIndices*numIndices*sizeof(*tmpValues));
4315: /* multiply constraints on the right */
4316: if (numFields) {
4317: for (f = 0; f < numFields; f++) {
4318: PetscInt oldOff = offsets[f];
4320: for (p = 0; p < numPoints; p++) {
4321: PetscInt cStart = newPointOffsets[f][p];
4322: PetscInt b = points[2 * p];
4323: PetscInt c, r, k;
4324: PetscInt dof;
4326: PetscSectionGetFieldDof(section,b,f,&dof);
4327: if (!dof) {
4328: continue;
4329: }
4330: if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4331: PetscInt nCols = newPointOffsets[f][p+1]-cStart;
4332: const PetscScalar *mat = pointMat[f] + pointMatOffsets[f][p];
4334: for (r = 0; r < numIndices; r++) {
4335: for (c = 0; c < nCols; c++) {
4336: for (k = 0; k < dof; k++) {
4337: tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4338: }
4339: }
4340: }
4341: }
4342: else {
4343: /* copy this column as is */
4344: for (r = 0; r < numIndices; r++) {
4345: for (c = 0; c < dof; c++) {
4346: tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4347: }
4348: }
4349: }
4350: oldOff += dof;
4351: }
4352: }
4353: }
4354: else {
4355: PetscInt oldOff = 0;
4356: for (p = 0; p < numPoints; p++) {
4357: PetscInt cStart = newPointOffsets[0][p];
4358: PetscInt b = points[2 * p];
4359: PetscInt c, r, k;
4360: PetscInt dof;
4362: PetscSectionGetDof(section,b,&dof);
4363: if (!dof) {
4364: continue;
4365: }
4366: if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4367: PetscInt nCols = newPointOffsets[0][p+1]-cStart;
4368: const PetscScalar *mat = pointMat[0] + pointMatOffsets[0][p];
4370: for (r = 0; r < numIndices; r++) {
4371: for (c = 0; c < nCols; c++) {
4372: for (k = 0; k < dof; k++) {
4373: tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4374: }
4375: }
4376: }
4377: }
4378: else {
4379: /* copy this column as is */
4380: for (r = 0; r < numIndices; r++) {
4381: for (c = 0; c < dof; c++) {
4382: tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4383: }
4384: }
4385: }
4386: oldOff += dof;
4387: }
4388: }
4390: if (multiplyLeft) {
4391: DMGetWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4392: PetscMemzero(newValues,newNumIndices*newNumIndices*sizeof(*newValues));
4393: /* multiply constraints transpose on the left */
4394: if (numFields) {
4395: for (f = 0; f < numFields; f++) {
4396: PetscInt oldOff = offsets[f];
4398: for (p = 0; p < numPoints; p++) {
4399: PetscInt rStart = newPointOffsets[f][p];
4400: PetscInt b = points[2 * p];
4401: PetscInt c, r, k;
4402: PetscInt dof;
4404: PetscSectionGetFieldDof(section,b,f,&dof);
4405: if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4406: PetscInt nRows = newPointOffsets[f][p+1]-rStart;
4407: const PetscScalar *PETSC_RESTRICT mat = pointMat[f] + pointMatOffsets[f][p];
4409: for (r = 0; r < nRows; r++) {
4410: for (c = 0; c < newNumIndices; c++) {
4411: for (k = 0; k < dof; k++) {
4412: newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4413: }
4414: }
4415: }
4416: }
4417: else {
4418: /* copy this row as is */
4419: for (r = 0; r < dof; r++) {
4420: for (c = 0; c < newNumIndices; c++) {
4421: newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4422: }
4423: }
4424: }
4425: oldOff += dof;
4426: }
4427: }
4428: }
4429: else {
4430: PetscInt oldOff = 0;
4432: for (p = 0; p < numPoints; p++) {
4433: PetscInt rStart = newPointOffsets[0][p];
4434: PetscInt b = points[2 * p];
4435: PetscInt c, r, k;
4436: PetscInt dof;
4438: PetscSectionGetDof(section,b,&dof);
4439: if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4440: PetscInt nRows = newPointOffsets[0][p+1]-rStart;
4441: const PetscScalar *PETSC_RESTRICT mat = pointMat[0] + pointMatOffsets[0][p];
4443: for (r = 0; r < nRows; r++) {
4444: for (c = 0; c < newNumIndices; c++) {
4445: for (k = 0; k < dof; k++) {
4446: newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4447: }
4448: }
4449: }
4450: }
4451: else {
4452: /* copy this row as is */
4453: for (r = 0; r < dof; r++) {
4454: for (c = 0; c < newNumIndices; c++) {
4455: newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4456: }
4457: }
4458: }
4459: oldOff += dof;
4460: }
4461: }
4463: DMRestoreWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4464: }
4465: else {
4466: newValues = tmpValues;
4467: }
4468: }
4470: /* clean up */
4471: DMRestoreWorkArray(dm,maxDof,PETSC_INT,&indices);
4472: DMRestoreWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4474: if (numFields) {
4475: for (f = 0; f < numFields; f++) {
4476: DMRestoreWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4477: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4478: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4479: }
4480: }
4481: else {
4482: DMRestoreWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4483: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4484: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[0]);
4485: }
4486: ISRestoreIndices(aIS,&anchors);
4488: /* output */
4489: if (outPoints) {
4490: *outPoints = newPoints;
4491: }
4492: else {
4493: DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4494: }
4495: if (outValues) {
4496: *outValues = newValues;
4497: }
4498: else {
4499: DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4500: }
4501: for (f = 0; f <= numFields; f++) {
4502: offsets[f] = newOffsets[f];
4503: }
4504: return(0);
4505: }
4509: PetscErrorCode DMPlexGetClosureIndices(DM dm, PetscSection section, PetscSection globalSection, PetscInt point, PetscInt *numIndices, PetscInt **indices, PetscInt *outOffsets)
4510: {
4511: PetscSection clSection;
4512: IS clPoints;
4513: const PetscInt *clp;
4514: PetscInt *points = NULL, *pointsNew;
4515: PetscInt numPoints, numPointsNew;
4516: PetscInt offsets[32];
4517: PetscInt Nf, Nind, NindNew, off, globalOff, f, p;
4518: PetscErrorCode ierr;
4526: PetscSectionGetNumFields(section, &Nf);
4527: if (Nf > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", Nf);
4528: PetscMemzero(offsets, 32 * sizeof(PetscInt));
4529: /* Get points in closure */
4530: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4531: if (!clPoints) {
4532: PetscInt pStart, pEnd, q;
4534: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4535: /* Compress out points not in the section */
4536: PetscSectionGetChart(section, &pStart, &pEnd);
4537: for (p = 0, q = 0; p < numPoints*2; p += 2) {
4538: if ((points[p] >= pStart) && (points[p] < pEnd)) {
4539: points[q*2] = points[p];
4540: points[q*2+1] = points[p+1];
4541: ++q;
4542: }
4543: }
4544: numPoints = q;
4545: } else {
4546: PetscInt dof, off;
4548: PetscSectionGetDof(clSection, point, &dof);
4549: numPoints = dof/2;
4550: PetscSectionGetOffset(clSection, point, &off);
4551: ISGetIndices(clPoints, &clp);
4552: points = (PetscInt *) &clp[off];
4553: }
4554: /* Get number of indices and indices per field */
4555: for (p = 0, Nind = 0; p < numPoints*2; p += 2) {
4556: PetscInt dof, fdof;
4558: PetscSectionGetDof(section, points[p], &dof);
4559: for (f = 0; f < Nf; ++f) {
4560: PetscSectionGetFieldDof(section, points[p], f, &fdof);
4561: offsets[f+1] += fdof;
4562: }
4563: Nind += dof;
4564: }
4565: for (f = 1; f < Nf; ++f) offsets[f+1] += offsets[f];
4566: if (Nf && offsets[Nf] != Nind) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[Nf], Nind);
4567: /* Correct for hanging node constraints */
4568: {
4569: DMPlexAnchorsModifyMat(dm, section, numPoints, Nind, points, NULL, &numPointsNew, &NindNew, &pointsNew, NULL, offsets, PETSC_TRUE);
4570: if (numPointsNew) {
4571: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4572: else {ISRestoreIndices(clPoints, &clp);}
4573: numPoints = numPointsNew;
4574: Nind = NindNew;
4575: points = pointsNew;
4576: }
4577: }
4578: /* Calculate indices */
4579: DMGetWorkArray(dm, Nind, PETSC_INT, indices);
4580: if (Nf) {
4581: if (outOffsets) {
4582: PetscInt f;
4584: for (f = 0; f <= Nf; f++) {
4585: outOffsets[f] = offsets[f];
4586: }
4587: }
4588: for (p = 0; p < numPoints*2; p += 2) {
4589: PetscInt o = points[p+1];
4590: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4591: indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, *indices);
4592: }
4593: } else {
4594: for (p = 0, off = 0; p < numPoints*2; p += 2) {
4595: PetscInt o = points[p+1];
4596: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4597: indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, *indices);
4598: }
4599: }
4600: /* Cleanup points */
4601: if (numPointsNew) {
4602: DMRestoreWorkArray(dm, 2*numPointsNew, PETSC_INT, &pointsNew);
4603: } else {
4604: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4605: else {ISRestoreIndices(clPoints, &clp);}
4606: }
4607: if (numIndices) *numIndices = Nind;
4608: return(0);
4609: }
4613: PetscErrorCode DMPlexRestoreClosureIndices(DM dm, PetscSection section, PetscSection globalSection, PetscInt point, PetscInt *numIndices, PetscInt **indices,PetscInt *outOffsets)
4614: {
4620: DMRestoreWorkArray(dm, 0, PETSC_INT, indices);
4621: return(0);
4622: }
4626: /*@C
4627: DMPlexMatSetClosure - Set an array of the values on the closure of 'point'
4629: Not collective
4631: Input Parameters:
4632: + dm - The DM
4633: . section - The section describing the layout in v, or NULL to use the default section
4634: . globalSection - The section describing the layout in v, or NULL to use the default global section
4635: . A - The matrix
4636: . point - The sieve point in the DM
4637: . values - The array of values
4638: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
4640: Fortran Notes:
4641: This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
4643: Level: intermediate
4645: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
4646: @*/
4647: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4648: {
4649: DM_Plex *mesh = (DM_Plex*) dm->data;
4650: PetscSection clSection;
4651: IS clPoints;
4652: PetscInt *points = NULL, *newPoints;
4653: const PetscInt *clp;
4654: PetscInt *indices;
4655: PetscInt offsets[32];
4656: PetscInt numFields, numPoints, newNumPoints, numIndices, newNumIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
4657: PetscScalar *newValues;
4658: PetscErrorCode ierr;
4662: if (!section) {DMGetDefaultSection(dm, §ion);}
4664: if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
4667: PetscSectionGetNumFields(section, &numFields);
4668: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4669: PetscMemzero(offsets, 32 * sizeof(PetscInt));
4670: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4671: if (!clPoints) {
4672: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4673: /* Compress out points not in the section */
4674: PetscSectionGetChart(section, &pStart, &pEnd);
4675: for (p = 0, q = 0; p < numPoints*2; p += 2) {
4676: if ((points[p] >= pStart) && (points[p] < pEnd)) {
4677: points[q*2] = points[p];
4678: points[q*2+1] = points[p+1];
4679: ++q;
4680: }
4681: }
4682: numPoints = q;
4683: } else {
4684: PetscInt dof, off;
4686: PetscSectionGetDof(clSection, point, &dof);
4687: numPoints = dof/2;
4688: PetscSectionGetOffset(clSection, point, &off);
4689: ISGetIndices(clPoints, &clp);
4690: points = (PetscInt *) &clp[off];
4691: }
4692: for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
4693: PetscInt fdof;
4695: PetscSectionGetDof(section, points[p], &dof);
4696: for (f = 0; f < numFields; ++f) {
4697: PetscSectionGetFieldDof(section, points[p], f, &fdof);
4698: offsets[f+1] += fdof;
4699: }
4700: numIndices += dof;
4701: }
4702: for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];
4704: if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
4705: DMPlexAnchorsModifyMat(dm,section,numPoints,numIndices,points,values,&newNumPoints,&newNumIndices,&newPoints,&newValues,offsets,PETSC_TRUE);
4706: if (newNumPoints) {
4707: if (!clPoints) {
4708: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4709: } else {
4710: ISRestoreIndices(clPoints, &clp);
4711: }
4712: numPoints = newNumPoints;
4713: numIndices = newNumIndices;
4714: points = newPoints;
4715: values = newValues;
4716: }
4717: DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
4718: if (numFields) {
4719: for (p = 0; p < numPoints*2; p += 2) {
4720: PetscInt o = points[p+1];
4721: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4722: indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
4723: }
4724: } else {
4725: for (p = 0, off = 0; p < numPoints*2; p += 2) {
4726: PetscInt o = points[p+1];
4727: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4728: indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
4729: }
4730: }
4731: if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
4732: MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
4733: if (mesh->printFEM > 1) {
4734: PetscInt i;
4735: PetscPrintf(PETSC_COMM_SELF, " Indices:");
4736: for (i = 0; i < numIndices; ++i) {PetscPrintf(PETSC_COMM_SELF, " %d", indices[i]);}
4737: PetscPrintf(PETSC_COMM_SELF, "\n");
4738: }
4739: if (ierr) {
4740: PetscMPIInt rank;
4741: PetscErrorCode ierr2;
4743: ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4744: ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4745: ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
4746: ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
4747:
4748: }
4749: if (newNumPoints) {
4750: DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4751: DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4752: }
4753: else {
4754: if (!clPoints) {
4755: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4756: } else {
4757: ISRestoreIndices(clPoints, &clp);
4758: }
4759: }
4760: DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
4761: return(0);
4762: }
4766: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4767: {
4768: DM_Plex *mesh = (DM_Plex*) dmf->data;
4769: PetscInt *fpoints = NULL, *ftotpoints = NULL;
4770: PetscInt *cpoints = NULL;
4771: PetscInt *findices, *cindices;
4772: PetscInt foffsets[32], coffsets[32];
4773: CellRefiner cellRefiner;
4774: PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4775: PetscErrorCode ierr;
4780: if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4782: if (!csection) {DMGetDefaultSection(dmc, &csection);}
4784: if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4786: if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4789: PetscSectionGetNumFields(fsection, &numFields);
4790: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4791: PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4792: PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4793: /* Column indices */
4794: DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4795: maxFPoints = numCPoints;
4796: /* Compress out points not in the section */
4797: /* TODO: Squeeze out points with 0 dof as well */
4798: PetscSectionGetChart(csection, &pStart, &pEnd);
4799: for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4800: if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4801: cpoints[q*2] = cpoints[p];
4802: cpoints[q*2+1] = cpoints[p+1];
4803: ++q;
4804: }
4805: }
4806: numCPoints = q;
4807: for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4808: PetscInt fdof;
4810: PetscSectionGetDof(csection, cpoints[p], &dof);
4811: if (!dof) continue;
4812: for (f = 0; f < numFields; ++f) {
4813: PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4814: coffsets[f+1] += fdof;
4815: }
4816: numCIndices += dof;
4817: }
4818: for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4819: /* Row indices */
4820: DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4821: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4822: DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4823: for (r = 0, q = 0; r < numSubcells; ++r) {
4824: /* TODO Map from coarse to fine cells */
4825: DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4826: /* Compress out points not in the section */
4827: PetscSectionGetChart(fsection, &pStart, &pEnd);
4828: for (p = 0; p < numFPoints*2; p += 2) {
4829: if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4830: PetscSectionGetDof(fsection, fpoints[p], &dof);
4831: if (!dof) continue;
4832: for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4833: if (s < q) continue;
4834: ftotpoints[q*2] = fpoints[p];
4835: ftotpoints[q*2+1] = fpoints[p+1];
4836: ++q;
4837: }
4838: }
4839: DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4840: }
4841: numFPoints = q;
4842: for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4843: PetscInt fdof;
4845: PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4846: if (!dof) continue;
4847: for (f = 0; f < numFields; ++f) {
4848: PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4849: foffsets[f+1] += fdof;
4850: }
4851: numFIndices += dof;
4852: }
4853: for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];
4855: if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4856: if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4857: DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4858: DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4859: if (numFields) {
4860: for (p = 0; p < numFPoints*2; p += 2) {
4861: PetscInt o = ftotpoints[p+1];
4862: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4863: indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4864: }
4865: for (p = 0; p < numCPoints*2; p += 2) {
4866: PetscInt o = cpoints[p+1];
4867: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4868: indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4869: }
4870: } else {
4871: for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4872: PetscInt o = ftotpoints[p+1];
4873: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4874: indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4875: }
4876: for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4877: PetscInt o = cpoints[p+1];
4878: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4879: indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4880: }
4881: }
4882: if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
4883: MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
4884: if (ierr) {
4885: PetscMPIInt rank;
4886: PetscErrorCode ierr2;
4888: ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4889: ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4890: ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
4891: ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
4892: ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
4893:
4894: }
4895: DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4896: DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4897: DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4898: DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4899: return(0);
4900: }
4904: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
4905: {
4906: PetscInt *fpoints = NULL, *ftotpoints = NULL;
4907: PetscInt *cpoints = NULL;
4908: PetscInt foffsets[32], coffsets[32];
4909: CellRefiner cellRefiner;
4910: PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4916: if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4918: if (!csection) {DMGetDefaultSection(dmc, &csection);}
4920: if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4922: if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4924: PetscSectionGetNumFields(fsection, &numFields);
4925: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4926: PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4927: PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4928: /* Column indices */
4929: DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4930: maxFPoints = numCPoints;
4931: /* Compress out points not in the section */
4932: /* TODO: Squeeze out points with 0 dof as well */
4933: PetscSectionGetChart(csection, &pStart, &pEnd);
4934: for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4935: if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4936: cpoints[q*2] = cpoints[p];
4937: cpoints[q*2+1] = cpoints[p+1];
4938: ++q;
4939: }
4940: }
4941: numCPoints = q;
4942: for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4943: PetscInt fdof;
4945: PetscSectionGetDof(csection, cpoints[p], &dof);
4946: if (!dof) continue;
4947: for (f = 0; f < numFields; ++f) {
4948: PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4949: coffsets[f+1] += fdof;
4950: }
4951: numCIndices += dof;
4952: }
4953: for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4954: /* Row indices */
4955: DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4956: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4957: DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4958: for (r = 0, q = 0; r < numSubcells; ++r) {
4959: /* TODO Map from coarse to fine cells */
4960: DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4961: /* Compress out points not in the section */
4962: PetscSectionGetChart(fsection, &pStart, &pEnd);
4963: for (p = 0; p < numFPoints*2; p += 2) {
4964: if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4965: PetscSectionGetDof(fsection, fpoints[p], &dof);
4966: if (!dof) continue;
4967: for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4968: if (s < q) continue;
4969: ftotpoints[q*2] = fpoints[p];
4970: ftotpoints[q*2+1] = fpoints[p+1];
4971: ++q;
4972: }
4973: }
4974: DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4975: }
4976: numFPoints = q;
4977: for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4978: PetscInt fdof;
4980: PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4981: if (!dof) continue;
4982: for (f = 0; f < numFields; ++f) {
4983: PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4984: foffsets[f+1] += fdof;
4985: }
4986: numFIndices += dof;
4987: }
4988: for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];
4990: if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4991: if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4992: if (numFields) {
4993: for (p = 0; p < numFPoints*2; p += 2) {
4994: PetscInt o = ftotpoints[p+1];
4995: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4996: indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4997: }
4998: for (p = 0; p < numCPoints*2; p += 2) {
4999: PetscInt o = cpoints[p+1];
5000: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5001: indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5002: }
5003: } else {
5004: for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5005: PetscInt o = ftotpoints[p+1];
5006: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5007: indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5008: }
5009: for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5010: PetscInt o = cpoints[p+1];
5011: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5012: indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5013: }
5014: }
5015: DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5016: DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5017: return(0);
5018: }
5022: /*@
5023: DMPlexGetHybridBounds - Get the first mesh point of each dimension which is a hybrid
5025: Input Parameter:
5026: . dm - The DMPlex object
5028: Output Parameters:
5029: + cMax - The first hybrid cell
5030: . fMax - The first hybrid face
5031: . eMax - The first hybrid edge
5032: - vMax - The first hybrid vertex
5034: Level: developer
5036: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
5037: @*/
5038: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
5039: {
5040: DM_Plex *mesh = (DM_Plex*) dm->data;
5041: PetscInt dim;
5046: DMGetDimension(dm, &dim);
5047: if (cMax) *cMax = mesh->hybridPointMax[dim];
5048: if (fMax) *fMax = mesh->hybridPointMax[dim-1];
5049: if (eMax) *eMax = mesh->hybridPointMax[1];
5050: if (vMax) *vMax = mesh->hybridPointMax[0];
5051: return(0);
5052: }
5056: /*@
5057: DMPlexSetHybridBounds - Set the first mesh point of each dimension which is a hybrid
5059: Input Parameters:
5060: . dm - The DMPlex object
5061: . cMax - The first hybrid cell
5062: . fMax - The first hybrid face
5063: . eMax - The first hybrid edge
5064: - vMax - The first hybrid vertex
5066: Level: developer
5068: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5069: @*/
5070: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5071: {
5072: DM_Plex *mesh = (DM_Plex*) dm->data;
5073: PetscInt dim;
5078: DMGetDimension(dm, &dim);
5079: if (cMax >= 0) mesh->hybridPointMax[dim] = cMax;
5080: if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5081: if (eMax >= 0) mesh->hybridPointMax[1] = eMax;
5082: if (vMax >= 0) mesh->hybridPointMax[0] = vMax;
5083: return(0);
5084: }
5088: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5089: {
5090: DM_Plex *mesh = (DM_Plex*) dm->data;
5095: *cellHeight = mesh->vtkCellHeight;
5096: return(0);
5097: }
5101: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5102: {
5103: DM_Plex *mesh = (DM_Plex*) dm->data;
5107: mesh->vtkCellHeight = cellHeight;
5108: return(0);
5109: }
5113: /* We can easily have a form that takes an IS instead */
5114: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5115: {
5116: PetscSection section, globalSection;
5117: PetscInt *numbers, p;
5121: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
5122: PetscSectionSetChart(section, pStart, pEnd);
5123: for (p = pStart; p < pEnd; ++p) {
5124: PetscSectionSetDof(section, p, 1);
5125: }
5126: PetscSectionSetUp(section);
5127: PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
5128: PetscMalloc1(pEnd - pStart, &numbers);
5129: for (p = pStart; p < pEnd; ++p) {
5130: PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5131: if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5132: else numbers[p-pStart] += shift;
5133: }
5134: ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5135: if (globalSize) {
5136: PetscLayout layout;
5137: PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5138: PetscLayoutGetSize(layout, globalSize);
5139: PetscLayoutDestroy(&layout);
5140: }
5141: PetscSectionDestroy(§ion);
5142: PetscSectionDestroy(&globalSection);
5143: return(0);
5144: }
5148: PetscErrorCode DMPlexCreateCellNumbering_Internal(DM dm, PetscBool includeHybrid, IS *globalCellNumbers)
5149: {
5150: PetscInt cellHeight, cStart, cEnd, cMax;
5154: DMPlexGetVTKCellHeight(dm, &cellHeight);
5155: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5156: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5157: if (cMax >= 0 && !includeHybrid) cEnd = PetscMin(cEnd, cMax);
5158: DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, globalCellNumbers);
5159: return(0);
5160: }
5164: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5165: {
5166: DM_Plex *mesh = (DM_Plex*) dm->data;
5171: if (!mesh->globalCellNumbers) {DMPlexCreateCellNumbering_Internal(dm, PETSC_FALSE, &mesh->globalCellNumbers);}
5172: *globalCellNumbers = mesh->globalCellNumbers;
5173: return(0);
5174: }
5178: PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM dm, PetscBool includeHybrid, IS *globalVertexNumbers)
5179: {
5180: PetscInt vStart, vEnd, vMax;
5185: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5186: DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5187: if (vMax >= 0 && !includeHybrid) vEnd = PetscMin(vEnd, vMax);
5188: DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, globalVertexNumbers);
5189: return(0);
5190: }
5194: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5195: {
5196: DM_Plex *mesh = (DM_Plex*) dm->data;
5201: if (!mesh->globalVertexNumbers) {DMPlexCreateVertexNumbering_Internal(dm, PETSC_FALSE, &mesh->globalVertexNumbers);}
5202: *globalVertexNumbers = mesh->globalVertexNumbers;
5203: return(0);
5204: }
5208: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5209: {
5210: IS nums[4];
5211: PetscInt depths[4];
5212: PetscInt depth, d, shift = 0;
5217: DMPlexGetDepth(dm, &depth);
5218: /* For unstratified meshes use dim instead of depth */
5219: if (depth < 0) {DMGetDimension(dm, &depth);}
5220: depths[0] = depth; depths[1] = 0;
5221: for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5222: for (d = 0; d <= depth; ++d) {
5223: PetscInt pStart, pEnd, gsize;
5225: DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5226: DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5227: shift += gsize;
5228: }
5229: ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5230: for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5231: return(0);
5232: }
5236: /*@
5237: DMPlexCheckSymmetry - Check that the adjacency information in the mesh is symmetric.
5239: Input Parameters:
5240: + dm - The DMPlex object
5242: Note: This is a useful diagnostic when creating meshes programmatically.
5244: Level: developer
5246: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5247: @*/
5248: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5249: {
5250: PetscSection coneSection, supportSection;
5251: const PetscInt *cone, *support;
5252: PetscInt coneSize, c, supportSize, s;
5253: PetscInt pStart, pEnd, p, csize, ssize;
5254: PetscErrorCode ierr;
5258: DMPlexGetConeSection(dm, &coneSection);
5259: DMPlexGetSupportSection(dm, &supportSection);
5260: /* Check that point p is found in the support of its cone points, and vice versa */
5261: DMPlexGetChart(dm, &pStart, &pEnd);
5262: for (p = pStart; p < pEnd; ++p) {
5263: DMPlexGetConeSize(dm, p, &coneSize);
5264: DMPlexGetCone(dm, p, &cone);
5265: for (c = 0; c < coneSize; ++c) {
5266: PetscBool dup = PETSC_FALSE;
5267: PetscInt d;
5268: for (d = c-1; d >= 0; --d) {
5269: if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5270: }
5271: DMPlexGetSupportSize(dm, cone[c], &supportSize);
5272: DMPlexGetSupport(dm, cone[c], &support);
5273: for (s = 0; s < supportSize; ++s) {
5274: if (support[s] == p) break;
5275: }
5276: if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5277: PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5278: for (s = 0; s < coneSize; ++s) {
5279: PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
5280: }
5281: PetscPrintf(PETSC_COMM_SELF, "\n");
5282: PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
5283: for (s = 0; s < supportSize; ++s) {
5284: PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
5285: }
5286: PetscPrintf(PETSC_COMM_SELF, "\n");
5287: if (dup) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
5288: else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
5289: }
5290: }
5291: DMPlexGetSupportSize(dm, p, &supportSize);
5292: DMPlexGetSupport(dm, p, &support);
5293: for (s = 0; s < supportSize; ++s) {
5294: DMPlexGetConeSize(dm, support[s], &coneSize);
5295: DMPlexGetCone(dm, support[s], &cone);
5296: for (c = 0; c < coneSize; ++c) {
5297: if (cone[c] == p) break;
5298: }
5299: if (c >= coneSize) {
5300: PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
5301: for (c = 0; c < supportSize; ++c) {
5302: PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
5303: }
5304: PetscPrintf(PETSC_COMM_SELF, "\n");
5305: PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
5306: for (c = 0; c < coneSize; ++c) {
5307: PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
5308: }
5309: PetscPrintf(PETSC_COMM_SELF, "\n");
5310: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
5311: }
5312: }
5313: }
5314: PetscSectionGetStorageSize(coneSection, &csize);
5315: PetscSectionGetStorageSize(supportSection, &ssize);
5316: if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
5317: return(0);
5318: }
5322: /*@
5323: DMPlexCheckSkeleton - Check that each cell has the correct number of vertices
5325: Input Parameters:
5326: + dm - The DMPlex object
5327: . isSimplex - Are the cells simplices or tensor products
5328: - cellHeight - Normally 0
5330: Note: This is a useful diagnostic when creating meshes programmatically.
5332: Level: developer
5334: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
5335: @*/
5336: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5337: {
5338: PetscInt dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;
5343: DMGetDimension(dm, &dim);
5344: switch (dim) {
5345: case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
5346: case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
5347: case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
5348: default:
5349: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
5350: }
5351: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5352: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5353: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5354: cMax = cMax >= 0 ? cMax : cEnd;
5355: for (c = cStart; c < cMax; ++c) {
5356: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
5358: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5359: for (cl = 0; cl < closureSize*2; cl += 2) {
5360: const PetscInt p = closure[cl];
5361: if ((p >= vStart) && (p < vEnd)) ++coneSize;
5362: }
5363: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5364: if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d vertices != %d", c, coneSize, numCorners);
5365: }
5366: for (c = cMax; c < cEnd; ++c) {
5367: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
5369: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5370: for (cl = 0; cl < closureSize*2; cl += 2) {
5371: const PetscInt p = closure[cl];
5372: if ((p >= vStart) && (p < vEnd)) ++coneSize;
5373: }
5374: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5375: if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has %d vertices > %d", c, coneSize, numHybridCorners);
5376: }
5377: return(0);
5378: }
5382: /*@
5383: DMPlexCheckFaces - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type
5385: Input Parameters:
5386: + dm - The DMPlex object
5387: . isSimplex - Are the cells simplices or tensor products
5388: - cellHeight - Normally 0
5390: Note: This is a useful diagnostic when creating meshes programmatically.
5392: Level: developer
5394: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
5395: @*/
5396: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5397: {
5398: PetscInt pMax[4];
5399: PetscInt dim, vStart, vEnd, cStart, cEnd, c, h;
5404: DMGetDimension(dm, &dim);
5405: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5406: DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
5407: for (h = cellHeight; h < dim; ++h) {
5408: DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
5409: for (c = cStart; c < cEnd; ++c) {
5410: const PetscInt *cone, *ornt, *faces;
5411: PetscInt numFaces, faceSize, coneSize,f;
5412: PetscInt *closure = NULL, closureSize, cl, numCorners = 0;
5414: if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
5415: DMPlexGetConeSize(dm, c, &coneSize);
5416: DMPlexGetCone(dm, c, &cone);
5417: DMPlexGetConeOrientation(dm, c, &ornt);
5418: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5419: for (cl = 0; cl < closureSize*2; cl += 2) {
5420: const PetscInt p = closure[cl];
5421: if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
5422: }
5423: DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
5424: if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
5425: for (f = 0; f < numFaces; ++f) {
5426: PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;
5428: DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
5429: for (cl = 0; cl < fclosureSize*2; cl += 2) {
5430: const PetscInt p = fclosure[cl];
5431: if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
5432: }
5433: 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);
5434: for (v = 0; v < fnumCorners; ++v) {
5435: 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]);
5436: }
5437: DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
5438: }
5439: DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
5440: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5441: }
5442: }
5443: return(0);
5444: }
5448: /* Pointwise interpolation
5449: Just code FEM for now
5450: u^f = I u^c
5451: sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
5452: u^f_i = sum_j psi^f_i I phi^c_j u^c_j
5453: I_{ij} = psi^f_i phi^c_j
5454: */
5455: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
5456: {
5457: PetscSection gsc, gsf;
5458: PetscInt m, n;
5459: void *ctx;
5460: DM cdm;
5461: PetscBool regular;
5465: DMGetDefaultGlobalSection(dmFine, &gsf);
5466: PetscSectionGetConstrainedStorageSize(gsf, &m);
5467: DMGetDefaultGlobalSection(dmCoarse, &gsc);
5468: PetscSectionGetConstrainedStorageSize(gsc, &n);
5470: MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
5471: MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
5472: MatSetType(*interpolation, dmCoarse->mattype);
5473: DMGetApplicationContext(dmFine, &ctx);
5475: DMGetCoarseDM(dmFine, &cdm);
5476: DMPlexGetRegularRefinement(dmFine, ®ular);
5477: if (regular && cdm == dmCoarse) {DMPlexComputeInterpolatorNested(dmCoarse, dmFine, *interpolation, ctx);}
5478: else {DMPlexComputeInterpolatorGeneral(dmCoarse, dmFine, *interpolation, ctx);}
5479: MatViewFromOptions(*interpolation, NULL, "-interp_mat_view");
5480: /* Use naive scaling */
5481: DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
5482: return(0);
5483: }
5487: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat)
5488: {
5490: VecScatter ctx;
5493: DMPlexComputeInjectorFEM(dmCoarse, dmFine, &ctx, NULL);
5494: MatCreateScatter(PetscObjectComm((PetscObject)ctx), ctx, mat);
5495: VecScatterDestroy(&ctx);
5496: return(0);
5497: }
5501: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
5502: {
5503: PetscSection section;
5504: IS *bcPoints, *bcComps;
5505: PetscBool *isFE;
5506: PetscInt *bcFields, *numComp, *numDof;
5507: PetscInt depth, dim, numBd, numBC = 0, numFields, bd, bc = 0, f;
5508: PetscInt cStart, cEnd, cEndInterior;
5512: DMGetNumFields(dm, &numFields);
5513: if (!numFields) return(0);
5514: /* FE and FV boundary conditions are handled slightly differently */
5515: PetscMalloc1(numFields, &isFE);
5516: for (f = 0; f < numFields; ++f) {
5517: PetscObject obj;
5518: PetscClassId id;
5520: DMGetField(dm, f, &obj);
5521: PetscObjectGetClassId(obj, &id);
5522: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
5523: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
5524: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
5525: }
5526: /* Allocate boundary point storage for FEM boundaries */
5527: DMPlexGetDepth(dm, &depth);
5528: DMGetDimension(dm, &dim);
5529: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
5530: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
5531: DMGetNumBoundary(dm, &numBd);
5532: for (bd = 0; bd < numBd; ++bd) {
5533: PetscInt field;
5534: PetscBool isEssential;
5536: DMGetBoundary(dm, bd, &isEssential, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL, NULL);
5537: if (isFE[field] && isEssential) ++numBC;
5538: }
5539: /* Add ghost cell boundaries for FVM */
5540: for (f = 0; f < numFields; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
5541: PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);
5542: /* Constrain ghost cells for FV */
5543: for (f = 0; f < numFields; ++f) {
5544: PetscInt *newidx, c;
5546: if (isFE[f] || cEndInterior < 0) continue;
5547: PetscMalloc1(cEnd-cEndInterior,&newidx);
5548: for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
5549: bcFields[bc] = f;
5550: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5551: }
5552: /* Handle FEM Dirichlet boundaries */
5553: for (bd = 0; bd < numBd; ++bd) {
5554: const char *bdLabel;
5555: DMLabel label;
5556: const PetscInt *comps;
5557: const PetscInt *values;
5558: PetscInt bd2, field, numComps, numValues;
5559: PetscBool isEssential, duplicate = PETSC_FALSE;
5561: DMGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
5562: if (!isFE[field]) continue;
5563: DMGetLabel(dm, bdLabel, &label);
5564: /* Only want to modify label once */
5565: for (bd2 = 0; bd2 < bd; ++bd2) {
5566: const char *bdname;
5567: DMGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
5568: PetscStrcmp(bdname, bdLabel, &duplicate);
5569: if (duplicate) break;
5570: }
5571: if (!duplicate && (isFE[field])) {
5572: /* don't complete cells, which are just present to give orientation to the boundary */
5573: DMPlexLabelComplete(dm, label);
5574: }
5575: /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
5576: if (isEssential) {
5577: PetscInt *newidx;
5578: PetscInt n, newn = 0, p, v;
5580: bcFields[bc] = field;
5581: if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
5582: for (v = 0; v < numValues; ++v) {
5583: IS tmp;
5584: const PetscInt *idx;
5586: DMGetStratumIS(dm, bdLabel, values[v], &tmp);
5587: if (!tmp) continue;
5588: ISGetLocalSize(tmp, &n);
5589: ISGetIndices(tmp, &idx);
5590: if (isFE[field]) {
5591: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
5592: } else {
5593: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
5594: }
5595: ISRestoreIndices(tmp, &idx);
5596: ISDestroy(&tmp);
5597: }
5598: PetscMalloc1(newn,&newidx);
5599: newn = 0;
5600: for (v = 0; v < numValues; ++v) {
5601: IS tmp;
5602: const PetscInt *idx;
5604: DMGetStratumIS(dm, bdLabel, values[v], &tmp);
5605: if (!tmp) continue;
5606: ISGetLocalSize(tmp, &n);
5607: ISGetIndices(tmp, &idx);
5608: if (isFE[field]) {
5609: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
5610: } else {
5611: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
5612: }
5613: ISRestoreIndices(tmp, &idx);
5614: ISDestroy(&tmp);
5615: }
5616: ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5617: }
5618: }
5619: /* Handle discretization */
5620: PetscCalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
5621: for (f = 0; f < numFields; ++f) {
5622: PetscObject obj;
5624: DMGetField(dm, f, &obj);
5625: if (isFE[f]) {
5626: PetscFE fe = (PetscFE) obj;
5627: const PetscInt *numFieldDof;
5628: PetscInt d;
5630: PetscFEGetNumComponents(fe, &numComp[f]);
5631: PetscFEGetNumDof(fe, &numFieldDof);
5632: for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
5633: } else {
5634: PetscFV fv = (PetscFV) obj;
5636: PetscFVGetNumComponents(fv, &numComp[f]);
5637: numDof[f*(dim+1)+dim] = numComp[f];
5638: }
5639: }
5640: for (f = 0; f < numFields; ++f) {
5641: PetscInt d;
5642: for (d = 1; d < dim; ++d) {
5643: 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.");
5644: }
5645: }
5646: DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, §ion);
5647: for (f = 0; f < numFields; ++f) {
5648: PetscFE fe;
5649: const char *name;
5651: DMGetField(dm, f, (PetscObject *) &fe);
5652: PetscObjectGetName((PetscObject) fe, &name);
5653: PetscSectionSetFieldName(section, f, name);
5654: }
5655: DMSetDefaultSection(dm, section);
5656: PetscSectionDestroy(§ion);
5657: for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
5658: PetscFree3(bcFields,bcPoints,bcComps);
5659: PetscFree2(numComp,numDof);
5660: PetscFree(isFE);
5661: return(0);
5662: }
5666: /*@
5667: DMPlexGetRegularRefinement - Get the flag indicating that this mesh was obtained by regular refinement from its coarse mesh
5669: Input Parameter:
5670: . dm - The DMPlex object
5672: Output Parameter:
5673: . regular - The flag
5675: Level: intermediate
5677: .seealso: DMPlexSetRegularRefinement()
5678: @*/
5679: PetscErrorCode DMPlexGetRegularRefinement(DM dm, PetscBool *regular)
5680: {
5684: *regular = ((DM_Plex *) dm->data)->regularRefinement;
5685: return(0);
5686: }
5690: /*@
5691: DMPlexSetRegularRefinement - Set the flag indicating that this mesh was obtained by regular refinement from its coarse mesh
5693: Input Parameters:
5694: + dm - The DMPlex object
5695: - regular - The flag
5697: Level: intermediate
5699: .seealso: DMPlexGetRegularRefinement()
5700: @*/
5701: PetscErrorCode DMPlexSetRegularRefinement(DM dm, PetscBool regular)
5702: {
5705: ((DM_Plex *) dm->data)->regularRefinement = regular;
5706: return(0);
5707: }
5709: /* anchors */
5712: /*@
5713: DMPlexGetAnchors - Get the layout of the anchor (point-to-point) constraints. Typically, the user will not have to
5714: call DMPlexGetAnchors() directly: if there are anchors, then DMPlexGetAnchors() is called during DMGetConstraints().
5716: not collective
5718: Input Parameters:
5719: . dm - The DMPlex object
5721: Output Parameters:
5722: + anchorSection - If not NULL, set to the section describing which points anchor the constrained points.
5723: - anchorIS - If not NULL, set to the list of anchors indexed by anchorSection
5726: Level: intermediate
5728: .seealso: DMPlexSetAnchors(), DMGetConstraints(), DMSetConstraints()
5729: @*/
5730: PetscErrorCode DMPlexGetAnchors(DM dm, PetscSection *anchorSection, IS *anchorIS)
5731: {
5732: DM_Plex *plex = (DM_Plex *)dm->data;
5737: if (!plex->anchorSection && !plex->anchorIS && plex->createanchors) {(*plex->createanchors)(dm);}
5738: if (anchorSection) *anchorSection = plex->anchorSection;
5739: if (anchorIS) *anchorIS = plex->anchorIS;
5740: return(0);
5741: }
5745: /*@
5746: DMPlexSetAnchors - Set the layout of the local anchor (point-to-point) constraints. Unlike boundary conditions,
5747: when a point's degrees of freedom in a section are constrained to an outside value, the anchor constraints set a
5748: point's degrees of freedom to be a linear combination of other points' degrees of freedom.
5750: After specifying the layout of constraints with DMPlexSetAnchors(), one specifies the constraints by calling
5751: DMGetConstraints() and filling in the entries in the constraint matrix.
5753: collective on dm
5755: Input Parameters:
5756: + dm - The DMPlex object
5757: . 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).
5758: - anchorIS - The list of all anchor points. Must have a local communicator (PETSC_COMM_SELF or derivative).
5760: The reference counts of anchorSection and anchorIS are incremented.
5762: Level: intermediate
5764: .seealso: DMPlexGetAnchors(), DMGetConstraints(), DMSetConstraints()
5765: @*/
5766: PetscErrorCode DMPlexSetAnchors(DM dm, PetscSection anchorSection, IS anchorIS)
5767: {
5768: DM_Plex *plex = (DM_Plex *)dm->data;
5769: PetscMPIInt result;
5774: if (anchorSection) {
5776: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorSection),&result);
5777: if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor section must have local communicator");
5778: }
5779: if (anchorIS) {
5781: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorIS),&result);
5782: if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor IS must have local communicator");
5783: }
5785: PetscObjectReference((PetscObject)anchorSection);
5786: PetscSectionDestroy(&plex->anchorSection);
5787: plex->anchorSection = anchorSection;
5789: PetscObjectReference((PetscObject)anchorIS);
5790: ISDestroy(&plex->anchorIS);
5791: plex->anchorIS = anchorIS;
5793: #if defined(PETSC_USE_DEBUG)
5794: if (anchorIS && anchorSection) {
5795: PetscInt size, a, pStart, pEnd;
5796: const PetscInt *anchors;
5798: PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5799: ISGetLocalSize(anchorIS,&size);
5800: ISGetIndices(anchorIS,&anchors);
5801: for (a = 0; a < size; a++) {
5802: PetscInt p;
5804: p = anchors[a];
5805: if (p >= pStart && p < pEnd) {
5806: PetscInt dof;
5808: PetscSectionGetDof(anchorSection,p,&dof);
5809: if (dof) {
5810: PetscErrorCode ierr2;
5812: ierr2 = ISRestoreIndices(anchorIS,&anchors);CHKERRQ(ierr2);
5813: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point %d cannot be constrained and an anchor",p);
5814: }
5815: }
5816: }
5817: ISRestoreIndices(anchorIS,&anchors);
5818: }
5819: #endif
5820: /* reset the generic constraints */
5821: DMSetDefaultConstraints(dm,NULL,NULL);
5822: return(0);
5823: }
5827: static PetscErrorCode DMPlexCreateConstraintSection_Anchors(DM dm, PetscSection section, PetscSection *cSec)
5828: {
5829: PetscSection anchorSection;
5830: PetscInt pStart, pEnd, sStart, sEnd, p, dof, numFields, f;
5835: DMPlexGetAnchors(dm,&anchorSection,NULL);
5836: PetscSectionCreate(PETSC_COMM_SELF,cSec);
5837: PetscSectionGetNumFields(section,&numFields);
5838: if (numFields) {
5839: PetscInt f;
5840: PetscSectionSetNumFields(*cSec,numFields);
5842: for (f = 0; f < numFields; f++) {
5843: PetscInt numComp;
5845: PetscSectionGetFieldComponents(section,f,&numComp);
5846: PetscSectionSetFieldComponents(*cSec,f,numComp);
5847: }
5848: }
5849: PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5850: PetscSectionGetChart(section,&sStart,&sEnd);
5851: pStart = PetscMax(pStart,sStart);
5852: pEnd = PetscMin(pEnd,sEnd);
5853: pEnd = PetscMax(pStart,pEnd);
5854: PetscSectionSetChart(*cSec,pStart,pEnd);
5855: for (p = pStart; p < pEnd; p++) {
5856: PetscSectionGetDof(anchorSection,p,&dof);
5857: if (dof) {
5858: PetscSectionGetDof(section,p,&dof);
5859: PetscSectionSetDof(*cSec,p,dof);
5860: for (f = 0; f < numFields; f++) {
5861: PetscSectionGetFieldDof(section,p,f,&dof);
5862: PetscSectionSetFieldDof(*cSec,p,f,dof);
5863: }
5864: }
5865: }
5866: PetscSectionSetUp(*cSec);
5867: return(0);
5868: }
5872: static PetscErrorCode DMPlexCreateConstraintMatrix_Anchors(DM dm, PetscSection section, PetscSection cSec, Mat *cMat)
5873: {
5874: PetscSection aSec;
5875: PetscInt pStart, pEnd, p, dof, aDof, aOff, off, nnz, annz, m, n, q, a, offset, *i, *j;
5876: const PetscInt *anchors;
5877: PetscInt numFields, f;
5878: IS aIS;
5883: PetscSectionGetStorageSize(cSec, &m);
5884: PetscSectionGetStorageSize(section, &n);
5885: MatCreate(PETSC_COMM_SELF,cMat);
5886: MatSetSizes(*cMat,m,n,m,n);
5887: MatSetType(*cMat,MATSEQAIJ);
5888: DMPlexGetAnchors(dm,&aSec,&aIS);
5889: ISGetIndices(aIS,&anchors);
5890: /* cSec will be a subset of aSec and section */
5891: PetscSectionGetChart(cSec,&pStart,&pEnd);
5892: PetscMalloc1(m+1,&i);
5893: i[0] = 0;
5894: PetscSectionGetNumFields(section,&numFields);
5895: for (p = pStart; p < pEnd; p++) {
5896: PetscInt rDof, rOff, r;
5898: PetscSectionGetDof(aSec,p,&rDof);
5899: if (!rDof) continue;
5900: PetscSectionGetOffset(aSec,p,&rOff);
5901: if (numFields) {
5902: for (f = 0; f < numFields; f++) {
5903: annz = 0;
5904: for (r = 0; r < rDof; r++) {
5905: a = anchors[rOff + r];
5906: PetscSectionGetFieldDof(section,a,f,&aDof);
5907: annz += aDof;
5908: }
5909: PetscSectionGetFieldDof(cSec,p,f,&dof);
5910: PetscSectionGetFieldOffset(cSec,p,f,&off);
5911: for (q = 0; q < dof; q++) {
5912: i[off + q + 1] = i[off + q] + annz;
5913: }
5914: }
5915: }
5916: else {
5917: annz = 0;
5918: for (q = 0; q < dof; q++) {
5919: a = anchors[off + q];
5920: PetscSectionGetDof(section,a,&aDof);
5921: annz += aDof;
5922: }
5923: PetscSectionGetDof(cSec,p,&dof);
5924: PetscSectionGetOffset(cSec,p,&off);
5925: for (q = 0; q < dof; q++) {
5926: i[off + q + 1] = i[off + q] + annz;
5927: }
5928: }
5929: }
5930: nnz = i[m];
5931: PetscMalloc1(nnz,&j);
5932: offset = 0;
5933: for (p = pStart; p < pEnd; p++) {
5934: if (numFields) {
5935: for (f = 0; f < numFields; f++) {
5936: PetscSectionGetFieldDof(cSec,p,f,&dof);
5937: for (q = 0; q < dof; q++) {
5938: PetscInt rDof, rOff, r;
5939: PetscSectionGetDof(aSec,p,&rDof);
5940: PetscSectionGetOffset(aSec,p,&rOff);
5941: for (r = 0; r < rDof; r++) {
5942: PetscInt s;
5944: a = anchors[rOff + r];
5945: PetscSectionGetFieldDof(section,a,f,&aDof);
5946: PetscSectionGetFieldOffset(section,a,f,&aOff);
5947: for (s = 0; s < aDof; s++) {
5948: j[offset++] = aOff + s;
5949: }
5950: }
5951: }
5952: }
5953: }
5954: else {
5955: PetscSectionGetDof(cSec,p,&dof);
5956: for (q = 0; q < dof; q++) {
5957: PetscInt rDof, rOff, r;
5958: PetscSectionGetDof(aSec,p,&rDof);
5959: PetscSectionGetOffset(aSec,p,&rOff);
5960: for (r = 0; r < rDof; r++) {
5961: PetscInt s;
5963: a = anchors[rOff + r];
5964: PetscSectionGetDof(section,a,&aDof);
5965: PetscSectionGetOffset(section,a,&aOff);
5966: for (s = 0; s < aDof; s++) {
5967: j[offset++] = aOff + s;
5968: }
5969: }
5970: }
5971: }
5972: }
5973: MatSeqAIJSetPreallocationCSR(*cMat,i,j,NULL);
5974: PetscFree(i);
5975: PetscFree(j);
5976: ISRestoreIndices(aIS,&anchors);
5977: return(0);
5978: }
5982: PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm)
5983: {
5984: DM_Plex *plex = (DM_Plex *)dm->data;
5985: PetscSection anchorSection, section, cSec;
5986: Mat cMat;
5991: DMPlexGetAnchors(dm,&anchorSection,NULL);
5992: if (anchorSection) {
5993: PetscDS ds;
5994: PetscInt nf;
5996: DMGetDefaultSection(dm,§ion);
5997: DMPlexCreateConstraintSection_Anchors(dm,section,&cSec);
5998: DMPlexCreateConstraintMatrix_Anchors(dm,section,cSec,&cMat);
5999: DMGetDS(dm,&ds);
6000: PetscDSGetNumFields(ds,&nf);
6001: if (nf && plex->computeanchormatrix) {(*plex->computeanchormatrix)(dm,section,cSec,cMat);}
6002: DMSetDefaultConstraints(dm,cSec,cMat);
6003: PetscSectionDestroy(&cSec);
6004: MatDestroy(&cMat);
6005: }
6006: return(0);
6007: }