Actual source code: plex.c
petsc-3.6.4 2016-04-12
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petsc/private/isimpl.h>
3: #include <petscsf.h>
4: #include <petscds.h>
6: /* Logging support */
7: PetscLogEvent DMPLEX_Interpolate, PETSCPARTITIONER_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeOverlap, DMPLEX_DistributeField, DMPLEX_DistributeData, DMPLEX_Migrate, DMPLEX_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM, DMPLEX_InterpolatorFEM, DMPLEX_InjectorFEM, DMPLEX_IntegralFEM, DMPLEX_CreateGmsh;
9: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
10: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
11: PETSC_EXTERN PetscErrorCode VecLoad_Default(Vec, PetscViewer);
15: PetscErrorCode DMPlexGetFieldType_Internal(DM dm, PetscSection section, PetscInt field, PetscInt *sStart, PetscInt *sEnd, PetscViewerVTKFieldType *ft)
16: {
17: PetscInt dim, pStart, pEnd, vStart, vEnd, cStart, cEnd, cEndInterior, vdof = 0, cdof = 0;
21: *ft = PETSC_VTK_POINT_FIELD;
22: DMGetDimension(dm, &dim);
23: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
24: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
25: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
26: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
27: PetscSectionGetChart(section, &pStart, &pEnd);
28: if (field >= 0) {
29: if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetFieldDof(section, vStart, field, &vdof);}
30: if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetFieldDof(section, cStart, field, &cdof);}
31: } else {
32: if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetDof(section, vStart, &vdof);}
33: if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetDof(section, cStart, &cdof);}
34: }
35: if (vdof) {
36: *sStart = vStart;
37: *sEnd = vEnd;
38: if (vdof == dim) *ft = PETSC_VTK_POINT_VECTOR_FIELD;
39: else *ft = PETSC_VTK_POINT_FIELD;
40: } else if (cdof) {
41: *sStart = cStart;
42: *sEnd = cEnd;
43: if (cdof == dim) *ft = PETSC_VTK_CELL_VECTOR_FIELD;
44: else *ft = PETSC_VTK_CELL_FIELD;
45: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Could not classify input Vec for VTK");
46: return(0);
47: }
51: PetscErrorCode VecView_Plex_Local(Vec v, PetscViewer viewer)
52: {
53: DM dm;
54: PetscBool isvtk, ishdf5, isseq;
58: VecGetDM(v, &dm);
59: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
60: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
61: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
62: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
63: if (isvtk || ishdf5) {
64: PetscInt numFields;
65: PetscBool fem = PETSC_FALSE;
67: DMGetNumFields(dm, &numFields);
68: if (numFields) {
69: PetscObject fe;
71: DMGetField(dm, 0, &fe);
72: if (fe->classid == PETSCFE_CLASSID) fem = PETSC_TRUE;
73: }
74: if (fem) {DMPlexInsertBoundaryValues(dm, v, 0.0, NULL, NULL, NULL);}
75: }
76: if (isvtk) {
77: PetscSection section;
78: PetscViewerVTKFieldType ft;
79: PetscInt pStart, pEnd;
81: DMGetDefaultSection(dm, §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 VecLoad_Plex_Local(Vec v, PetscViewer viewer)
140: {
141: DM dm;
142: PetscBool ishdf5;
146: VecGetDM(v, &dm);
147: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
148: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
149: if (ishdf5) {
150: DM dmBC;
151: Vec gv;
152: const char *name;
154: DMGetOutputDM(dm, &dmBC);
155: DMGetGlobalVector(dmBC, &gv);
156: PetscObjectGetName((PetscObject) v, &name);
157: PetscObjectSetName((PetscObject) gv, name);
158: VecLoad_Default(gv, viewer);
159: DMGlobalToLocalBegin(dmBC, gv, INSERT_VALUES, v);
160: DMGlobalToLocalEnd(dmBC, gv, INSERT_VALUES, v);
161: DMRestoreGlobalVector(dmBC, &gv);
162: } else {
163: VecLoad_Default(v, viewer);
164: }
165: return(0);
166: }
170: PetscErrorCode VecLoad_Plex(Vec v, PetscViewer viewer)
171: {
172: DM dm;
173: PetscBool ishdf5;
177: VecGetDM(v, &dm);
178: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
179: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
180: if (ishdf5) {
181: #if defined(PETSC_HAVE_HDF5)
182: VecLoad_Plex_HDF5(v, viewer);
183: #else
184: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
185: #endif
186: } else {
187: VecLoad_Default(v, viewer);
188: }
189: return(0);
190: }
194: PetscErrorCode DMPlexView_Ascii_Geometry(DM dm, PetscViewer viewer)
195: {
196: PetscSection coordSection;
197: Vec coordinates;
198: DMLabel depthLabel;
199: const char *name[4];
200: const PetscScalar *a;
201: PetscInt dim, pStart, pEnd, cStart, cEnd, c;
202: PetscErrorCode ierr;
205: DMGetDimension(dm, &dim);
206: DMGetCoordinatesLocal(dm, &coordinates);
207: DMGetCoordinateSection(dm, &coordSection);
208: DMPlexGetDepthLabel(dm, &depthLabel);
209: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
210: PetscSectionGetChart(coordSection, &pStart, &pEnd);
211: VecGetArrayRead(coordinates, &a);
212: name[0] = "vertex";
213: name[1] = "edge";
214: name[dim-1] = "face";
215: name[dim] = "cell";
216: for (c = cStart; c < cEnd; ++c) {
217: PetscInt *closure = NULL;
218: PetscInt closureSize, cl;
220: PetscViewerASCIIPrintf(viewer, "Geometry for cell %d:\n", c);
221: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
222: PetscViewerASCIIPushTab(viewer);
223: for (cl = 0; cl < closureSize*2; cl += 2) {
224: PetscInt point = closure[cl], depth, dof, off, d, p;
226: if ((point < pStart) || (point >= pEnd)) continue;
227: PetscSectionGetDof(coordSection, point, &dof);
228: if (!dof) continue;
229: DMLabelGetValue(depthLabel, point, &depth);
230: PetscSectionGetOffset(coordSection, point, &off);
231: PetscViewerASCIIPrintf(viewer, "%s %d coords:", name[depth], point);
232: for (p = 0; p < dof/dim; ++p) {
233: PetscViewerASCIIPrintf(viewer, " (");
234: for (d = 0; d < dim; ++d) {
235: if (d > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
236: PetscViewerASCIIPrintf(viewer, "%g", PetscRealPart(a[off+p*dim+d]));
237: }
238: PetscViewerASCIIPrintf(viewer, ")");
239: }
240: PetscViewerASCIIPrintf(viewer, "\n");
241: }
242: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
243: PetscViewerASCIIPopTab(viewer);
244: }
245: VecRestoreArrayRead(coordinates, &a);
246: return(0);
247: }
251: PetscErrorCode DMPlexView_Ascii(DM dm, PetscViewer viewer)
252: {
253: DM_Plex *mesh = (DM_Plex*) dm->data;
254: DM cdm;
255: DMLabel markers;
256: PetscSection coordSection;
257: Vec coordinates;
258: PetscViewerFormat format;
259: PetscErrorCode ierr;
262: DMGetCoordinateDM(dm, &cdm);
263: DMGetDefaultSection(cdm, &coordSection);
264: DMGetCoordinatesLocal(dm, &coordinates);
265: PetscViewerGetFormat(viewer, &format);
266: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
267: const char *name;
268: PetscInt maxConeSize, maxSupportSize;
269: PetscInt pStart, pEnd, p;
270: PetscMPIInt rank, size;
272: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
273: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
274: PetscObjectGetName((PetscObject) dm, &name);
275: DMPlexGetChart(dm, &pStart, &pEnd);
276: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
277: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
278: PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
279: PetscViewerASCIISynchronizedPrintf(viewer, "Max sizes cone: %D support: %D\n", maxConeSize, maxSupportSize);
280: PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
281: PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
282: for (p = pStart; p < pEnd; ++p) {
283: PetscInt dof, off, s;
285: PetscSectionGetDof(mesh->supportSection, p, &dof);
286: PetscSectionGetOffset(mesh->supportSection, p, &off);
287: for (s = off; s < off+dof; ++s) {
288: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D ----> %D\n", rank, p, mesh->supports[s]);
289: }
290: }
291: PetscViewerFlush(viewer);
292: PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
293: for (p = pStart; p < pEnd; ++p) {
294: PetscInt dof, off, c;
296: PetscSectionGetDof(mesh->coneSection, p, &dof);
297: PetscSectionGetOffset(mesh->coneSection, p, &off);
298: for (c = off; c < off+dof; ++c) {
299: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %D <---- %D (%D)\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]);
300: }
301: }
302: PetscViewerFlush(viewer);
303: PetscSectionGetChart(coordSection, &pStart, NULL);
304: if (pStart >= 0) {PetscSectionVecView(coordSection, coordinates, viewer);}
305: DMPlexGetLabel(dm, "marker", &markers);
306: DMLabelView(markers,viewer);
307: if (size > 1) {
308: PetscSF sf;
310: DMGetPointSF(dm, &sf);
311: PetscSFView(sf, viewer);
312: }
313: PetscViewerFlush(viewer);
314: } else if (format == PETSC_VIEWER_ASCII_LATEX) {
315: const char *name, *color;
316: const char *defcolors[3] = {"gray", "orange", "green"};
317: const char *deflcolors[4] = {"blue", "cyan", "red", "magenta"};
318: PetscReal scale = 2.0;
319: PetscBool useNumbers = PETSC_TRUE, useLabels, useColors;
320: double tcoords[3];
321: PetscScalar *coords;
322: PetscInt numLabels, l, numColors, numLColors, dim, depth, cStart, cEnd, c, vStart, vEnd, v, eStart = 0, eEnd = 0, e, p;
323: PetscMPIInt rank, size;
324: char **names, **colors, **lcolors;
326: DMGetDimension(dm, &dim);
327: DMPlexGetDepth(dm, &depth);
328: DMPlexGetNumLabels(dm, &numLabels);
329: numLabels = PetscMax(numLabels, 10);
330: numColors = 10;
331: numLColors = 10;
332: PetscCalloc3(numLabels, &names, numColors, &colors, numLColors, &lcolors);
333: PetscOptionsGetReal(((PetscObject) viewer)->prefix, "-dm_plex_view_scale", &scale, NULL);
334: PetscOptionsGetBool(((PetscObject) viewer)->prefix, "-dm_plex_view_numbers", &useNumbers, NULL);
335: PetscOptionsGetStringArray(((PetscObject) viewer)->prefix, "-dm_plex_view_labels", names, &numLabels, &useLabels);
336: if (!useLabels) numLabels = 0;
337: PetscOptionsGetStringArray(((PetscObject) viewer)->prefix, "-dm_plex_view_colors", colors, &numColors, &useColors);
338: if (!useColors) {
339: numColors = 3;
340: for (c = 0; c < numColors; ++c) {PetscStrallocpy(defcolors[c], &colors[c]);}
341: }
342: PetscOptionsGetStringArray(((PetscObject) viewer)->prefix, "-dm_plex_view_lcolors", lcolors, &numLColors, &useColors);
343: if (!useColors) {
344: numLColors = 4;
345: for (c = 0; c < numLColors; ++c) {PetscStrallocpy(deflcolors[c], &lcolors[c]);}
346: }
347: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
348: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
349: PetscObjectGetName((PetscObject) dm, &name);
350: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
351: PetscViewerASCIIPrintf(viewer, "\
352: \\documentclass[tikz]{standalone}\n\n\
353: \\usepackage{pgflibraryshapes}\n\
354: \\usetikzlibrary{backgrounds}\n\
355: \\usetikzlibrary{arrows}\n\
356: \\begin{document}\n");
357: if (size > 1) {
358: PetscViewerASCIIPrintf(viewer, "%s for process ", name);
359: for (p = 0; p < size; ++p) {
360: if (p > 0 && p == size-1) {
361: PetscViewerASCIIPrintf(viewer, ", and ", colors[p%numColors], p);
362: } else if (p > 0) {
363: PetscViewerASCIIPrintf(viewer, ", ", colors[p%numColors], p);
364: }
365: PetscViewerASCIIPrintf(viewer, "{\\textcolor{%s}%D}", colors[p%numColors], p);
366: }
367: PetscViewerASCIIPrintf(viewer, ".\n\n\n");
368: }
369: PetscViewerASCIIPrintf(viewer, "\\begin{tikzpicture}[scale = %g,font=\\fontsize{8}{8}\\selectfont]\n", 1.0);
370: /* Plot vertices */
371: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
372: VecGetArray(coordinates, &coords);
373: for (v = vStart; v < vEnd; ++v) {
374: PetscInt off, dof, d;
375: PetscBool isLabeled = PETSC_FALSE;
377: PetscSectionGetDof(coordSection, v, &dof);
378: PetscSectionGetOffset(coordSection, v, &off);
379: PetscViewerASCIISynchronizedPrintf(viewer, "\\path (");
380: if (PetscUnlikely(dof > 3)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"coordSection vertex %D has dof %D > 3",v,dof);
381: for (d = 0; d < dof; ++d) {
382: tcoords[d] = (double) (scale*PetscRealPart(coords[off+d]));
383: tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
384: }
385: /* Rotate coordinates since PGF makes z point out of the page instead of up */
386: if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
387: for (d = 0; d < dof; ++d) {
388: if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
389: PetscViewerASCIISynchronizedPrintf(viewer, "%g", tcoords[d]);
390: }
391: color = colors[rank%numColors];
392: for (l = 0; l < numLabels; ++l) {
393: PetscInt val;
394: DMPlexGetLabelValue(dm, names[l], v, &val);
395: if (val >= 0) {color = lcolors[l%numLColors]; isLabeled = PETSC_TRUE; break;}
396: }
397: if (useNumbers) {
398: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [draw,shape=circle,color=%s] {%D};\n", v, rank, color, v);
399: } else {
400: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [fill,inner sep=%dpt,shape=circle,color=%s] {};\n", v, rank, !isLabeled ? 1 : 2, color);
401: }
402: }
403: VecRestoreArray(coordinates, &coords);
404: PetscViewerFlush(viewer);
405: /* Plot edges */
406: if (depth > 1) {DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);}
407: if (dim < 3 && useNumbers) {
408: VecGetArray(coordinates, &coords);
409: PetscViewerASCIIPrintf(viewer, "\\path\n");
410: for (e = eStart; e < eEnd; ++e) {
411: const PetscInt *cone;
412: PetscInt coneSize, offA, offB, dof, d;
414: DMPlexGetConeSize(dm, e, &coneSize);
415: if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %d cone should have two vertices, not %d", e, coneSize);
416: DMPlexGetCone(dm, e, &cone);
417: PetscSectionGetDof(coordSection, cone[0], &dof);
418: PetscSectionGetOffset(coordSection, cone[0], &offA);
419: PetscSectionGetOffset(coordSection, cone[1], &offB);
420: PetscViewerASCIISynchronizedPrintf(viewer, "(");
421: for (d = 0; d < dof; ++d) {
422: tcoords[d] = (double) (scale*PetscRealPart(coords[offA+d]+coords[offB+d]));
423: tcoords[d] = PetscAbsReal(tcoords[d]) < 1e-10 ? 0.0 : tcoords[d];
424: }
425: /* Rotate coordinates since PGF makes z point out of the page instead of up */
426: if (dim == 3) {PetscReal tmp = tcoords[1]; tcoords[1] = tcoords[2]; tcoords[2] = -tmp;}
427: for (d = 0; d < dof; ++d) {
428: if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
429: PetscViewerASCIISynchronizedPrintf(viewer, "%g", tcoords[d]);
430: }
431: color = colors[rank%numColors];
432: for (l = 0; l < numLabels; ++l) {
433: PetscInt val;
434: DMPlexGetLabelValue(dm, names[l], v, &val);
435: if (val >= 0) {color = lcolors[l%numLColors]; break;}
436: }
437: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%d) [draw,shape=circle,color=%s] {%D} --\n", e, rank, color, e);
438: }
439: VecRestoreArray(coordinates, &coords);
440: PetscViewerFlush(viewer);
441: PetscViewerASCIIPrintf(viewer, "(0,0);\n");
442: }
443: /* Plot cells */
444: if (dim == 3 || !useNumbers) {
445: for (e = eStart; e < eEnd; ++e) {
446: const PetscInt *cone;
448: color = colors[rank%numColors];
449: for (l = 0; l < numLabels; ++l) {
450: PetscInt val;
451: DMPlexGetLabelValue(dm, names[l], e, &val);
452: if (val >= 0) {color = lcolors[l%numLColors]; break;}
453: }
454: DMPlexGetCone(dm, e, &cone);
455: PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] (%D_%d) -- (%D_%d);\n", color, cone[0], rank, cone[1], rank);
456: }
457: } else {
458: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
459: for (c = cStart; c < cEnd; ++c) {
460: PetscInt *closure = NULL;
461: PetscInt closureSize, firstPoint = -1;
463: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
464: PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] ", colors[rank%numColors]);
465: for (p = 0; p < closureSize*2; p += 2) {
466: const PetscInt point = closure[p];
468: if ((point < vStart) || (point >= vEnd)) continue;
469: if (firstPoint >= 0) {PetscViewerASCIISynchronizedPrintf(viewer, " -- ");}
470: PetscViewerASCIISynchronizedPrintf(viewer, "(%D_%d)", point, rank);
471: if (firstPoint < 0) firstPoint = point;
472: }
473: /* Why doesn't this work? PetscViewerASCIISynchronizedPrintf(viewer, " -- cycle;\n"); */
474: PetscViewerASCIISynchronizedPrintf(viewer, " -- (%D_%d);\n", firstPoint, rank);
475: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
476: }
477: }
478: PetscViewerFlush(viewer);
479: PetscViewerASCIIPrintf(viewer, "\\end{tikzpicture}\n");
480: PetscViewerASCIIPrintf(viewer, "\\end{document}\n", name);
481: for (l = 0; l < numLabels; ++l) {PetscFree(names[l]);}
482: for (c = 0; c < numColors; ++c) {PetscFree(colors[c]);}
483: for (c = 0; c < numLColors; ++c) {PetscFree(lcolors[c]);}
484: PetscFree3(names, colors, lcolors);
485: } else {
486: MPI_Comm comm;
487: PetscInt *sizes, *hybsizes;
488: PetscInt locDepth, depth, dim, d, pMax[4];
489: PetscInt pStart, pEnd, p;
490: PetscInt numLabels, l;
491: const char *name;
492: PetscMPIInt size;
494: PetscObjectGetComm((PetscObject)dm,&comm);
495: MPI_Comm_size(comm, &size);
496: DMGetDimension(dm, &dim);
497: PetscObjectGetName((PetscObject) dm, &name);
498: if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dim);}
499: else {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dim);}
500: DMPlexGetDepth(dm, &locDepth);
501: MPI_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
502: DMPlexGetHybridBounds(dm, &pMax[depth], depth > 0 ? &pMax[depth-1] : NULL, &pMax[1], &pMax[0]);
503: PetscMalloc2(size,&sizes,size,&hybsizes);
504: if (depth == 1) {
505: DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd);
506: pEnd = pEnd - pStart;
507: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
508: PetscViewerASCIIPrintf(viewer, " %D-cells:", 0);
509: for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
510: PetscViewerASCIIPrintf(viewer, "\n");
511: DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd);
512: pEnd = pEnd - pStart;
513: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
514: PetscViewerASCIIPrintf(viewer, " %D-cells:", dim);
515: for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
516: PetscViewerASCIIPrintf(viewer, "\n");
517: } else {
518: for (d = 0; d <= dim; d++) {
519: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
520: pEnd -= pStart;
521: pMax[d] -= pStart;
522: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
523: MPI_Gather(&pMax[d], 1, MPIU_INT, hybsizes, 1, MPIU_INT, 0, comm);
524: PetscViewerASCIIPrintf(viewer, " %D-cells:", d);
525: for (p = 0; p < size; ++p) {
526: if (hybsizes[p] >= 0) {PetscViewerASCIIPrintf(viewer, " %D (%D)", sizes[p], sizes[p] - hybsizes[p]);}
527: else {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
528: }
529: PetscViewerASCIIPrintf(viewer, "\n");
530: }
531: }
532: PetscFree2(sizes,hybsizes);
533: DMPlexGetNumLabels(dm, &numLabels);
534: if (numLabels) {PetscViewerASCIIPrintf(viewer, "Labels:\n");}
535: for (l = 0; l < numLabels; ++l) {
536: DMLabel label;
537: const char *name;
538: IS valueIS;
539: const PetscInt *values;
540: PetscInt numValues, v;
542: DMPlexGetLabelName(dm, l, &name);
543: DMPlexGetLabel(dm, name, &label);
544: DMLabelGetNumValues(label, &numValues);
545: PetscViewerASCIIPrintf(viewer, " %s: %d strata of sizes (", name, numValues);
546: DMLabelGetValueIS(label, &valueIS);
547: ISGetIndices(valueIS, &values);
548: for (v = 0; v < numValues; ++v) {
549: PetscInt size;
551: DMLabelGetStratumSize(label, values[v], &size);
552: if (v > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
553: PetscViewerASCIIPrintf(viewer, "%d", size);
554: }
555: PetscViewerASCIIPrintf(viewer, ")\n");
556: ISRestoreIndices(valueIS, &values);
557: ISDestroy(&valueIS);
558: }
559: }
560: return(0);
561: }
565: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
566: {
567: PetscBool iascii, ishdf5, isvtk;
573: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
574: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
575: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
576: if (iascii) {
577: DMPlexView_Ascii(dm, viewer);
578: } else if (ishdf5) {
579: #if defined(PETSC_HAVE_HDF5)
580: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
581: DMPlexView_HDF5(dm, viewer);
582: PetscViewerPopFormat(viewer);
583: #else
584: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
585: #endif
586: }
587: else if (isvtk) {
588: DMPlexVTKWriteAll((PetscObject) dm,viewer);
589: }
590: return(0);
591: }
595: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
596: {
597: PetscBool isbinary, ishdf5;
603: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
604: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
605: if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");}
606: else if (ishdf5) {
607: #if defined(PETSC_HAVE_HDF5)
608: DMPlexLoad_HDF5(dm, viewer);
609: #else
610: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
611: #endif
612: }
613: return(0);
614: }
618: static PetscErrorCode BoundaryDestroy(DMBoundary *boundary)
619: {
620: DMBoundary b, next;
624: if (!boundary) return(0);
625: b = *boundary;
626: *boundary = NULL;
627: for (; b; b = next) {
628: next = b->next;
629: PetscFree(b->comps);
630: PetscFree(b->ids);
631: PetscFree(b->name);
632: PetscFree(b->labelname);
633: PetscFree(b);
634: }
635: return(0);
636: }
640: PetscErrorCode DMDestroy_Plex(DM dm)
641: {
642: DM_Plex *mesh = (DM_Plex*) dm->data;
643: PlexLabel next = mesh->labels;
647: if (--mesh->refct > 0) return(0);
648: PetscSectionDestroy(&mesh->coneSection);
649: PetscFree(mesh->cones);
650: PetscFree(mesh->coneOrientations);
651: PetscSectionDestroy(&mesh->supportSection);
652: PetscFree(mesh->supports);
653: PetscFree(mesh->facesTmp);
654: PetscFree(mesh->tetgenOpts);
655: PetscFree(mesh->triangleOpts);
656: PetscPartitionerDestroy(&mesh->partitioner);
657: while (next) {
658: PlexLabel tmp = next->next;
660: DMLabelDestroy(&next->label);
661: PetscFree(next);
662: next = tmp;
663: }
664: DMDestroy(&mesh->coarseMesh);
665: DMLabelDestroy(&mesh->subpointMap);
666: ISDestroy(&mesh->globalVertexNumbers);
667: ISDestroy(&mesh->globalCellNumbers);
668: BoundaryDestroy(&mesh->boundary);
669: PetscSectionDestroy(&mesh->anchorSection);
670: ISDestroy(&mesh->anchorIS);
671: PetscSectionDestroy(&mesh->parentSection);
672: PetscFree(mesh->parents);
673: PetscFree(mesh->childIDs);
674: PetscSectionDestroy(&mesh->childSection);
675: PetscFree(mesh->children);
676: DMDestroy(&mesh->referenceTree);
677: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
678: PetscFree(mesh);
679: return(0);
680: }
684: PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
685: {
686: PetscSection sectionGlobal;
687: PetscInt bs = -1;
688: PetscInt localSize;
689: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock;
691: MatType mtype;
692: ISLocalToGlobalMapping ltog;
695: MatInitializePackage();
696: mtype = dm->mattype;
697: DMGetDefaultGlobalSection(dm, §ionGlobal);
698: /* PetscSectionGetStorageSize(sectionGlobal, &localSize); */
699: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
700: MatCreate(PetscObjectComm((PetscObject)dm), J);
701: MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
702: MatSetType(*J, mtype);
703: MatSetFromOptions(*J);
704: PetscStrcmp(mtype, MATSHELL, &isShell);
705: PetscStrcmp(mtype, MATBAIJ, &isBlock);
706: PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
707: PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
708: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
709: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
710: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
711: if (!isShell) {
712: PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
713: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin;
715: if (bs < 0) {
716: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
717: PetscInt pStart, pEnd, p, dof, cdof;
719: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
720: for (p = pStart; p < pEnd; ++p) {
721: PetscSectionGetDof(sectionGlobal, p, &dof);
722: PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
723: if (dof-cdof) {
724: if (bs < 0) {
725: bs = dof-cdof;
726: } else if (bs != dof-cdof) {
727: /* Layout does not admit a pointwise block size */
728: bs = 1;
729: break;
730: }
731: }
732: }
733: /* Must have same blocksize on all procs (some might have no points) */
734: bsLocal = bs;
735: MPI_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
736: bsLocal = bs < 0 ? bsMax : bs;
737: MPI_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
738: if (bsMin != bsMax) {
739: bs = 1;
740: } else {
741: bs = bsMax;
742: }
743: } else {
744: bs = 1;
745: }
746: }
747: PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
748: DMPlexPreallocateOperator(dm, bs, dnz, onz, dnzu, onzu, *J, fillMatrix);
749: PetscFree4(dnz, onz, dnzu, onzu);
751: /* Set localtoglobalmapping on the matrix for MatSetValuesLocal() to work */
752: DMGetLocalToGlobalMapping(dm,<og);
753: MatSetLocalToGlobalMapping(*J,ltog,ltog);
754: }
755: return(0);
756: }
760: /*@
761: DMPlexGetChart - Return the interval for all mesh points [pStart, pEnd)
763: Not collective
765: Input Parameter:
766: . mesh - The DMPlex
768: Output Parameters:
769: + pStart - The first mesh point
770: - pEnd - The upper bound for mesh points
772: Level: beginner
774: .seealso: DMPlexCreate(), DMPlexSetChart()
775: @*/
776: PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
777: {
778: DM_Plex *mesh = (DM_Plex*) dm->data;
783: PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
784: return(0);
785: }
789: /*@
790: DMPlexSetChart - Set the interval for all mesh points [pStart, pEnd)
792: Not collective
794: Input Parameters:
795: + mesh - The DMPlex
796: . pStart - The first mesh point
797: - pEnd - The upper bound for mesh points
799: Output Parameters:
801: Level: beginner
803: .seealso: DMPlexCreate(), DMPlexGetChart()
804: @*/
805: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
806: {
807: DM_Plex *mesh = (DM_Plex*) dm->data;
812: PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
813: PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
814: return(0);
815: }
819: /*@
820: DMPlexGetConeSize - Return the number of in-edges for this point in the Sieve DAG
822: Not collective
824: Input Parameters:
825: + mesh - The DMPlex
826: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
828: Output Parameter:
829: . size - The cone size for point p
831: Level: beginner
833: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
834: @*/
835: PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size)
836: {
837: DM_Plex *mesh = (DM_Plex*) dm->data;
843: PetscSectionGetDof(mesh->coneSection, p, size);
844: return(0);
845: }
849: /*@
850: DMPlexSetConeSize - Set the number of in-edges for this point in the Sieve DAG
852: Not collective
854: Input Parameters:
855: + mesh - The DMPlex
856: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
857: - size - The cone size for point p
859: Output Parameter:
861: Note:
862: This should be called after DMPlexSetChart().
864: Level: beginner
866: .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart()
867: @*/
868: PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size)
869: {
870: DM_Plex *mesh = (DM_Plex*) dm->data;
875: PetscSectionSetDof(mesh->coneSection, p, size);
877: mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
878: return(0);
879: }
883: /*@
884: DMPlexAddConeSize - Add the given number of in-edges to this point in the Sieve DAG
886: Not collective
888: Input Parameters:
889: + mesh - The DMPlex
890: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
891: - size - The additional cone size for point p
893: Output Parameter:
895: Note:
896: This should be called after DMPlexSetChart().
898: Level: beginner
900: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart()
901: @*/
902: PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size)
903: {
904: DM_Plex *mesh = (DM_Plex*) dm->data;
905: PetscInt csize;
910: PetscSectionAddDof(mesh->coneSection, p, size);
911: PetscSectionGetDof(mesh->coneSection, p, &csize);
913: mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
914: return(0);
915: }
919: /*@C
920: DMPlexGetCone - Return the points on the in-edges for this point in the Sieve DAG
922: Not collective
924: Input Parameters:
925: + mesh - The DMPlex
926: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
928: Output Parameter:
929: . cone - An array of points which are on the in-edges for point p
931: Level: beginner
933: Fortran Notes:
934: Since it returns an array, this routine is only available in Fortran 90, and you must
935: include petsc.h90 in your code.
937: You must also call DMPlexRestoreCone() after you finish using the returned array.
939: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart()
940: @*/
941: PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
942: {
943: DM_Plex *mesh = (DM_Plex*) dm->data;
944: PetscInt off;
950: PetscSectionGetOffset(mesh->coneSection, p, &off);
951: *cone = &mesh->cones[off];
952: return(0);
953: }
957: /*@
958: DMPlexSetCone - Set the points on the in-edges for this point in the Sieve DAG
960: Not collective
962: Input Parameters:
963: + mesh - The DMPlex
964: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
965: - cone - An array of points which are on the in-edges for point p
967: Output Parameter:
969: Note:
970: This should be called after all calls to DMPlexSetConeSize() and DMSetUp().
972: Level: beginner
974: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
975: @*/
976: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
977: {
978: DM_Plex *mesh = (DM_Plex*) dm->data;
979: PetscInt pStart, pEnd;
980: PetscInt dof, off, c;
985: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
986: PetscSectionGetDof(mesh->coneSection, p, &dof);
988: PetscSectionGetOffset(mesh->coneSection, p, &off);
989: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
990: for (c = 0; c < dof; ++c) {
991: if ((cone[c] < pStart) || (cone[c] >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D, %D)", cone[c], pStart, pEnd);
992: mesh->cones[off+c] = cone[c];
993: }
994: return(0);
995: }
999: /*@C
1000: DMPlexGetConeOrientation - Return the orientations on the in-edges for this point in the Sieve DAG
1002: Not collective
1004: Input Parameters:
1005: + mesh - The DMPlex
1006: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1008: Output Parameter:
1009: . coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1010: integer giving the prescription for cone traversal. If it is negative, the cone is
1011: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1012: the index of the cone point on which to start.
1014: Level: beginner
1016: Fortran Notes:
1017: Since it returns an array, this routine is only available in Fortran 90, and you must
1018: include petsc.h90 in your code.
1020: You must also call DMPlexRestoreConeOrientation() after you finish using the returned array.
1022: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
1023: @*/
1024: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
1025: {
1026: DM_Plex *mesh = (DM_Plex*) dm->data;
1027: PetscInt off;
1032: #if defined(PETSC_USE_DEBUG)
1033: {
1034: PetscInt dof;
1035: PetscSectionGetDof(mesh->coneSection, p, &dof);
1037: }
1038: #endif
1039: PetscSectionGetOffset(mesh->coneSection, p, &off);
1041: *coneOrientation = &mesh->coneOrientations[off];
1042: return(0);
1043: }
1047: /*@
1048: DMPlexSetConeOrientation - Set the orientations on the in-edges for this point in the Sieve DAG
1050: Not collective
1052: Input Parameters:
1053: + mesh - The DMPlex
1054: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1055: - coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1056: integer giving the prescription for cone traversal. If it is negative, the cone is
1057: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1058: the index of the cone point on which to start.
1060: Output Parameter:
1062: Note:
1063: This should be called after all calls to DMPlexSetConeSize() and DMSetUp().
1065: Level: beginner
1067: .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
1068: @*/
1069: PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
1070: {
1071: DM_Plex *mesh = (DM_Plex*) dm->data;
1072: PetscInt pStart, pEnd;
1073: PetscInt dof, off, c;
1078: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1079: PetscSectionGetDof(mesh->coneSection, p, &dof);
1081: PetscSectionGetOffset(mesh->coneSection, p, &off);
1082: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1083: for (c = 0; c < dof; ++c) {
1084: PetscInt cdof, o = coneOrientation[c];
1086: PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1087: if (o && ((o < -(cdof+1)) || (o >= cdof))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone orientation %D is not in the valid range [%D. %D)", o, -(cdof+1), cdof);
1088: mesh->coneOrientations[off+c] = o;
1089: }
1090: return(0);
1091: }
1095: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
1096: {
1097: DM_Plex *mesh = (DM_Plex*) dm->data;
1098: PetscInt pStart, pEnd;
1099: PetscInt dof, off;
1104: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1105: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1106: if ((conePoint < pStart) || (conePoint >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D, %D)", conePoint, pStart, pEnd);
1107: PetscSectionGetDof(mesh->coneSection, p, &dof);
1108: PetscSectionGetOffset(mesh->coneSection, p, &off);
1109: if ((conePos < 0) || (conePos >= dof)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone position %D of point %D is not in the valid range [0, %D)", conePos, p, dof);
1110: mesh->cones[off+conePos] = conePoint;
1111: return(0);
1112: }
1116: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1117: {
1118: DM_Plex *mesh = (DM_Plex*) dm->data;
1119: PetscInt pStart, pEnd;
1120: PetscInt dof, off;
1125: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1126: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1127: PetscSectionGetDof(mesh->coneSection, p, &dof);
1128: PetscSectionGetOffset(mesh->coneSection, p, &off);
1129: if ((conePos < 0) || (conePos >= dof)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone position %D of point %D is not in the valid range [0, %D)", conePos, p, dof);
1130: mesh->coneOrientations[off+conePos] = coneOrientation;
1131: return(0);
1132: }
1136: /*@
1137: DMPlexGetSupportSize - Return the number of out-edges for this point in the Sieve DAG
1139: Not collective
1141: Input Parameters:
1142: + mesh - The DMPlex
1143: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1145: Output Parameter:
1146: . size - The support size for point p
1148: Level: beginner
1150: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1151: @*/
1152: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1153: {
1154: DM_Plex *mesh = (DM_Plex*) dm->data;
1160: PetscSectionGetDof(mesh->supportSection, p, size);
1161: return(0);
1162: }
1166: /*@
1167: DMPlexSetSupportSize - Set the number of out-edges for this point in the Sieve DAG
1169: Not collective
1171: Input Parameters:
1172: + mesh - The DMPlex
1173: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1174: - size - The support size for point p
1176: Output Parameter:
1178: Note:
1179: This should be called after DMPlexSetChart().
1181: Level: beginner
1183: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1184: @*/
1185: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1186: {
1187: DM_Plex *mesh = (DM_Plex*) dm->data;
1192: PetscSectionSetDof(mesh->supportSection, p, size);
1194: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1195: return(0);
1196: }
1200: /*@C
1201: DMPlexGetSupport - Return the points on the out-edges for this point in the Sieve DAG
1203: Not collective
1205: Input Parameters:
1206: + mesh - The DMPlex
1207: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1209: Output Parameter:
1210: . support - An array of points which are on the out-edges for point p
1212: Level: beginner
1214: Fortran Notes:
1215: Since it returns an array, this routine is only available in Fortran 90, and you must
1216: include petsc.h90 in your code.
1218: You must also call DMPlexRestoreSupport() after you finish using the returned array.
1220: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1221: @*/
1222: PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1223: {
1224: DM_Plex *mesh = (DM_Plex*) dm->data;
1225: PetscInt off;
1231: PetscSectionGetOffset(mesh->supportSection, p, &off);
1232: *support = &mesh->supports[off];
1233: return(0);
1234: }
1238: /*@
1239: DMPlexSetSupport - Set the points on the out-edges for this point in the Sieve DAG
1241: Not collective
1243: Input Parameters:
1244: + mesh - The DMPlex
1245: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1246: - support - An array of points which are on the in-edges for point p
1248: Output Parameter:
1250: Note:
1251: This should be called after all calls to DMPlexSetSupportSize() and DMSetUp().
1253: Level: beginner
1255: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1256: @*/
1257: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1258: {
1259: DM_Plex *mesh = (DM_Plex*) dm->data;
1260: PetscInt pStart, pEnd;
1261: PetscInt dof, off, c;
1266: PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1267: PetscSectionGetDof(mesh->supportSection, p, &dof);
1269: PetscSectionGetOffset(mesh->supportSection, p, &off);
1270: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1271: for (c = 0; c < dof; ++c) {
1272: if ((support[c] < pStart) || (support[c] >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support point %D is not in the valid range [%D, %D)", support[c], pStart, pEnd);
1273: mesh->supports[off+c] = support[c];
1274: }
1275: return(0);
1276: }
1280: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1281: {
1282: DM_Plex *mesh = (DM_Plex*) dm->data;
1283: PetscInt pStart, pEnd;
1284: PetscInt dof, off;
1289: PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1290: PetscSectionGetDof(mesh->supportSection, p, &dof);
1291: PetscSectionGetOffset(mesh->supportSection, p, &off);
1292: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1293: if ((supportPoint < pStart) || (supportPoint >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support point %D is not in the valid range [%D, %D)", supportPoint, pStart, pEnd);
1294: if (supportPos >= dof) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support position %D of point %D is not in the valid range [0, %D)", supportPos, p, dof);
1295: mesh->supports[off+supportPos] = supportPoint;
1296: return(0);
1297: }
1301: /*@C
1302: DMPlexGetTransitiveClosure - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1304: Not collective
1306: Input Parameters:
1307: + mesh - The DMPlex
1308: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1309: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1310: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
1312: Output Parameters:
1313: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1314: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
1316: Note:
1317: If using internal storage (points is NULL on input), each call overwrites the last output.
1319: Fortran Notes:
1320: Since it returns an array, this routine is only available in Fortran 90, and you must
1321: include petsc.h90 in your code.
1323: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1325: Level: beginner
1327: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1328: @*/
1329: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1330: {
1331: DM_Plex *mesh = (DM_Plex*) dm->data;
1332: PetscInt *closure, *fifo;
1333: const PetscInt *tmp = NULL, *tmpO = NULL;
1334: PetscInt tmpSize, t;
1335: PetscInt depth = 0, maxSize;
1336: PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0;
1337: PetscErrorCode ierr;
1341: DMPlexGetDepth(dm, &depth);
1342: /* This is only 1-level */
1343: if (useCone) {
1344: DMPlexGetConeSize(dm, p, &tmpSize);
1345: DMPlexGetCone(dm, p, &tmp);
1346: DMPlexGetConeOrientation(dm, p, &tmpO);
1347: } else {
1348: DMPlexGetSupportSize(dm, p, &tmpSize);
1349: DMPlexGetSupport(dm, p, &tmp);
1350: }
1351: if (depth == 1) {
1352: if (*points) {
1353: closure = *points;
1354: } else {
1355: maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1356: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1357: }
1358: closure[0] = p; closure[1] = 0;
1359: for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1360: closure[closureSize] = tmp[t];
1361: closure[closureSize+1] = tmpO ? tmpO[t] : 0;
1362: }
1363: if (numPoints) *numPoints = closureSize/2;
1364: if (points) *points = closure;
1365: return(0);
1366: }
1367: {
1368: PetscInt c, coneSeries, s,supportSeries;
1370: c = mesh->maxConeSize;
1371: coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1372: s = mesh->maxSupportSize;
1373: supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1374: maxSize = 2*PetscMax(coneSeries,supportSeries);
1375: }
1376: DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1377: if (*points) {
1378: closure = *points;
1379: } else {
1380: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1381: }
1382: closure[0] = p; closure[1] = 0;
1383: for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1384: const PetscInt cp = tmp[t];
1385: const PetscInt co = tmpO ? tmpO[t] : 0;
1387: closure[closureSize] = cp;
1388: closure[closureSize+1] = co;
1389: fifo[fifoSize] = cp;
1390: fifo[fifoSize+1] = co;
1391: }
1392: /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1393: while (fifoSize - fifoStart) {
1394: const PetscInt q = fifo[fifoStart];
1395: const PetscInt o = fifo[fifoStart+1];
1396: const PetscInt rev = o >= 0 ? 0 : 1;
1397: const PetscInt off = rev ? -(o+1) : o;
1399: if (useCone) {
1400: DMPlexGetConeSize(dm, q, &tmpSize);
1401: DMPlexGetCone(dm, q, &tmp);
1402: DMPlexGetConeOrientation(dm, q, &tmpO);
1403: } else {
1404: DMPlexGetSupportSize(dm, q, &tmpSize);
1405: DMPlexGetSupport(dm, q, &tmp);
1406: tmpO = NULL;
1407: }
1408: for (t = 0; t < tmpSize; ++t) {
1409: const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize;
1410: const PetscInt cp = tmp[i];
1411: /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1412: /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1413: const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1414: PetscInt co = tmpO ? tmpO[i] : 0;
1415: PetscInt c;
1417: if (rev) {
1418: PetscInt childSize, coff;
1419: DMPlexGetConeSize(dm, cp, &childSize);
1420: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1421: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1422: }
1423: /* Check for duplicate */
1424: for (c = 0; c < closureSize; c += 2) {
1425: if (closure[c] == cp) break;
1426: }
1427: if (c == closureSize) {
1428: closure[closureSize] = cp;
1429: closure[closureSize+1] = co;
1430: fifo[fifoSize] = cp;
1431: fifo[fifoSize+1] = co;
1432: closureSize += 2;
1433: fifoSize += 2;
1434: }
1435: }
1436: fifoStart += 2;
1437: }
1438: if (numPoints) *numPoints = closureSize/2;
1439: if (points) *points = closure;
1440: DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1441: return(0);
1442: }
1446: /*@C
1447: DMPlexGetTransitiveClosure_Internal - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG with a specified initial orientation
1449: Not collective
1451: Input Parameters:
1452: + mesh - The DMPlex
1453: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1454: . orientation - The orientation of the point
1455: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1456: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
1458: Output Parameters:
1459: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1460: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
1462: Note:
1463: If using internal storage (points is NULL on input), each call overwrites the last output.
1465: Fortran Notes:
1466: Since it returns an array, this routine is only available in Fortran 90, and you must
1467: include petsc.h90 in your code.
1469: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1471: Level: beginner
1473: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1474: @*/
1475: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1476: {
1477: DM_Plex *mesh = (DM_Plex*) dm->data;
1478: PetscInt *closure, *fifo;
1479: const PetscInt *tmp = NULL, *tmpO = NULL;
1480: PetscInt tmpSize, t;
1481: PetscInt depth = 0, maxSize;
1482: PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0;
1483: PetscErrorCode ierr;
1487: DMPlexGetDepth(dm, &depth);
1488: /* This is only 1-level */
1489: if (useCone) {
1490: DMPlexGetConeSize(dm, p, &tmpSize);
1491: DMPlexGetCone(dm, p, &tmp);
1492: DMPlexGetConeOrientation(dm, p, &tmpO);
1493: } else {
1494: DMPlexGetSupportSize(dm, p, &tmpSize);
1495: DMPlexGetSupport(dm, p, &tmp);
1496: }
1497: if (depth == 1) {
1498: if (*points) {
1499: closure = *points;
1500: } else {
1501: maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1502: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1503: }
1504: closure[0] = p; closure[1] = ornt;
1505: for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1506: const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1507: closure[closureSize] = tmp[i];
1508: closure[closureSize+1] = tmpO ? tmpO[i] : 0;
1509: }
1510: if (numPoints) *numPoints = closureSize/2;
1511: if (points) *points = closure;
1512: return(0);
1513: }
1514: {
1515: PetscInt c, coneSeries, s,supportSeries;
1517: c = mesh->maxConeSize;
1518: coneSeries = (c > 1) ? ((PetscPowInt(c,depth+1)-1)/(c-1)) : depth+1;
1519: s = mesh->maxSupportSize;
1520: supportSeries = (s > 1) ? ((PetscPowInt(s,depth+1)-1)/(s-1)) : depth+1;
1521: maxSize = 2*PetscMax(coneSeries,supportSeries);
1522: }
1523: DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1524: if (*points) {
1525: closure = *points;
1526: } else {
1527: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1528: }
1529: closure[0] = p; closure[1] = ornt;
1530: for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1531: const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1532: const PetscInt cp = tmp[i];
1533: PetscInt co = tmpO ? tmpO[i] : 0;
1535: if (ornt < 0) {
1536: PetscInt childSize, coff;
1537: DMPlexGetConeSize(dm, cp, &childSize);
1538: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1539: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1540: }
1541: closure[closureSize] = cp;
1542: closure[closureSize+1] = co;
1543: fifo[fifoSize] = cp;
1544: fifo[fifoSize+1] = co;
1545: }
1546: /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1547: while (fifoSize - fifoStart) {
1548: const PetscInt q = fifo[fifoStart];
1549: const PetscInt o = fifo[fifoStart+1];
1550: const PetscInt rev = o >= 0 ? 0 : 1;
1551: const PetscInt off = rev ? -(o+1) : o;
1553: if (useCone) {
1554: DMPlexGetConeSize(dm, q, &tmpSize);
1555: DMPlexGetCone(dm, q, &tmp);
1556: DMPlexGetConeOrientation(dm, q, &tmpO);
1557: } else {
1558: DMPlexGetSupportSize(dm, q, &tmpSize);
1559: DMPlexGetSupport(dm, q, &tmp);
1560: tmpO = NULL;
1561: }
1562: for (t = 0; t < tmpSize; ++t) {
1563: const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize;
1564: const PetscInt cp = tmp[i];
1565: /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1566: /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1567: const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1568: PetscInt co = tmpO ? tmpO[i] : 0;
1569: PetscInt c;
1571: if (rev) {
1572: PetscInt childSize, coff;
1573: DMPlexGetConeSize(dm, cp, &childSize);
1574: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1575: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1576: }
1577: /* Check for duplicate */
1578: for (c = 0; c < closureSize; c += 2) {
1579: if (closure[c] == cp) break;
1580: }
1581: if (c == closureSize) {
1582: closure[closureSize] = cp;
1583: closure[closureSize+1] = co;
1584: fifo[fifoSize] = cp;
1585: fifo[fifoSize+1] = co;
1586: closureSize += 2;
1587: fifoSize += 2;
1588: }
1589: }
1590: fifoStart += 2;
1591: }
1592: if (numPoints) *numPoints = closureSize/2;
1593: if (points) *points = closure;
1594: DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1595: return(0);
1596: }
1600: /*@C
1601: DMPlexRestoreTransitiveClosure - Restore the array of points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1603: Not collective
1605: Input Parameters:
1606: + mesh - The DMPlex
1607: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1608: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1609: . numPoints - The number of points in the closure, so points[] is of size 2*numPoints, zeroed on exit
1610: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...], zeroed on exit
1612: Note:
1613: If not using internal storage (points is not NULL on input), this call is unnecessary
1615: Fortran Notes:
1616: Since it returns an array, this routine is only available in Fortran 90, and you must
1617: include petsc.h90 in your code.
1619: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1621: Level: beginner
1623: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1624: @*/
1625: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1626: {
1633: DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1634: if (numPoints) *numPoints = 0;
1635: return(0);
1636: }
1640: /*@
1641: DMPlexGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG
1643: Not collective
1645: Input Parameter:
1646: . mesh - The DMPlex
1648: Output Parameters:
1649: + maxConeSize - The maximum number of in-edges
1650: - maxSupportSize - The maximum number of out-edges
1652: Level: beginner
1654: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1655: @*/
1656: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1657: {
1658: DM_Plex *mesh = (DM_Plex*) dm->data;
1662: if (maxConeSize) *maxConeSize = mesh->maxConeSize;
1663: if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1664: return(0);
1665: }
1669: PetscErrorCode DMSetUp_Plex(DM dm)
1670: {
1671: DM_Plex *mesh = (DM_Plex*) dm->data;
1672: PetscInt size;
1677: PetscSectionSetUp(mesh->coneSection);
1678: PetscSectionGetStorageSize(mesh->coneSection, &size);
1679: PetscMalloc1(size, &mesh->cones);
1680: PetscCalloc1(size, &mesh->coneOrientations);
1681: if (mesh->maxSupportSize) {
1682: PetscSectionSetUp(mesh->supportSection);
1683: PetscSectionGetStorageSize(mesh->supportSection, &size);
1684: PetscMalloc1(size, &mesh->supports);
1685: }
1686: return(0);
1687: }
1691: PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1692: {
1696: if (subdm) {DMClone(dm, subdm);}
1697: DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1698: return(0);
1699: }
1703: /*@
1704: DMPlexSymmetrize - Creates support (out-edge) information from cone (in-edge) inoformation
1706: Not collective
1708: Input Parameter:
1709: . mesh - The DMPlex
1711: Output Parameter:
1713: Note:
1714: This should be called after all calls to DMPlexSetCone()
1716: Level: beginner
1718: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1719: @*/
1720: PetscErrorCode DMPlexSymmetrize(DM dm)
1721: {
1722: DM_Plex *mesh = (DM_Plex*) dm->data;
1723: PetscInt *offsets;
1724: PetscInt supportSize;
1725: PetscInt pStart, pEnd, p;
1730: if (mesh->supports) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Supports were already setup in this DMPlex");
1731: /* Calculate support sizes */
1732: DMPlexGetChart(dm, &pStart, &pEnd);
1733: for (p = pStart; p < pEnd; ++p) {
1734: PetscInt dof, off, c;
1736: PetscSectionGetDof(mesh->coneSection, p, &dof);
1737: PetscSectionGetOffset(mesh->coneSection, p, &off);
1738: for (c = off; c < off+dof; ++c) {
1739: PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1740: }
1741: }
1742: for (p = pStart; p < pEnd; ++p) {
1743: PetscInt dof;
1745: PetscSectionGetDof(mesh->supportSection, p, &dof);
1747: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1748: }
1749: PetscSectionSetUp(mesh->supportSection);
1750: /* Calculate supports */
1751: PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1752: PetscMalloc1(supportSize, &mesh->supports);
1753: PetscCalloc1(pEnd - pStart, &offsets);
1754: for (p = pStart; p < pEnd; ++p) {
1755: PetscInt dof, off, c;
1757: PetscSectionGetDof(mesh->coneSection, p, &dof);
1758: PetscSectionGetOffset(mesh->coneSection, p, &off);
1759: for (c = off; c < off+dof; ++c) {
1760: const PetscInt q = mesh->cones[c];
1761: PetscInt offS;
1763: PetscSectionGetOffset(mesh->supportSection, q, &offS);
1765: mesh->supports[offS+offsets[q]] = p;
1766: ++offsets[q];
1767: }
1768: }
1769: PetscFree(offsets);
1770: return(0);
1771: }
1775: /*@
1776: DMPlexStratify - The Sieve DAG for most topologies is a graded poset (http://en.wikipedia.org/wiki/Graded_poset), and
1777: can be illustrated by Hasse Diagram (a http://en.wikipedia.org/wiki/Hasse_diagram). The strata group all points of the
1778: same grade, and this function calculates the strata. This grade can be seen as the height (or depth) of the point in
1779: the DAG.
1781: Not collective
1783: Input Parameter:
1784: . mesh - The DMPlex
1786: Output Parameter:
1788: Notes:
1789: Concretely, DMPlexStratify() creates a new label named "depth" containing the dimension of each element: 0 for vertices,
1790: 1 for edges, and so on. The depth label can be accessed through DMPlexGetDepthLabel() or DMPlexGetDepthStratum(), or
1791: manually via DMPlexGetLabel(). The height is defined implicitly by height = maxDimension - depth, and can be accessed
1792: via DMPlexGetHeightStratum(). For example, cells have height 0 and faces have height 1.
1794: DMPlexStratify() should be called after all calls to DMPlexSymmetrize()
1796: Level: beginner
1798: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1799: @*/
1800: PetscErrorCode DMPlexStratify(DM dm)
1801: {
1802: DMLabel label;
1803: PetscInt pStart, pEnd, p;
1804: PetscInt numRoots = 0, numLeaves = 0;
1809: PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1810: /* Calculate depth */
1811: DMPlexGetChart(dm, &pStart, &pEnd);
1812: DMPlexCreateLabel(dm, "depth");
1813: DMPlexGetDepthLabel(dm, &label);
1814: /* Initialize roots and count leaves */
1815: for (p = pStart; p < pEnd; ++p) {
1816: PetscInt coneSize, supportSize;
1818: DMPlexGetConeSize(dm, p, &coneSize);
1819: DMPlexGetSupportSize(dm, p, &supportSize);
1820: if (!coneSize && supportSize) {
1821: ++numRoots;
1822: DMLabelSetValue(label, p, 0);
1823: } else if (!supportSize && coneSize) {
1824: ++numLeaves;
1825: } else if (!supportSize && !coneSize) {
1826: /* Isolated points */
1827: DMLabelSetValue(label, p, 0);
1828: }
1829: }
1830: if (numRoots + numLeaves == (pEnd - pStart)) {
1831: for (p = pStart; p < pEnd; ++p) {
1832: PetscInt coneSize, supportSize;
1834: DMPlexGetConeSize(dm, p, &coneSize);
1835: DMPlexGetSupportSize(dm, p, &supportSize);
1836: if (!supportSize && coneSize) {
1837: DMLabelSetValue(label, p, 1);
1838: }
1839: }
1840: } else {
1841: IS pointIS;
1842: PetscInt numPoints = 0, level = 0;
1844: DMLabelGetStratumIS(label, level, &pointIS);
1845: if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1846: while (numPoints) {
1847: const PetscInt *points;
1848: const PetscInt newLevel = level+1;
1850: ISGetIndices(pointIS, &points);
1851: for (p = 0; p < numPoints; ++p) {
1852: const PetscInt point = points[p];
1853: const PetscInt *support;
1854: PetscInt supportSize, s;
1856: DMPlexGetSupportSize(dm, point, &supportSize);
1857: DMPlexGetSupport(dm, point, &support);
1858: for (s = 0; s < supportSize; ++s) {
1859: DMLabelSetValue(label, support[s], newLevel);
1860: }
1861: }
1862: ISRestoreIndices(pointIS, &points);
1863: ++level;
1864: ISDestroy(&pointIS);
1865: DMLabelGetStratumIS(label, level, &pointIS);
1866: if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1867: else {numPoints = 0;}
1868: }
1869: ISDestroy(&pointIS);
1870: }
1871: PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1872: return(0);
1873: }
1877: /*@C
1878: DMPlexGetJoin - Get an array for the join of the set of points
1880: Not Collective
1882: Input Parameters:
1883: + dm - The DMPlex object
1884: . numPoints - The number of input points for the join
1885: - points - The input points
1887: Output Parameters:
1888: + numCoveredPoints - The number of points in the join
1889: - coveredPoints - The points in the join
1891: Level: intermediate
1893: Note: Currently, this is restricted to a single level join
1895: Fortran Notes:
1896: Since it returns an array, this routine is only available in Fortran 90, and you must
1897: include petsc.h90 in your code.
1899: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1901: .keywords: mesh
1902: .seealso: DMPlexRestoreJoin(), DMPlexGetMeet()
1903: @*/
1904: PetscErrorCode DMPlexGetJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1905: {
1906: DM_Plex *mesh = (DM_Plex*) dm->data;
1907: PetscInt *join[2];
1908: PetscInt joinSize, i = 0;
1909: PetscInt dof, off, p, c, m;
1917: DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[0]);
1918: DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1]);
1919: /* Copy in support of first point */
1920: PetscSectionGetDof(mesh->supportSection, points[0], &dof);
1921: PetscSectionGetOffset(mesh->supportSection, points[0], &off);
1922: for (joinSize = 0; joinSize < dof; ++joinSize) {
1923: join[i][joinSize] = mesh->supports[off+joinSize];
1924: }
1925: /* Check each successive support */
1926: for (p = 1; p < numPoints; ++p) {
1927: PetscInt newJoinSize = 0;
1929: PetscSectionGetDof(mesh->supportSection, points[p], &dof);
1930: PetscSectionGetOffset(mesh->supportSection, points[p], &off);
1931: for (c = 0; c < dof; ++c) {
1932: const PetscInt point = mesh->supports[off+c];
1934: for (m = 0; m < joinSize; ++m) {
1935: if (point == join[i][m]) {
1936: join[1-i][newJoinSize++] = point;
1937: break;
1938: }
1939: }
1940: }
1941: joinSize = newJoinSize;
1942: i = 1-i;
1943: }
1944: *numCoveredPoints = joinSize;
1945: *coveredPoints = join[i];
1946: DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
1947: return(0);
1948: }
1952: /*@C
1953: DMPlexRestoreJoin - Restore an array for the join of the set of points
1955: Not Collective
1957: Input Parameters:
1958: + dm - The DMPlex object
1959: . numPoints - The number of input points for the join
1960: - points - The input points
1962: Output Parameters:
1963: + numCoveredPoints - The number of points in the join
1964: - coveredPoints - The points in the join
1966: Fortran Notes:
1967: Since it returns an array, this routine is only available in Fortran 90, and you must
1968: include petsc.h90 in your code.
1970: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1972: Level: intermediate
1974: .keywords: mesh
1975: .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet()
1976: @*/
1977: PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1978: {
1986: DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
1987: if (numCoveredPoints) *numCoveredPoints = 0;
1988: return(0);
1989: }
1993: /*@C
1994: DMPlexGetFullJoin - Get an array for the join of the set of points
1996: Not Collective
1998: Input Parameters:
1999: + dm - The DMPlex object
2000: . numPoints - The number of input points for the join
2001: - points - The input points
2003: Output Parameters:
2004: + numCoveredPoints - The number of points in the join
2005: - coveredPoints - The points in the join
2007: Fortran Notes:
2008: Since it returns an array, this routine is only available in Fortran 90, and you must
2009: include petsc.h90 in your code.
2011: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2013: Level: intermediate
2015: .keywords: mesh
2016: .seealso: DMPlexGetJoin(), DMPlexRestoreJoin(), DMPlexGetMeet()
2017: @*/
2018: PetscErrorCode DMPlexGetFullJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2019: {
2020: DM_Plex *mesh = (DM_Plex*) dm->data;
2021: PetscInt *offsets, **closures;
2022: PetscInt *join[2];
2023: PetscInt depth = 0, maxSize, joinSize = 0, i = 0;
2024: PetscInt p, d, c, m, ms;
2033: DMPlexGetDepth(dm, &depth);
2034: PetscCalloc1(numPoints, &closures);
2035: DMGetWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2036: ms = mesh->maxSupportSize;
2037: maxSize = (ms > 1) ? ((PetscPowInt(ms,depth+1)-1)/(ms-1)) : depth + 1;
2038: DMGetWorkArray(dm, maxSize, PETSC_INT, &join[0]);
2039: DMGetWorkArray(dm, maxSize, PETSC_INT, &join[1]);
2041: for (p = 0; p < numPoints; ++p) {
2042: PetscInt closureSize;
2044: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closures[p]);
2046: offsets[p*(depth+2)+0] = 0;
2047: for (d = 0; d < depth+1; ++d) {
2048: PetscInt pStart, pEnd, i;
2050: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
2051: for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) {
2052: if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2053: offsets[p*(depth+2)+d+1] = i;
2054: break;
2055: }
2056: }
2057: if (i == closureSize) offsets[p*(depth+2)+d+1] = i;
2058: }
2059: if (offsets[p*(depth+2)+depth+1] != closureSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Total size of closure %D should be %D", offsets[p*(depth+2)+depth+1], closureSize);
2060: }
2061: for (d = 0; d < depth+1; ++d) {
2062: PetscInt dof;
2064: /* Copy in support of first point */
2065: dof = offsets[d+1] - offsets[d];
2066: for (joinSize = 0; joinSize < dof; ++joinSize) {
2067: join[i][joinSize] = closures[0][(offsets[d]+joinSize)*2];
2068: }
2069: /* Check each successive cone */
2070: for (p = 1; p < numPoints && joinSize; ++p) {
2071: PetscInt newJoinSize = 0;
2073: dof = offsets[p*(depth+2)+d+1] - offsets[p*(depth+2)+d];
2074: for (c = 0; c < dof; ++c) {
2075: const PetscInt point = closures[p][(offsets[p*(depth+2)+d]+c)*2];
2077: for (m = 0; m < joinSize; ++m) {
2078: if (point == join[i][m]) {
2079: join[1-i][newJoinSize++] = point;
2080: break;
2081: }
2082: }
2083: }
2084: joinSize = newJoinSize;
2085: i = 1-i;
2086: }
2087: if (joinSize) break;
2088: }
2089: *numCoveredPoints = joinSize;
2090: *coveredPoints = join[i];
2091: for (p = 0; p < numPoints; ++p) {
2092: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, NULL, &closures[p]);
2093: }
2094: PetscFree(closures);
2095: DMRestoreWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
2096: DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
2097: return(0);
2098: }
2102: /*@C
2103: DMPlexGetMeet - Get an array for the meet of the set of points
2105: Not Collective
2107: Input Parameters:
2108: + dm - The DMPlex object
2109: . numPoints - The number of input points for the meet
2110: - points - The input points
2112: Output Parameters:
2113: + numCoveredPoints - The number of points in the meet
2114: - coveredPoints - The points in the meet
2116: Level: intermediate
2118: Note: Currently, this is restricted to a single level meet
2120: Fortran Notes:
2121: Since it returns an array, this routine is only available in Fortran 90, and you must
2122: include petsc.h90 in your code.
2124: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2126: .keywords: mesh
2127: .seealso: DMPlexRestoreMeet(), DMPlexGetJoin()
2128: @*/
2129: PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2130: {
2131: DM_Plex *mesh = (DM_Plex*) dm->data;
2132: PetscInt *meet[2];
2133: PetscInt meetSize, i = 0;
2134: PetscInt dof, off, p, c, m;
2142: DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[0]);
2143: DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1]);
2144: /* Copy in cone of first point */
2145: PetscSectionGetDof(mesh->coneSection, points[0], &dof);
2146: PetscSectionGetOffset(mesh->coneSection, points[0], &off);
2147: for (meetSize = 0; meetSize < dof; ++meetSize) {
2148: meet[i][meetSize] = mesh->cones[off+meetSize];
2149: }
2150: /* Check each successive cone */
2151: for (p = 1; p < numPoints; ++p) {
2152: PetscInt newMeetSize = 0;
2154: PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2155: PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2156: for (c = 0; c < dof; ++c) {
2157: const PetscInt point = mesh->cones[off+c];
2159: for (m = 0; m < meetSize; ++m) {
2160: if (point == meet[i][m]) {
2161: meet[1-i][newMeetSize++] = point;
2162: break;
2163: }
2164: }
2165: }
2166: meetSize = newMeetSize;
2167: i = 1-i;
2168: }
2169: *numCoveringPoints = meetSize;
2170: *coveringPoints = meet[i];
2171: DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2172: return(0);
2173: }
2177: /*@C
2178: DMPlexRestoreMeet - Restore an array for the meet of the set of points
2180: Not Collective
2182: Input Parameters:
2183: + dm - The DMPlex object
2184: . numPoints - The number of input points for the meet
2185: - points - The input points
2187: Output Parameters:
2188: + numCoveredPoints - The number of points in the meet
2189: - coveredPoints - The points in the meet
2191: Level: intermediate
2193: Fortran Notes:
2194: Since it returns an array, this routine is only available in Fortran 90, and you must
2195: include petsc.h90 in your code.
2197: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2199: .keywords: mesh
2200: .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin()
2201: @*/
2202: PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2203: {
2211: DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2212: if (numCoveredPoints) *numCoveredPoints = 0;
2213: return(0);
2214: }
2218: /*@C
2219: DMPlexGetFullMeet - Get an array for the meet of the set of points
2221: Not Collective
2223: Input Parameters:
2224: + dm - The DMPlex object
2225: . numPoints - The number of input points for the meet
2226: - points - The input points
2228: Output Parameters:
2229: + numCoveredPoints - The number of points in the meet
2230: - coveredPoints - The points in the meet
2232: Level: intermediate
2234: Fortran Notes:
2235: Since it returns an array, this routine is only available in Fortran 90, and you must
2236: include petsc.h90 in your code.
2238: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2240: .keywords: mesh
2241: .seealso: DMPlexGetMeet(), DMPlexRestoreMeet(), DMPlexGetJoin()
2242: @*/
2243: PetscErrorCode DMPlexGetFullMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2244: {
2245: DM_Plex *mesh = (DM_Plex*) dm->data;
2246: PetscInt *offsets, **closures;
2247: PetscInt *meet[2];
2248: PetscInt height = 0, maxSize, meetSize = 0, i = 0;
2249: PetscInt p, h, c, m, mc;
2258: DMPlexGetDepth(dm, &height);
2259: PetscMalloc1(numPoints, &closures);
2260: DMGetWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2261: mc = mesh->maxConeSize;
2262: maxSize = (mc > 1) ? ((PetscPowInt(mc,height+1)-1)/(mc-1)) : height + 1;
2263: DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[0]);
2264: DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[1]);
2266: for (p = 0; p < numPoints; ++p) {
2267: PetscInt closureSize;
2269: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closures[p]);
2271: offsets[p*(height+2)+0] = 0;
2272: for (h = 0; h < height+1; ++h) {
2273: PetscInt pStart, pEnd, i;
2275: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
2276: for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) {
2277: if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2278: offsets[p*(height+2)+h+1] = i;
2279: break;
2280: }
2281: }
2282: if (i == closureSize) offsets[p*(height+2)+h+1] = i;
2283: }
2284: if (offsets[p*(height+2)+height+1] != closureSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Total size of closure %D should be %D", offsets[p*(height+2)+height+1], closureSize);
2285: }
2286: for (h = 0; h < height+1; ++h) {
2287: PetscInt dof;
2289: /* Copy in cone of first point */
2290: dof = offsets[h+1] - offsets[h];
2291: for (meetSize = 0; meetSize < dof; ++meetSize) {
2292: meet[i][meetSize] = closures[0][(offsets[h]+meetSize)*2];
2293: }
2294: /* Check each successive cone */
2295: for (p = 1; p < numPoints && meetSize; ++p) {
2296: PetscInt newMeetSize = 0;
2298: dof = offsets[p*(height+2)+h+1] - offsets[p*(height+2)+h];
2299: for (c = 0; c < dof; ++c) {
2300: const PetscInt point = closures[p][(offsets[p*(height+2)+h]+c)*2];
2302: for (m = 0; m < meetSize; ++m) {
2303: if (point == meet[i][m]) {
2304: meet[1-i][newMeetSize++] = point;
2305: break;
2306: }
2307: }
2308: }
2309: meetSize = newMeetSize;
2310: i = 1-i;
2311: }
2312: if (meetSize) break;
2313: }
2314: *numCoveredPoints = meetSize;
2315: *coveredPoints = meet[i];
2316: for (p = 0; p < numPoints; ++p) {
2317: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, NULL, &closures[p]);
2318: }
2319: PetscFree(closures);
2320: DMRestoreWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2321: DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2322: return(0);
2323: }
2327: /*@C
2328: DMPlexEqual - Determine if two DMs have the same topology
2330: Not Collective
2332: Input Parameters:
2333: + dmA - A DMPlex object
2334: - dmB - A DMPlex object
2336: Output Parameters:
2337: . equal - PETSC_TRUE if the topologies are identical
2339: Level: intermediate
2341: Notes:
2342: We are not solving graph isomorphism, so we do not permutation.
2344: .keywords: mesh
2345: .seealso: DMPlexGetCone()
2346: @*/
2347: PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal)
2348: {
2349: PetscInt depth, depthB, pStart, pEnd, pStartB, pEndB, p;
2357: *equal = PETSC_FALSE;
2358: DMPlexGetDepth(dmA, &depth);
2359: DMPlexGetDepth(dmB, &depthB);
2360: if (depth != depthB) return(0);
2361: DMPlexGetChart(dmA, &pStart, &pEnd);
2362: DMPlexGetChart(dmB, &pStartB, &pEndB);
2363: if ((pStart != pStartB) || (pEnd != pEndB)) return(0);
2364: for (p = pStart; p < pEnd; ++p) {
2365: const PetscInt *cone, *coneB, *ornt, *orntB, *support, *supportB;
2366: PetscInt coneSize, coneSizeB, c, supportSize, supportSizeB, s;
2368: DMPlexGetConeSize(dmA, p, &coneSize);
2369: DMPlexGetCone(dmA, p, &cone);
2370: DMPlexGetConeOrientation(dmA, p, &ornt);
2371: DMPlexGetConeSize(dmB, p, &coneSizeB);
2372: DMPlexGetCone(dmB, p, &coneB);
2373: DMPlexGetConeOrientation(dmB, p, &orntB);
2374: if (coneSize != coneSizeB) return(0);
2375: for (c = 0; c < coneSize; ++c) {
2376: if (cone[c] != coneB[c]) return(0);
2377: if (ornt[c] != orntB[c]) return(0);
2378: }
2379: DMPlexGetSupportSize(dmA, p, &supportSize);
2380: DMPlexGetSupport(dmA, p, &support);
2381: DMPlexGetSupportSize(dmB, p, &supportSizeB);
2382: DMPlexGetSupport(dmB, p, &supportB);
2383: if (supportSize != supportSizeB) return(0);
2384: for (s = 0; s < supportSize; ++s) {
2385: if (support[s] != supportB[s]) return(0);
2386: }
2387: }
2388: *equal = PETSC_TRUE;
2389: return(0);
2390: }
2394: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2395: {
2396: MPI_Comm comm;
2400: PetscObjectGetComm((PetscObject)dm,&comm);
2402: switch (cellDim) {
2403: case 0:
2404: *numFaceVertices = 0;
2405: break;
2406: case 1:
2407: *numFaceVertices = 1;
2408: break;
2409: case 2:
2410: switch (numCorners) {
2411: case 3: /* triangle */
2412: *numFaceVertices = 2; /* Edge has 2 vertices */
2413: break;
2414: case 4: /* quadrilateral */
2415: *numFaceVertices = 2; /* Edge has 2 vertices */
2416: break;
2417: case 6: /* quadratic triangle, tri and quad cohesive Lagrange cells */
2418: *numFaceVertices = 3; /* Edge has 3 vertices */
2419: break;
2420: case 9: /* quadratic quadrilateral, quadratic quad cohesive Lagrange cells */
2421: *numFaceVertices = 3; /* Edge has 3 vertices */
2422: break;
2423: default:
2424: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim);
2425: }
2426: break;
2427: case 3:
2428: switch (numCorners) {
2429: case 4: /* tetradehdron */
2430: *numFaceVertices = 3; /* Face has 3 vertices */
2431: break;
2432: case 6: /* tet cohesive cells */
2433: *numFaceVertices = 4; /* Face has 4 vertices */
2434: break;
2435: case 8: /* hexahedron */
2436: *numFaceVertices = 4; /* Face has 4 vertices */
2437: break;
2438: case 9: /* tet cohesive Lagrange cells */
2439: *numFaceVertices = 6; /* Face has 6 vertices */
2440: break;
2441: case 10: /* quadratic tetrahedron */
2442: *numFaceVertices = 6; /* Face has 6 vertices */
2443: break;
2444: case 12: /* hex cohesive Lagrange cells */
2445: *numFaceVertices = 6; /* Face has 6 vertices */
2446: break;
2447: case 18: /* quadratic tet cohesive Lagrange cells */
2448: *numFaceVertices = 6; /* Face has 6 vertices */
2449: break;
2450: case 27: /* quadratic hexahedron, quadratic hex cohesive Lagrange cells */
2451: *numFaceVertices = 9; /* Face has 9 vertices */
2452: break;
2453: default:
2454: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim);
2455: }
2456: break;
2457: default:
2458: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid cell dimension %d", cellDim);
2459: }
2460: return(0);
2461: }
2465: /*@
2466: DMPlexLocalizeCoordinate - If a mesh is periodic (a torus with lengths L_i, some of which can be infinite), project the coordinate onto [0, L_i) in each dimension.
2468: Input Parameters:
2469: + dm - The DM
2470: - in - The input coordinate point (dim numbers)
2472: Output Parameter:
2473: . out - The localized coordinate point
2475: Level: developer
2477: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeAddCoordinate()
2478: @*/
2479: PetscErrorCode DMPlexLocalizeCoordinate(DM dm, const PetscScalar in[], PetscScalar out[])
2480: {
2481: PetscInt dim, d;
2485: DMGetCoordinateDim(dm, &dim);
2486: if (!dm->maxCell) {
2487: for (d = 0; d < dim; ++d) out[d] = in[d];
2488: } else {
2489: for (d = 0; d < dim; ++d) {
2490: out[d] = in[d] - dm->L[d]*floor(PetscRealPart(in[d])/dm->L[d]);
2491: }
2492: }
2493: return(0);
2494: }
2498: /*
2499: DMPlexLocalizeCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.
2501: Input Parameters:
2502: + dm - The DM
2503: . dim - The spatial dimension
2504: . anchor - The anchor point, the input point can be no more than maxCell away from it
2505: - in - The input coordinate point (dim numbers)
2507: Output Parameter:
2508: . out - The localized coordinate point
2510: Level: developer
2512: Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell
2514: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeAddCoordinate()
2515: */
2516: PetscErrorCode DMPlexLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2517: {
2518: PetscInt d;
2521: if (!dm->maxCell) {
2522: for (d = 0; d < dim; ++d) out[d] = in[d];
2523: } else {
2524: for (d = 0; d < dim; ++d) {
2525: if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2526: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2527: } else {
2528: out[d] = in[d];
2529: }
2530: }
2531: }
2532: return(0);
2533: }
2536: PetscErrorCode DMPlexLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
2537: {
2538: PetscInt d;
2541: if (!dm->maxCell) {
2542: for (d = 0; d < dim; ++d) out[d] = in[d];
2543: } else {
2544: for (d = 0; d < dim; ++d) {
2545: if (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d]) {
2546: out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
2547: } else {
2548: out[d] = in[d];
2549: }
2550: }
2551: }
2552: return(0);
2553: }
2557: /*
2558: DMPlexLocalizeAddCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.
2560: Input Parameters:
2561: + dm - The DM
2562: . dim - The spatial dimension
2563: . anchor - The anchor point, the input point can be no more than maxCell away from it
2564: . in - The input coordinate delta (dim numbers)
2565: - out - The input coordinate point (dim numbers)
2567: Output Parameter:
2568: . out - The localized coordinate in + out
2570: Level: developer
2572: Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell
2574: .seealso: DMPlexLocalizeCoordinates(), DMPlexLocalizeCoordinate()
2575: */
2576: PetscErrorCode DMPlexLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
2577: {
2578: PetscInt d;
2581: if (!dm->maxCell) {
2582: for (d = 0; d < dim; ++d) out[d] += in[d];
2583: } else {
2584: for (d = 0; d < dim; ++d) {
2585: if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
2586: out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
2587: } else {
2588: out[d] += in[d];
2589: }
2590: }
2591: }
2592: return(0);
2593: }
2597: /*@
2598: DMPlexLocalizeCoordinates - If a mesh is periodic, create local coordinates for each cell
2600: Input Parameter:
2601: . dm - The DM
2603: Level: developer
2605: .seealso: DMPlexLocalizeCoordinate(), DMPlexLocalizeAddCoordinate()
2606: @*/
2607: PetscErrorCode DMPlexLocalizeCoordinates(DM dm)
2608: {
2609: PetscSection coordSection, cSection;
2610: Vec coordinates, cVec;
2611: PetscScalar *coords, *coords2, *anchor;
2612: PetscInt Nc, cStart, cEnd, c, vStart, vEnd, v, dof, d, off, off2, bs, coordSize;
2617: if (!dm->maxCell) return(0);
2618: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2619: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2620: DMGetCoordinatesLocal(dm, &coordinates);
2621: DMGetCoordinateSection(dm, &coordSection);
2622: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
2623: PetscSectionSetNumFields(cSection, 1);
2624: PetscSectionGetFieldComponents(coordSection, 0, &Nc);
2625: PetscSectionSetFieldComponents(cSection, 0, Nc);
2626: PetscSectionSetChart(cSection, cStart, vEnd);
2627: for (v = vStart; v < vEnd; ++v) {
2628: PetscSectionGetDof(coordSection, v, &dof);
2629: PetscSectionSetDof(cSection, v, dof);
2630: PetscSectionSetFieldDof(cSection, v, 0, dof);
2631: }
2632: for (c = cStart; c < cEnd; ++c) {
2633: DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, NULL);
2634: PetscSectionSetDof(cSection, c, dof);
2635: PetscSectionSetFieldDof(cSection, c, 0, dof);
2636: }
2637: PetscSectionSetUp(cSection);
2638: PetscSectionGetStorageSize(cSection, &coordSize);
2639: VecCreate(PetscObjectComm((PetscObject) dm), &cVec);
2640: VecGetBlockSize(coordinates, &bs);
2641: VecSetBlockSize(cVec, bs);
2642: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
2643: VecSetType(cVec,VECSTANDARD);
2644: VecGetArray(coordinates, &coords);
2645: VecGetArray(cVec, &coords2);
2646: for (v = vStart; v < vEnd; ++v) {
2647: PetscSectionGetDof(coordSection, v, &dof);
2648: PetscSectionGetOffset(coordSection, v, &off);
2649: PetscSectionGetOffset(cSection, v, &off2);
2650: for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
2651: }
2652: DMGetWorkArray(dm, 3, PETSC_SCALAR, &anchor);
2653: for (c = cStart; c < cEnd; ++c) {
2654: PetscScalar *cellCoords = NULL;
2655: PetscInt b;
2657: DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2658: PetscSectionGetOffset(cSection, c, &off2);
2659: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
2660: for (d = 0; d < dof/bs; ++d) {DMPlexLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
2661: DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
2662: }
2663: DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);
2664: VecRestoreArray(coordinates, &coords);
2665: VecRestoreArray(cVec, &coords2);
2666: DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
2667: DMSetCoordinatesLocal(dm, cVec);
2668: VecDestroy(&cVec);
2669: PetscSectionDestroy(&cSection);
2670: return(0);
2671: }
2675: /*@
2676: DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point
2678: Not Collective
2680: Input Parameter:
2681: . dm - The DMPlex object
2683: Output Parameter:
2684: . depthLabel - The DMLabel recording point depth
2686: Level: developer
2688: .keywords: mesh, points
2689: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2690: @*/
2691: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
2692: {
2693: DM_Plex *mesh = (DM_Plex*) dm->data;
2699: if (!mesh->depthLabel) {DMPlexGetLabel(dm, "depth", &mesh->depthLabel);}
2700: *depthLabel = mesh->depthLabel;
2701: return(0);
2702: }
2706: /*@
2707: DMPlexGetDepth - Get the depth of the DAG representing this mesh
2709: Not Collective
2711: Input Parameter:
2712: . dm - The DMPlex object
2714: Output Parameter:
2715: . depth - The number of strata (breadth first levels) in the DAG
2717: Level: developer
2719: .keywords: mesh, points
2720: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
2721: @*/
2722: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
2723: {
2724: DMLabel label;
2725: PetscInt d = 0;
2731: DMPlexGetDepthLabel(dm, &label);
2732: if (label) {DMLabelGetNumValues(label, &d);}
2733: *depth = d-1;
2734: return(0);
2735: }
2739: /*@
2740: DMPlexGetDepthStratum - Get the bounds [start, end) for all points at a certain depth.
2742: Not Collective
2744: Input Parameters:
2745: + dm - The DMPlex object
2746: - stratumValue - The requested depth
2748: Output Parameters:
2749: + start - The first point at this depth
2750: - end - One beyond the last point at this depth
2752: Level: developer
2754: .keywords: mesh, points
2755: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
2756: @*/
2757: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2758: {
2759: DMLabel label;
2760: PetscInt pStart, pEnd;
2767: DMPlexGetChart(dm, &pStart, &pEnd);
2768: if (pStart == pEnd) return(0);
2769: if (stratumValue < 0) {
2770: if (start) *start = pStart;
2771: if (end) *end = pEnd;
2772: return(0);
2773: }
2774: DMPlexGetDepthLabel(dm, &label);
2775: if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2776: DMLabelGetStratumBounds(label, stratumValue, start, end);
2777: return(0);
2778: }
2782: /*@
2783: DMPlexGetHeightStratum - Get the bounds [start, end) for all points at a certain height.
2785: Not Collective
2787: Input Parameters:
2788: + dm - The DMPlex object
2789: - stratumValue - The requested height
2791: Output Parameters:
2792: + start - The first point at this height
2793: - end - One beyond the last point at this height
2795: Level: developer
2797: .keywords: mesh, points
2798: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
2799: @*/
2800: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
2801: {
2802: DMLabel label;
2803: PetscInt depth, pStart, pEnd;
2810: DMPlexGetChart(dm, &pStart, &pEnd);
2811: if (pStart == pEnd) return(0);
2812: if (stratumValue < 0) {
2813: if (start) *start = pStart;
2814: if (end) *end = pEnd;
2815: return(0);
2816: }
2817: DMPlexGetDepthLabel(dm, &label);
2818: if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
2819: DMLabelGetNumValues(label, &depth);
2820: DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
2821: return(0);
2822: }
2826: /* Set the number of dof on each point and separate by fields */
2827: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
2828: {
2829: PetscInt *pMax;
2830: PetscInt depth, pStart = 0, pEnd = 0;
2831: PetscInt Nf, p, d, dep, f;
2832: PetscBool *isFE;
2836: PetscMalloc1(numFields, &isFE);
2837: DMGetNumFields(dm, &Nf);
2838: for (f = 0; f < numFields; ++f) {
2839: PetscObject obj;
2840: PetscClassId id;
2842: isFE[f] = PETSC_FALSE;
2843: if (f >= Nf) continue;
2844: DMGetField(dm, f, &obj);
2845: PetscObjectGetClassId(obj, &id);
2846: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
2847: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
2848: }
2849: PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
2850: if (numFields > 0) {
2851: PetscSectionSetNumFields(*section, numFields);
2852: if (numComp) {
2853: for (f = 0; f < numFields; ++f) {
2854: PetscSectionSetFieldComponents(*section, f, numComp[f]);
2855: }
2856: }
2857: }
2858: DMPlexGetChart(dm, &pStart, &pEnd);
2859: PetscSectionSetChart(*section, pStart, pEnd);
2860: DMPlexGetDepth(dm, &depth);
2861: PetscMalloc1(depth+1,&pMax);
2862: DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
2863: for (dep = 0; dep <= depth; ++dep) {
2864: d = dim == depth ? dep : (!dep ? 0 : dim);
2865: DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
2866: pMax[dep] = pMax[dep] < 0 ? pEnd : pMax[dep];
2867: for (p = pStart; p < pEnd; ++p) {
2868: PetscInt tot = 0;
2870: for (f = 0; f < numFields; ++f) {
2871: if (isFE[f] && p >= pMax[dep]) continue;
2872: PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
2873: tot += numDof[f*(dim+1)+d];
2874: }
2875: PetscSectionSetDof(*section, p, tot);
2876: }
2877: }
2878: PetscFree(pMax);
2879: PetscFree(isFE);
2880: return(0);
2881: }
2885: /* Set the number of dof on each point and separate by fields
2886: If bcComps is NULL or the IS is NULL, constrain every dof on the point
2887: */
2888: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2889: {
2890: PetscInt numFields;
2891: PetscInt bc;
2892: PetscSection aSec;
2896: PetscSectionGetNumFields(section, &numFields);
2897: for (bc = 0; bc < numBC; ++bc) {
2898: PetscInt field = 0;
2899: const PetscInt *comp;
2900: const PetscInt *idx;
2901: PetscInt Nc = -1, n, i;
2903: if (numFields) field = bcField[bc];
2904: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2905: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2906: ISGetLocalSize(bcPoints[bc], &n);
2907: ISGetIndices(bcPoints[bc], &idx);
2908: for (i = 0; i < n; ++i) {
2909: const PetscInt p = idx[i];
2910: PetscInt numConst;
2912: if (numFields) {
2913: PetscSectionGetFieldDof(section, p, field, &numConst);
2914: } else {
2915: PetscSectionGetDof(section, p, &numConst);
2916: }
2917: /* If Nc < 0, constrain every dof on the point */
2918: if (Nc > 0) numConst = PetscMin(numConst, Nc);
2919: if (numFields) {PetscSectionAddFieldConstraintDof(section, p, field, numConst);}
2920: PetscSectionAddConstraintDof(section, p, numConst);
2921: }
2922: ISRestoreIndices(bcPoints[bc], &idx);
2923: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2924: }
2925: DMPlexGetAnchors(dm, &aSec, NULL);
2926: if (aSec) {
2927: PetscInt aStart, aEnd, a;
2929: PetscSectionGetChart(aSec, &aStart, &aEnd);
2930: for (a = aStart; a < aEnd; a++) {
2931: PetscInt dof, f;
2933: PetscSectionGetDof(aSec, a, &dof);
2934: if (dof) {
2935: /* if there are point-to-point constraints, then all dofs are constrained */
2936: PetscSectionGetDof(section, a, &dof);
2937: PetscSectionSetConstraintDof(section, a, dof);
2938: for (f = 0; f < numFields; f++) {
2939: PetscSectionGetFieldDof(section, a, f, &dof);
2940: PetscSectionSetFieldConstraintDof(section, a, f, dof);
2941: }
2942: }
2943: }
2944: }
2945: return(0);
2946: }
2950: /* Set the constrained field indices on each point
2951: If bcComps is NULL or the IS is NULL, constrain every dof on the point
2952: */
2953: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section)
2954: {
2955: PetscSection aSec;
2956: PetscInt *indices;
2957: PetscInt numFields, maxDof, pStart, pEnd, p, bc, f, d;
2961: PetscSectionGetNumFields(section, &numFields);
2962: if (!numFields) return(0);
2963: /* Initialize all field indices to -1 */
2964: PetscSectionGetChart(section, &pStart, &pEnd);
2965: PetscSectionGetMaxDof(section, &maxDof);
2966: PetscMalloc1(maxDof, &indices);
2967: for (d = 0; d < maxDof; ++d) indices[d] = -1;
2968: for (p = pStart; p < pEnd; ++p) for (f = 0; f < numFields; ++f) {PetscSectionSetFieldConstraintIndices(section, p, f, indices);}
2969: /* Handle BC constraints */
2970: for (bc = 0; bc < numBC; ++bc) {
2971: const PetscInt field = bcField[bc];
2972: const PetscInt *comp, *idx;
2973: PetscInt Nc = -1, n, i;
2975: if (bcComps && bcComps[bc]) {ISGetLocalSize(bcComps[bc], &Nc);}
2976: if (bcComps && bcComps[bc]) {ISGetIndices(bcComps[bc], &comp);}
2977: ISGetLocalSize(bcPoints[bc], &n);
2978: ISGetIndices(bcPoints[bc], &idx);
2979: for (i = 0; i < n; ++i) {
2980: const PetscInt p = idx[i];
2981: const PetscInt *find;
2982: PetscInt fcdof, c;
2984: PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);
2985: if (Nc < 0) {
2986: for (d = 0; d < fcdof; ++d) indices[d] = d;
2987: } else {
2988: PetscSectionGetFieldConstraintIndices(section, p, field, &find);
2989: for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];}
2990: for (c = 0; c < Nc; ++c) indices[d+c] = comp[c];
2991: PetscSortInt(d+Nc, indices);
2992: for (c = d+Nc; c < fcdof; ++c) indices[c] = -1;
2993: }
2994: PetscSectionSetFieldConstraintIndices(section, p, field, indices);
2995: }
2996: if (bcComps && bcComps[bc]) {ISRestoreIndices(bcComps[bc], &comp);}
2997: ISRestoreIndices(bcPoints[bc], &idx);
2998: }
2999: /* Handle anchors */
3000: DMPlexGetAnchors(dm, &aSec, NULL);
3001: if (aSec) {
3002: PetscInt aStart, aEnd, a;
3004: for (d = 0; d < maxDof; ++d) indices[d] = d;
3005: PetscSectionGetChart(aSec, &aStart, &aEnd);
3006: for (a = aStart; a < aEnd; a++) {
3007: PetscInt dof, fdof, f;
3009: PetscSectionGetDof(aSec, a, &dof);
3010: if (dof) {
3011: /* if there are point-to-point constraints, then all dofs are constrained */
3012: for (f = 0; f < numFields; f++) {
3013: PetscSectionGetFieldDof(section, a, f, &fdof);
3014: PetscSectionSetFieldConstraintIndices(section, a, f, indices);
3015: }
3016: }
3017: }
3018: }
3019: PetscFree(indices);
3020: return(0);
3021: }
3025: /* Set the constrained indices on each point */
3026: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
3027: {
3028: PetscInt *indices;
3029: PetscInt numFields, maxDof, pStart, pEnd, p, f, d;
3033: PetscSectionGetNumFields(section, &numFields);
3034: PetscSectionGetMaxDof(section, &maxDof);
3035: PetscSectionGetChart(section, &pStart, &pEnd);
3036: PetscMalloc1(maxDof, &indices);
3037: for (d = 0; d < maxDof; ++d) indices[d] = -1;
3038: for (p = pStart; p < pEnd; ++p) {
3039: PetscInt cdof, d;
3041: PetscSectionGetConstraintDof(section, p, &cdof);
3042: if (cdof) {
3043: if (numFields) {
3044: PetscInt numConst = 0, foff = 0;
3046: for (f = 0; f < numFields; ++f) {
3047: const PetscInt *find;
3048: PetscInt fcdof, fdof;
3050: PetscSectionGetFieldDof(section, p, f, &fdof);
3051: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
3052: /* Change constraint numbering from field component to local dof number */
3053: PetscSectionGetFieldConstraintIndices(section, p, f, &find);
3054: for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff;
3055: numConst += fcdof;
3056: foff += fdof;
3057: }
3058: if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
3059: } else {
3060: for (d = 0; d < cdof; ++d) indices[d] = d;
3061: }
3062: PetscSectionSetConstraintIndices(section, p, indices);
3063: }
3064: }
3065: PetscFree(indices);
3066: return(0);
3067: }
3071: /*@C
3072: DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided.
3074: Not Collective
3076: Input Parameters:
3077: + dm - The DMPlex object
3078: . dim - The spatial dimension of the problem
3079: . numFields - The number of fields in the problem
3080: . numComp - An array of size numFields that holds the number of components for each field
3081: . numDof - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
3082: . numBC - The number of boundary conditions
3083: . bcField - An array of size numBC giving the field number for each boundry condition
3084: . bcComps - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies
3085: . bcPoints - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies
3086: - perm - Optional permutation of the chart, or NULL
3088: Output Parameter:
3089: . section - The PetscSection object
3091: Notes: numDof[f*(dim+1)+d] gives the number of dof for field f on sieve points of dimension d. For instance, numDof[1] is the
3092: number of dof for field 0 on each edge.
3094: The chart permutation is the same one set using PetscSectionSetPermutation()
3096: Level: developer
3098: Fortran Notes:
3099: A Fortran 90 version is available as DMPlexCreateSectionF90()
3101: .keywords: mesh, elements
3102: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
3103: @*/
3104: PetscErrorCode DMPlexCreateSection(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], IS perm, PetscSection *section)
3105: {
3106: PetscSection aSec;
3110: DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
3111: DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);
3112: if (perm) {PetscSectionSetPermutation(*section, perm);}
3113: PetscSectionSetUp(*section);
3114: DMPlexGetAnchors(dm,&aSec,NULL);
3115: if (numBC || aSec) {
3116: DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);
3117: DMPlexCreateSectionBCIndices(dm, *section);
3118: }
3119: PetscSectionViewFromOptions(*section,NULL,"-section_view");
3120: return(0);
3121: }
3125: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
3126: {
3127: PetscSection section;
3131: DMClone(dm, cdm);
3132: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
3133: DMSetDefaultSection(*cdm, section);
3134: PetscSectionDestroy(§ion);
3135: return(0);
3136: }
3140: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
3141: {
3142: DM_Plex *mesh = (DM_Plex*) dm->data;
3146: if (section) *section = mesh->coneSection;
3147: return(0);
3148: }
3152: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
3153: {
3154: DM_Plex *mesh = (DM_Plex*) dm->data;
3158: if (section) *section = mesh->supportSection;
3159: return(0);
3160: }
3164: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
3165: {
3166: DM_Plex *mesh = (DM_Plex*) dm->data;
3170: if (cones) *cones = mesh->cones;
3171: return(0);
3172: }
3176: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
3177: {
3178: DM_Plex *mesh = (DM_Plex*) dm->data;
3182: if (coneOrientations) *coneOrientations = mesh->coneOrientations;
3183: return(0);
3184: }
3186: /******************************** FEM Support **********************************/
3190: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3191: {
3192: PetscScalar *array, *vArray;
3193: const PetscInt *cone, *coneO;
3194: PetscInt pStart, pEnd, p, numPoints, size = 0, offset = 0;
3195: PetscErrorCode ierr;
3198: PetscSectionGetChart(section, &pStart, &pEnd);
3199: DMPlexGetConeSize(dm, point, &numPoints);
3200: DMPlexGetCone(dm, point, &cone);
3201: DMPlexGetConeOrientation(dm, point, &coneO);
3202: if (!values || !*values) {
3203: if ((point >= pStart) && (point < pEnd)) {
3204: PetscInt dof;
3206: PetscSectionGetDof(section, point, &dof);
3207: size += dof;
3208: }
3209: for (p = 0; p < numPoints; ++p) {
3210: const PetscInt cp = cone[p];
3211: PetscInt dof;
3213: if ((cp < pStart) || (cp >= pEnd)) continue;
3214: PetscSectionGetDof(section, cp, &dof);
3215: size += dof;
3216: }
3217: if (!values) {
3218: if (csize) *csize = size;
3219: return(0);
3220: }
3221: DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
3222: } else {
3223: array = *values;
3224: }
3225: size = 0;
3226: VecGetArray(v, &vArray);
3227: if ((point >= pStart) && (point < pEnd)) {
3228: PetscInt dof, off, d;
3229: PetscScalar *varr;
3231: PetscSectionGetDof(section, point, &dof);
3232: PetscSectionGetOffset(section, point, &off);
3233: varr = &vArray[off];
3234: for (d = 0; d < dof; ++d, ++offset) {
3235: array[offset] = varr[d];
3236: }
3237: size += dof;
3238: }
3239: for (p = 0; p < numPoints; ++p) {
3240: const PetscInt cp = cone[p];
3241: PetscInt o = coneO[p];
3242: PetscInt dof, off, d;
3243: PetscScalar *varr;
3245: if ((cp < pStart) || (cp >= pEnd)) continue;
3246: PetscSectionGetDof(section, cp, &dof);
3247: PetscSectionGetOffset(section, cp, &off);
3248: varr = &vArray[off];
3249: if (o >= 0) {
3250: for (d = 0; d < dof; ++d, ++offset) {
3251: array[offset] = varr[d];
3252: }
3253: } else {
3254: for (d = dof-1; d >= 0; --d, ++offset) {
3255: array[offset] = varr[d];
3256: }
3257: }
3258: size += dof;
3259: }
3260: VecRestoreArray(v, &vArray);
3261: if (!*values) {
3262: if (csize) *csize = size;
3263: *values = array;
3264: } else {
3265: if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3266: *csize = size;
3267: }
3268: return(0);
3269: }
3273: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3274: {
3275: PetscInt offset = 0, p;
3279: *size = 0;
3280: for (p = 0; p < numPoints*2; p += 2) {
3281: const PetscInt point = points[p];
3282: const PetscInt o = points[p+1];
3283: PetscInt dof, off, d;
3284: const PetscScalar *varr;
3286: PetscSectionGetDof(section, point, &dof);
3287: PetscSectionGetOffset(section, point, &off);
3288: varr = &vArray[off];
3289: if (o >= 0) {
3290: for (d = 0; d < dof; ++d, ++offset) array[offset] = varr[d];
3291: } else {
3292: for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
3293: }
3294: }
3295: *size = offset;
3296: return(0);
3297: }
3301: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Fields_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], PetscInt numFields, const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
3302: {
3303: PetscInt offset = 0, f;
3307: *size = 0;
3308: for (f = 0; f < numFields; ++f) {
3309: PetscInt fcomp, p;
3311: PetscSectionGetFieldComponents(section, f, &fcomp);
3312: for (p = 0; p < numPoints*2; p += 2) {
3313: const PetscInt point = points[p];
3314: const PetscInt o = points[p+1];
3315: PetscInt fdof, foff, d, c;
3316: const PetscScalar *varr;
3318: PetscSectionGetFieldDof(section, point, f, &fdof);
3319: PetscSectionGetFieldOffset(section, point, f, &foff);
3320: varr = &vArray[foff];
3321: if (o >= 0) {
3322: for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
3323: } else {
3324: for (d = fdof/fcomp-1; d >= 0; --d) {
3325: for (c = 0; c < fcomp; ++c, ++offset) {
3326: array[offset] = varr[d*fcomp+c];
3327: }
3328: }
3329: }
3330: }
3331: }
3332: *size = offset;
3333: return(0);
3334: }
3338: /*@C
3339: DMPlexVecGetClosure - Get an array of the values on the closure of 'point'
3341: Not collective
3343: Input Parameters:
3344: + dm - The DM
3345: . section - The section describing the layout in v, or NULL to use the default section
3346: . v - The local vector
3347: - point - The sieve point in the DM
3349: Output Parameters:
3350: + csize - The number of values in the closure, or NULL
3351: - values - The array of values, which is a borrowed array and should not be freed
3353: Fortran Notes:
3354: Since it returns an array, this routine is only available in Fortran 90, and you must
3355: include petsc.h90 in your code.
3357: The csize argument is not present in the Fortran 90 binding since it is internal to the array.
3359: Level: intermediate
3361: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3362: @*/
3363: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3364: {
3365: PetscSection clSection;
3366: IS clPoints;
3367: PetscScalar *array, *vArray;
3368: PetscInt *points = NULL;
3369: const PetscInt *clp;
3370: PetscInt depth, numFields, numPoints, size;
3371: PetscErrorCode ierr;
3375: if (!section) {DMGetDefaultSection(dm, §ion);}
3378: DMPlexGetDepth(dm, &depth);
3379: PetscSectionGetNumFields(section, &numFields);
3380: if (depth == 1 && numFields < 2) {
3381: DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
3382: return(0);
3383: }
3384: /* Get points */
3385: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3386: if (!clPoints) {
3387: PetscInt pStart, pEnd, p, q;
3389: PetscSectionGetChart(section, &pStart, &pEnd);
3390: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3391: /* Compress out points not in the section */
3392: for (p = 0, q = 0; p < numPoints*2; p += 2) {
3393: if ((points[p] >= pStart) && (points[p] < pEnd)) {
3394: points[q*2] = points[p];
3395: points[q*2+1] = points[p+1];
3396: ++q;
3397: }
3398: }
3399: numPoints = q;
3400: } else {
3401: PetscInt dof, off;
3403: PetscSectionGetDof(clSection, point, &dof);
3404: PetscSectionGetOffset(clSection, point, &off);
3405: ISGetIndices(clPoints, &clp);
3406: numPoints = dof/2;
3407: points = (PetscInt *) &clp[off];
3408: }
3409: /* Get array */
3410: if (!values || !*values) {
3411: PetscInt asize = 0, dof, p;
3413: for (p = 0; p < numPoints*2; p += 2) {
3414: PetscSectionGetDof(section, points[p], &dof);
3415: asize += dof;
3416: }
3417: if (!values) {
3418: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3419: else {ISRestoreIndices(clPoints, &clp);}
3420: if (csize) *csize = asize;
3421: return(0);
3422: }
3423: DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
3424: } else {
3425: array = *values;
3426: }
3427: VecGetArray(v, &vArray);
3428: /* Get values */
3429: if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);}
3430: else {DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);}
3431: /* Cleanup points */
3432: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3433: else {ISRestoreIndices(clPoints, &clp);}
3434: /* Cleanup array */
3435: VecRestoreArray(v, &vArray);
3436: if (!*values) {
3437: if (csize) *csize = size;
3438: *values = array;
3439: } else {
3440: if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
3441: *csize = size;
3442: }
3443: return(0);
3444: }
3448: /*@C
3449: DMPlexVecRestoreClosure - Restore the array of the values on the closure of 'point'
3451: Not collective
3453: Input Parameters:
3454: + dm - The DM
3455: . section - The section describing the layout in v, or NULL to use the default section
3456: . v - The local vector
3457: . point - The sieve point in the DM
3458: . csize - The number of values in the closure, or NULL
3459: - values - The array of values, which is a borrowed array and should not be freed
3461: Fortran Notes:
3462: Since it returns an array, this routine is only available in Fortran 90, and you must
3463: include petsc.h90 in your code.
3465: The csize argument is not present in the Fortran 90 binding since it is internal to the array.
3467: Level: intermediate
3469: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
3470: @*/
3471: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
3472: {
3473: PetscInt size = 0;
3477: /* Should work without recalculating size */
3478: DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
3479: return(0);
3480: }
3482: PETSC_STATIC_INLINE void add (PetscScalar *x, PetscScalar y) {*x += y;}
3483: PETSC_STATIC_INLINE void insert(PetscScalar *x, PetscScalar y) {*x = y;}
3487: PETSC_STATIC_INLINE PetscErrorCode updatePoint_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3488: {
3489: PetscInt cdof; /* The number of constraints on this point */
3490: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3491: PetscScalar *a;
3492: PetscInt off, cind = 0, k;
3493: PetscErrorCode ierr;
3496: PetscSectionGetConstraintDof(section, point, &cdof);
3497: PetscSectionGetOffset(section, point, &off);
3498: a = &array[off];
3499: if (!cdof || setBC) {
3500: if (orientation >= 0) {
3501: for (k = 0; k < dof; ++k) {
3502: fuse(&a[k], values[k]);
3503: }
3504: } else {
3505: for (k = 0; k < dof; ++k) {
3506: fuse(&a[k], values[dof-k-1]);
3507: }
3508: }
3509: } else {
3510: PetscSectionGetConstraintIndices(section, point, &cdofs);
3511: if (orientation >= 0) {
3512: for (k = 0; k < dof; ++k) {
3513: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3514: fuse(&a[k], values[k]);
3515: }
3516: } else {
3517: for (k = 0; k < dof; ++k) {
3518: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3519: fuse(&a[k], values[dof-k-1]);
3520: }
3521: }
3522: }
3523: return(0);
3524: }
3528: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
3529: {
3530: PetscInt cdof; /* The number of constraints on this point */
3531: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3532: PetscScalar *a;
3533: PetscInt off, cind = 0, k;
3534: PetscErrorCode ierr;
3537: PetscSectionGetConstraintDof(section, point, &cdof);
3538: PetscSectionGetOffset(section, point, &off);
3539: a = &array[off];
3540: if (cdof) {
3541: PetscSectionGetConstraintIndices(section, point, &cdofs);
3542: if (orientation >= 0) {
3543: for (k = 0; k < dof; ++k) {
3544: if ((cind < cdof) && (k == cdofs[cind])) {
3545: fuse(&a[k], values[k]);
3546: ++cind;
3547: }
3548: }
3549: } else {
3550: for (k = 0; k < dof; ++k) {
3551: if ((cind < cdof) && (k == cdofs[cind])) {
3552: fuse(&a[k], values[dof-k-1]);
3553: ++cind;
3554: }
3555: }
3556: }
3557: }
3558: return(0);
3559: }
3563: PETSC_STATIC_INLINE PetscErrorCode updatePointFields_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, const PetscScalar values[], PetscInt *offset, PetscScalar array[])
3564: {
3565: PetscScalar *a;
3566: PetscInt fdof, foff, fcdof, foffset = *offset;
3567: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3568: PetscInt cind = 0, k, c;
3569: PetscErrorCode ierr;
3572: PetscSectionGetFieldDof(section, point, f, &fdof);
3573: PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3574: PetscSectionGetFieldOffset(section, point, f, &foff);
3575: a = &array[foff];
3576: if (!fcdof || setBC) {
3577: if (o >= 0) {
3578: for (k = 0; k < fdof; ++k) fuse(&a[k], values[foffset+k]);
3579: } else {
3580: for (k = fdof/fcomp-1; k >= 0; --k) {
3581: for (c = 0; c < fcomp; ++c) {
3582: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3583: }
3584: }
3585: }
3586: } else {
3587: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3588: if (o >= 0) {
3589: for (k = 0; k < fdof; ++k) {
3590: if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
3591: fuse(&a[k], values[foffset+k]);
3592: }
3593: } else {
3594: for (k = fdof/fcomp-1; k >= 0; --k) {
3595: for (c = 0; c < fcomp; ++c) {
3596: if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
3597: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3598: }
3599: }
3600: }
3601: }
3602: *offset += fdof;
3603: return(0);
3604: }
3608: PETSC_STATIC_INLINE PetscErrorCode updatePointFieldsBC_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), const PetscScalar values[], PetscInt *offset, PetscScalar array[])
3609: {
3610: PetscScalar *a;
3611: PetscInt fdof, foff, fcdof, foffset = *offset;
3612: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
3613: PetscInt cind = 0, k, c;
3614: PetscErrorCode ierr;
3617: PetscSectionGetFieldDof(section, point, f, &fdof);
3618: PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
3619: PetscSectionGetFieldOffset(section, point, f, &foff);
3620: a = &array[foff];
3621: if (fcdof) {
3622: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
3623: if (o >= 0) {
3624: for (k = 0; k < fdof; ++k) {
3625: if ((cind < fcdof) && (k == fcdofs[cind])) {
3626: fuse(&a[k], values[foffset+k]);
3627: ++cind;
3628: }
3629: }
3630: } else {
3631: for (k = fdof/fcomp-1; k >= 0; --k) {
3632: for (c = 0; c < fcomp; ++c) {
3633: if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {
3634: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
3635: ++cind;
3636: }
3637: }
3638: }
3639: }
3640: }
3641: *offset += fdof;
3642: return(0);
3643: }
3647: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3648: {
3649: PetscScalar *array;
3650: const PetscInt *cone, *coneO;
3651: PetscInt pStart, pEnd, p, numPoints, off, dof;
3652: PetscErrorCode ierr;
3655: PetscSectionGetChart(section, &pStart, &pEnd);
3656: DMPlexGetConeSize(dm, point, &numPoints);
3657: DMPlexGetCone(dm, point, &cone);
3658: DMPlexGetConeOrientation(dm, point, &coneO);
3659: VecGetArray(v, &array);
3660: for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
3661: const PetscInt cp = !p ? point : cone[p-1];
3662: const PetscInt o = !p ? 0 : coneO[p-1];
3664: if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
3665: PetscSectionGetDof(section, cp, &dof);
3666: /* ADD_VALUES */
3667: {
3668: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3669: PetscScalar *a;
3670: PetscInt cdof, coff, cind = 0, k;
3672: PetscSectionGetConstraintDof(section, cp, &cdof);
3673: PetscSectionGetOffset(section, cp, &coff);
3674: a = &array[coff];
3675: if (!cdof) {
3676: if (o >= 0) {
3677: for (k = 0; k < dof; ++k) {
3678: a[k] += values[off+k];
3679: }
3680: } else {
3681: for (k = 0; k < dof; ++k) {
3682: a[k] += values[off+dof-k-1];
3683: }
3684: }
3685: } else {
3686: PetscSectionGetConstraintIndices(section, cp, &cdofs);
3687: if (o >= 0) {
3688: for (k = 0; k < dof; ++k) {
3689: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3690: a[k] += values[off+k];
3691: }
3692: } else {
3693: for (k = 0; k < dof; ++k) {
3694: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
3695: a[k] += values[off+dof-k-1];
3696: }
3697: }
3698: }
3699: }
3700: }
3701: VecRestoreArray(v, &array);
3702: return(0);
3703: }
3707: /*@C
3708: DMPlexVecSetClosure - Set an array of the values on the closure of 'point'
3710: Not collective
3712: Input Parameters:
3713: + dm - The DM
3714: . section - The section describing the layout in v, or NULL to use the default section
3715: . v - The local vector
3716: . point - The sieve point in the DM
3717: . values - The array of values
3718: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
3720: Fortran Notes:
3721: This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
3723: Level: intermediate
3725: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
3726: @*/
3727: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
3728: {
3729: PetscSection clSection;
3730: IS clPoints;
3731: PetscScalar *array;
3732: PetscInt *points = NULL;
3733: const PetscInt *clp;
3734: PetscInt depth, numFields, numPoints, p;
3735: PetscErrorCode ierr;
3739: if (!section) {DMGetDefaultSection(dm, §ion);}
3742: DMPlexGetDepth(dm, &depth);
3743: PetscSectionGetNumFields(section, &numFields);
3744: if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
3745: DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
3746: return(0);
3747: }
3748: /* Get points */
3749: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3750: if (!clPoints) {
3751: PetscInt pStart, pEnd, q;
3753: PetscSectionGetChart(section, &pStart, &pEnd);
3754: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3755: /* Compress out points not in the section */
3756: for (p = 0, q = 0; p < numPoints*2; p += 2) {
3757: if ((points[p] >= pStart) && (points[p] < pEnd)) {
3758: points[q*2] = points[p];
3759: points[q*2+1] = points[p+1];
3760: ++q;
3761: }
3762: }
3763: numPoints = q;
3764: } else {
3765: PetscInt dof, off;
3767: PetscSectionGetDof(clSection, point, &dof);
3768: PetscSectionGetOffset(clSection, point, &off);
3769: ISGetIndices(clPoints, &clp);
3770: numPoints = dof/2;
3771: points = (PetscInt *) &clp[off];
3772: }
3773: /* Get array */
3774: VecGetArray(v, &array);
3775: /* Get values */
3776: if (numFields > 0) {
3777: PetscInt offset = 0, fcomp, f;
3778: for (f = 0; f < numFields; ++f) {
3779: PetscSectionGetFieldComponents(section, f, &fcomp);
3780: switch (mode) {
3781: case INSERT_VALUES:
3782: for (p = 0; p < numPoints*2; p += 2) {
3783: const PetscInt point = points[p];
3784: const PetscInt o = points[p+1];
3785: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3786: } break;
3787: case INSERT_ALL_VALUES:
3788: for (p = 0; p < numPoints*2; p += 2) {
3789: const PetscInt point = points[p];
3790: const PetscInt o = points[p+1];
3791: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3792: } break;
3793: case INSERT_BC_VALUES:
3794: for (p = 0; p < numPoints*2; p += 2) {
3795: const PetscInt point = points[p];
3796: const PetscInt o = points[p+1];
3797: updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3798: } break;
3799: case ADD_VALUES:
3800: for (p = 0; p < numPoints*2; p += 2) {
3801: const PetscInt point = points[p];
3802: const PetscInt o = points[p+1];
3803: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3804: } break;
3805: case ADD_ALL_VALUES:
3806: for (p = 0; p < numPoints*2; p += 2) {
3807: const PetscInt point = points[p];
3808: const PetscInt o = points[p+1];
3809: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3810: } break;
3811: default:
3812: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3813: }
3814: }
3815: } else {
3816: PetscInt dof, off;
3818: switch (mode) {
3819: case INSERT_VALUES:
3820: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3821: PetscInt o = points[p+1];
3822: PetscSectionGetDof(section, points[p], &dof);
3823: updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
3824: } break;
3825: case INSERT_ALL_VALUES:
3826: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3827: PetscInt o = points[p+1];
3828: PetscSectionGetDof(section, points[p], &dof);
3829: updatePoint_private(section, points[p], dof, insert, PETSC_TRUE, o, &values[off], array);
3830: } break;
3831: case INSERT_BC_VALUES:
3832: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3833: PetscInt o = points[p+1];
3834: PetscSectionGetDof(section, points[p], &dof);
3835: updatePointBC_private(section, points[p], dof, insert, o, &values[off], array);
3836: } break;
3837: case ADD_VALUES:
3838: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3839: PetscInt o = points[p+1];
3840: PetscSectionGetDof(section, points[p], &dof);
3841: updatePoint_private(section, points[p], dof, add, PETSC_FALSE, o, &values[off], array);
3842: } break;
3843: case ADD_ALL_VALUES:
3844: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
3845: PetscInt o = points[p+1];
3846: PetscSectionGetDof(section, points[p], &dof);
3847: updatePoint_private(section, points[p], dof, add, PETSC_TRUE, o, &values[off], array);
3848: } break;
3849: default:
3850: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3851: }
3852: }
3853: /* Cleanup points */
3854: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3855: else {ISRestoreIndices(clPoints, &clp);}
3856: /* Cleanup array */
3857: VecRestoreArray(v, &array);
3858: return(0);
3859: }
3863: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
3864: {
3865: PetscSection clSection;
3866: IS clPoints;
3867: PetscScalar *array;
3868: PetscInt *points = NULL;
3869: const PetscInt *clp;
3870: PetscInt numFields, numPoints, p;
3871: PetscInt offset = 0, fcomp, f;
3872: PetscErrorCode ierr;
3876: if (!section) {DMGetDefaultSection(dm, §ion);}
3879: PetscSectionGetNumFields(section, &numFields);
3880: /* Get points */
3881: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
3882: if (!clPoints) {
3883: PetscInt pStart, pEnd, q;
3885: PetscSectionGetChart(section, &pStart, &pEnd);
3886: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
3887: /* Compress out points not in the section */
3888: for (p = 0, q = 0; p < numPoints*2; p += 2) {
3889: if ((points[p] >= pStart) && (points[p] < pEnd)) {
3890: points[q*2] = points[p];
3891: points[q*2+1] = points[p+1];
3892: ++q;
3893: }
3894: }
3895: numPoints = q;
3896: } else {
3897: PetscInt dof, off;
3899: PetscSectionGetDof(clSection, point, &dof);
3900: PetscSectionGetOffset(clSection, point, &off);
3901: ISGetIndices(clPoints, &clp);
3902: numPoints = dof/2;
3903: points = (PetscInt *) &clp[off];
3904: }
3905: /* Get array */
3906: VecGetArray(v, &array);
3907: /* Get values */
3908: for (f = 0; f < numFields; ++f) {
3909: PetscSectionGetFieldComponents(section, f, &fcomp);
3910: if (!fieldActive[f]) {
3911: for (p = 0; p < numPoints*2; p += 2) {
3912: PetscInt fdof;
3913: PetscSectionGetFieldDof(section, points[p], f, &fdof);
3914: offset += fdof;
3915: }
3916: continue;
3917: }
3918: switch (mode) {
3919: case INSERT_VALUES:
3920: for (p = 0; p < numPoints*2; p += 2) {
3921: const PetscInt point = points[p];
3922: const PetscInt o = points[p+1];
3923: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
3924: } break;
3925: case INSERT_ALL_VALUES:
3926: for (p = 0; p < numPoints*2; p += 2) {
3927: const PetscInt point = points[p];
3928: const PetscInt o = points[p+1];
3929: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
3930: } break;
3931: case INSERT_BC_VALUES:
3932: for (p = 0; p < numPoints*2; p += 2) {
3933: const PetscInt point = points[p];
3934: const PetscInt o = points[p+1];
3935: updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
3936: } break;
3937: case ADD_VALUES:
3938: for (p = 0; p < numPoints*2; p += 2) {
3939: const PetscInt point = points[p];
3940: const PetscInt o = points[p+1];
3941: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
3942: } break;
3943: case ADD_ALL_VALUES:
3944: for (p = 0; p < numPoints*2; p += 2) {
3945: const PetscInt point = points[p];
3946: const PetscInt o = points[p+1];
3947: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
3948: } break;
3949: default:
3950: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
3951: }
3952: }
3953: /* Cleanup points */
3954: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
3955: else {ISRestoreIndices(clPoints, &clp);}
3956: /* Cleanup array */
3957: VecRestoreArray(v, &array);
3958: return(0);
3959: }
3963: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
3964: {
3965: PetscMPIInt rank;
3966: PetscInt i, j;
3970: MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
3971: PetscViewerASCIIPrintf(viewer, "[%D]mat for sieve point %D\n", rank, point);
3972: for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
3973: for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%d]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
3974: numCIndices = numCIndices ? numCIndices : numRIndices;
3975: for (i = 0; i < numRIndices; i++) {
3976: PetscViewerASCIIPrintf(viewer, "[%d]", rank);
3977: for (j = 0; j < numCIndices; j++) {
3978: #if defined(PETSC_USE_COMPLEX)
3979: PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
3980: #else
3981: PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
3982: #endif
3983: }
3984: PetscViewerASCIIPrintf(viewer, "\n");
3985: }
3986: return(0);
3987: }
3991: /* . off - The global offset of this point */
3992: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
3993: {
3994: PetscInt dof; /* The number of unknowns on this point */
3995: PetscInt cdof; /* The number of constraints on this point */
3996: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
3997: PetscInt cind = 0, k;
3998: PetscErrorCode ierr;
4001: PetscSectionGetDof(section, point, &dof);
4002: PetscSectionGetConstraintDof(section, point, &cdof);
4003: if (!cdof || setBC) {
4004: if (orientation >= 0) {
4005: for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
4006: } else {
4007: for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
4008: }
4009: } else {
4010: PetscSectionGetConstraintIndices(section, point, &cdofs);
4011: if (orientation >= 0) {
4012: for (k = 0; k < dof; ++k) {
4013: if ((cind < cdof) && (k == cdofs[cind])) {
4014: /* Insert check for returning constrained indices */
4015: indices[*loff+k] = -(off+k+1);
4016: ++cind;
4017: } else {
4018: indices[*loff+k] = off+k-cind;
4019: }
4020: }
4021: } else {
4022: for (k = 0; k < dof; ++k) {
4023: if ((cind < cdof) && (k == cdofs[cind])) {
4024: /* Insert check for returning constrained indices */
4025: indices[*loff+dof-k-1] = -(off+k+1);
4026: ++cind;
4027: } else {
4028: indices[*loff+dof-k-1] = off+k-cind;
4029: }
4030: }
4031: }
4032: }
4033: *loff += dof;
4034: return(0);
4035: }
4039: /* . off - The global offset of this point */
4040: PetscErrorCode indicesPointFields_private(PetscSection section, PetscInt point, PetscInt off, PetscInt foffs[], PetscBool setBC, PetscInt orientation, PetscInt indices[])
4041: {
4042: PetscInt numFields, foff, f;
4046: PetscSectionGetNumFields(section, &numFields);
4047: for (f = 0, foff = 0; f < numFields; ++f) {
4048: PetscInt fdof, fcomp, cfdof;
4049: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4050: PetscInt cind = 0, k, c;
4052: PetscSectionGetFieldComponents(section, f, &fcomp);
4053: PetscSectionGetFieldDof(section, point, f, &fdof);
4054: PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
4055: if (!cfdof || setBC) {
4056: if (orientation >= 0) {
4057: for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
4058: } else {
4059: for (k = fdof/fcomp-1; k >= 0; --k) {
4060: for (c = 0; c < fcomp; ++c) {
4061: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
4062: }
4063: }
4064: }
4065: } else {
4066: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
4067: if (orientation >= 0) {
4068: for (k = 0; k < fdof; ++k) {
4069: if ((cind < cfdof) && (k == fcdofs[cind])) {
4070: indices[foffs[f]+k] = -(off+foff+k+1);
4071: ++cind;
4072: } else {
4073: indices[foffs[f]+k] = off+foff+k-cind;
4074: }
4075: }
4076: } else {
4077: for (k = fdof/fcomp-1; k >= 0; --k) {
4078: for (c = 0; c < fcomp; ++c) {
4079: if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
4080: indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
4081: ++cind;
4082: } else {
4083: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
4084: }
4085: }
4086: }
4087: }
4088: }
4089: foff += fdof - cfdof;
4090: foffs[f] += fdof;
4091: }
4092: return(0);
4093: }
4097: static PetscErrorCode DMPlexAnchorsModifyMat(DM dm, PetscSection section, PetscInt numPoints, PetscInt numIndices, const PetscInt points[], const PetscScalar values[], PetscInt *outNumPoints, PetscInt *outNumIndices, PetscInt *outPoints[], PetscScalar *outValues[], PetscInt offsets[])
4098: {
4099: Mat cMat;
4100: PetscSection aSec, cSec;
4101: IS aIS;
4102: PetscInt aStart = -1, aEnd = -1;
4103: const PetscInt *anchors;
4104: PetscInt numFields, f, p, q, newP = 0;
4105: PetscInt newNumPoints = 0, newNumIndices = 0;
4106: PetscInt *newPoints, *indices, *newIndices;
4107: PetscInt maxAnchor, maxDof;
4108: PetscInt newOffsets[32];
4109: PetscInt *pointMatOffsets[32];
4110: PetscInt *newPointOffsets[32];
4111: PetscScalar *pointMat[32];
4112: PetscScalar *newValues,*tmpValues;
4113: PetscBool anyConstrained = PETSC_FALSE;
4114: PetscErrorCode ierr;
4119: PetscSectionGetNumFields(section, &numFields);
4121: DMPlexGetAnchors(dm,&aSec,&aIS);
4122: /* if there are point-to-point constraints */
4123: if (aSec) {
4124: PetscMemzero(newOffsets, 32 * sizeof(PetscInt));
4125: ISGetIndices(aIS,&anchors);
4126: PetscSectionGetChart(aSec,&aStart,&aEnd);
4127: /* figure out how many points are going to be in the new element matrix
4128: * (we allow double counting, because it's all just going to be summed
4129: * into the global matrix anyway) */
4130: for (p = 0; p < 2*numPoints; p+=2) {
4131: PetscInt b = points[p];
4132: PetscInt bDof = 0;
4134: if (b >= aStart && b < aEnd) {
4135: PetscSectionGetDof(aSec,b,&bDof);
4136: }
4137: if (bDof) {
4138: /* this point is constrained */
4139: /* it is going to be replaced by its anchors */
4140: PetscInt bOff, q;
4142: anyConstrained = PETSC_TRUE;
4143: newNumPoints += bDof;
4144: PetscSectionGetOffset(aSec,b,&bOff);
4145: for (q = 0; q < bDof; q++) {
4146: PetscInt a = anchors[bOff + q];
4147: PetscInt aDof;
4149: PetscSectionGetDof(section,a,&aDof);
4150: newNumIndices += aDof;
4151: for (f = 0; f < numFields; ++f) {
4152: PetscInt fDof;
4154: PetscSectionGetFieldDof(section, a, f, &fDof);
4155: newOffsets[f+1] += fDof;
4156: }
4157: }
4158: }
4159: else {
4160: /* this point is not constrained */
4161: newNumPoints++;
4162: PetscSectionGetDof(section,b,&bDof);
4163: newNumIndices += bDof;
4164: for (f = 0; f < numFields; ++f) {
4165: PetscInt fDof;
4167: PetscSectionGetFieldDof(section, b, f, &fDof);
4168: newOffsets[f+1] += fDof;
4169: }
4170: }
4171: }
4172: }
4173: if (!anyConstrained) {
4174: *outNumPoints = 0;
4175: *outNumIndices = 0;
4176: *outPoints = NULL;
4177: *outValues = NULL;
4178: if (aSec) {
4179: ISRestoreIndices(aIS,&anchors);
4180: }
4181: return(0);
4182: }
4184: for (f = 1; f < numFields; ++f) newOffsets[f+1] += newOffsets[f];
4186: if (numFields && newOffsets[numFields] != newNumIndices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", newOffsets[numFields], newNumIndices);
4188: DMGetDefaultConstraints(dm, &cSec, &cMat);
4190: /* output arrays */
4191: DMGetWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4192: DMGetWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4194: /* workspaces */
4195: DMGetWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4196: if (numFields) {
4197: for (f = 0; f < numFields; f++) {
4198: DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4199: DMGetWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4200: }
4201: }
4202: else {
4203: DMGetWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4204: DMGetWorkArray(dm,numPoints,PETSC_INT,&newPointOffsets[0]);
4205: }
4207: /* get workspaces for the point-to-point matrices */
4208: if (numFields) {
4209: for (p = 0; p < numPoints; p++) {
4210: PetscInt b = points[2*p];
4211: PetscInt bDof = 0;
4213: if (b >= aStart && b < aEnd) {
4214: PetscSectionGetDof(aSec, b, &bDof);
4215: }
4216: if (bDof) {
4217: for (f = 0; f < numFields; f++) {
4218: PetscInt fDof, q, bOff, allFDof = 0;
4220: PetscSectionGetFieldDof(section, b, f, &fDof);
4221: PetscSectionGetOffset(aSec, b, &bOff);
4222: for (q = 0; q < bDof; q++) {
4223: PetscInt a = anchors[bOff + q];
4224: PetscInt aFDof;
4226: PetscSectionGetFieldDof(section, a, f, &aFDof);
4227: allFDof += aFDof;
4228: }
4229: newPointOffsets[f][p+1] = allFDof;
4230: pointMatOffsets[f][p+1] = fDof * allFDof;
4231: }
4232: }
4233: else {
4234: for (f = 0; f < numFields; f++) {
4235: PetscInt fDof;
4237: PetscSectionGetFieldDof(section, b, f, &fDof);
4238: newPointOffsets[f][p+1] = fDof;
4239: pointMatOffsets[f][p+1] = 0;
4240: }
4241: }
4242: }
4243: for (f = 0; f < numFields; f++) {
4244: newPointOffsets[f][0] = 0;
4245: pointMatOffsets[f][0] = 0;
4246: for (p = 0; p < numPoints; p++) {
4247: newPointOffsets[f][p+1] += newPointOffsets[f][p];
4248: pointMatOffsets[f][p+1] += pointMatOffsets[f][p];
4249: }
4250: DMGetWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4251: }
4252: }
4253: else {
4254: for (p = 0; p < numPoints; p++) {
4255: PetscInt b = points[2*p];
4256: PetscInt bDof = 0;
4258: if (b >= aStart && b < aEnd) {
4259: PetscSectionGetDof(aSec, b, &bDof);
4260: }
4261: if (bDof) {
4262: PetscInt dof, bOff, q, allDof = 0;
4264: PetscSectionGetDof(section, b, &dof);
4265: PetscSectionGetOffset(aSec, b, &bOff);
4266: for (q = 0; q < bDof; q++) {
4267: PetscInt a = anchors[bOff + q], aDof;
4269: PetscSectionGetDof(section, a, &aDof);
4270: allDof += aDof;
4271: }
4272: newPointOffsets[0][p+1] = allDof;
4273: pointMatOffsets[0][p+1] = dof * allDof;
4274: }
4275: else {
4276: PetscInt dof;
4278: PetscSectionGetDof(section, b, &dof);
4279: newPointOffsets[0][p+1] = dof;
4280: pointMatOffsets[0][p+1] = 0;
4281: }
4282: }
4283: newPointOffsets[0][0] = 0;
4284: pointMatOffsets[0][0] = 0;
4285: for (p = 0; p < numPoints; p++) {
4286: newPointOffsets[0][p+1] += newPointOffsets[0][p];
4287: pointMatOffsets[0][p+1] += pointMatOffsets[0][p];
4288: }
4289: DMGetWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4290: }
4292: /* get the point-to-point matrices; construct newPoints */
4293: PetscSectionGetMaxDof(aSec, &maxAnchor);
4294: PetscSectionGetMaxDof(section, &maxDof);
4295: DMGetWorkArray(dm,maxDof,PETSC_INT,&indices);
4296: DMGetWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4297: if (numFields) {
4298: for (p = 0, newP = 0; p < numPoints; p++) {
4299: PetscInt b = points[2*p];
4300: PetscInt o = points[2*p+1];
4301: PetscInt bDof = 0;
4303: if (b >= aStart && b < aEnd) {
4304: PetscSectionGetDof(aSec, b, &bDof);
4305: }
4306: if (bDof) {
4307: PetscInt fStart[32], fEnd[32], fAnchorStart[32], fAnchorEnd[32], bOff, q;
4309: fStart[0] = 0;
4310: fEnd[0] = 0;
4311: for (f = 0; f < numFields; f++) {
4312: PetscInt fDof;
4314: PetscSectionGetFieldDof(cSec, b, f, &fDof);
4315: fStart[f+1] = fStart[f] + fDof;
4316: fEnd[f+1] = fStart[f+1];
4317: }
4318: PetscSectionGetOffset(cSec, b, &bOff);
4319: indicesPointFields_private(cSec, b, bOff, fEnd, PETSC_TRUE, o, indices);
4321: fAnchorStart[0] = 0;
4322: fAnchorEnd[0] = 0;
4323: for (f = 0; f < numFields; f++) {
4324: PetscInt fDof = newPointOffsets[f][p + 1] - newPointOffsets[f][p];
4326: fAnchorStart[f+1] = fAnchorStart[f] + fDof;
4327: fAnchorEnd[f+1] = fAnchorStart[f + 1];
4328: }
4329: PetscSectionGetOffset (aSec, b, &bOff);
4330: for (q = 0; q < bDof; q++) {
4331: PetscInt a = anchors[bOff + q], aOff;
4333: /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4334: newPoints[2*(newP + q)] = a;
4335: newPoints[2*(newP + q) + 1] = 0;
4336: PetscSectionGetOffset(section, a, &aOff);
4337: indicesPointFields_private(section, a, aOff, fAnchorEnd, PETSC_TRUE, 0, newIndices);
4338: }
4339: newP += bDof;
4341: /* get the point-to-point submatrix */
4342: for (f = 0; f < numFields; f++) {
4343: MatGetValues(cMat,fEnd[f]-fStart[f],indices + fStart[f],fAnchorEnd[f] - fAnchorStart[f],newIndices + fAnchorStart[f],pointMat[f] + pointMatOffsets[f][p]);
4344: }
4345: }
4346: else {
4347: newPoints[2 * newP] = b;
4348: newPoints[2 * newP + 1] = o;
4349: newP++;
4350: }
4351: }
4352: } else {
4353: for (p = 0; p < numPoints; p++) {
4354: PetscInt b = points[2*p];
4355: PetscInt o = points[2*p+1];
4356: PetscInt bDof = 0;
4358: if (b >= aStart && b < aEnd) {
4359: PetscSectionGetDof(aSec, b, &bDof);
4360: }
4361: if (bDof) {
4362: PetscInt bEnd = 0, bAnchorEnd = 0, bOff;
4364: PetscSectionGetOffset(cSec, b, &bOff);
4365: indicesPoint_private(cSec, b, bOff, &bEnd, PETSC_TRUE, o, indices);
4367: PetscSectionGetOffset (aSec, b, &bOff);
4368: for (q = 0; q < bDof; q++) {
4369: PetscInt a = anchors[bOff + q], aOff;
4371: /* we take the orientation of ap into account in the order that we constructed the indices above: the newly added points have no orientation */
4373: newPoints[2*(newP + q)] = a;
4374: newPoints[2*(newP + q) + 1] = 0;
4375: PetscSectionGetOffset(section, a, &aOff);
4376: indicesPoint_private(section, a, aOff, &bAnchorEnd, PETSC_TRUE, 0, newIndices);
4377: }
4378: newP += bDof;
4380: /* get the point-to-point submatrix */
4381: MatGetValues(cMat,bEnd,indices,bAnchorEnd,newIndices,pointMat[0] + pointMatOffsets[0][p]);
4382: }
4383: else {
4384: newPoints[2 * newP] = b;
4385: newPoints[2 * newP + 1] = o;
4386: newP++;
4387: }
4388: }
4389: }
4391: PetscMemzero(tmpValues,newNumIndices*numIndices*sizeof(*tmpValues));
4392: /* multiply constraints on the right */
4393: if (numFields) {
4394: for (f = 0; f < numFields; f++) {
4395: PetscInt oldOff = offsets[f];
4397: for (p = 0; p < numPoints; p++) {
4398: PetscInt cStart = newPointOffsets[f][p];
4399: PetscInt b = points[2 * p];
4400: PetscInt c, r, k;
4401: PetscInt dof;
4403: PetscSectionGetFieldDof(section,b,f,&dof);
4404: if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4405: PetscInt nCols = newPointOffsets[f][p+1]-cStart;
4406: const PetscScalar *mat = pointMat[f] + pointMatOffsets[f][p];
4408: for (r = 0; r < numIndices; r++) {
4409: for (c = 0; c < nCols; c++) {
4410: for (k = 0; k < dof; k++) {
4411: tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4412: }
4413: }
4414: }
4415: }
4416: else {
4417: /* copy this column as is */
4418: for (r = 0; r < numIndices; r++) {
4419: for (c = 0; c < dof; c++) {
4420: tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4421: }
4422: }
4423: }
4424: oldOff += dof;
4425: }
4426: }
4427: }
4428: else {
4429: PetscInt oldOff = 0;
4430: for (p = 0; p < numPoints; p++) {
4431: PetscInt cStart = newPointOffsets[0][p];
4432: PetscInt b = points[2 * p];
4433: PetscInt c, r, k;
4434: PetscInt dof;
4436: PetscSectionGetDof(section,b,&dof);
4437: if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4438: PetscInt nCols = newPointOffsets[0][p+1]-cStart;
4439: const PetscScalar *mat = pointMat[0] + pointMatOffsets[0][p];
4441: for (r = 0; r < numIndices; r++) {
4442: for (c = 0; c < nCols; c++) {
4443: for (k = 0; k < dof; k++) {
4444: tmpValues[r * newNumIndices + cStart + c] += mat[k * nCols + c] * values[r * numIndices + oldOff + k];
4445: }
4446: }
4447: }
4448: }
4449: else {
4450: /* copy this column as is */
4451: for (r = 0; r < numIndices; r++) {
4452: for (c = 0; c < dof; c++) {
4453: tmpValues[r * newNumIndices + cStart + c] = values[r * numIndices + oldOff + c];
4454: }
4455: }
4456: }
4457: oldOff += dof;
4458: }
4459: }
4461: PetscMemzero(newValues,newNumIndices*newNumIndices*sizeof(*newValues));
4462: /* multiply constraints transpose on the left */
4463: if (numFields) {
4464: for (f = 0; f < numFields; f++) {
4465: PetscInt oldOff = offsets[f];
4467: for (p = 0; p < numPoints; p++) {
4468: PetscInt rStart = newPointOffsets[f][p];
4469: PetscInt b = points[2 * p];
4470: PetscInt c, r, k;
4471: PetscInt dof;
4473: PetscSectionGetFieldDof(section,b,f,&dof);
4474: if (pointMatOffsets[f][p] < pointMatOffsets[f][p + 1]) {
4475: PetscInt nRows = newPointOffsets[f][p+1]-rStart;
4476: const PetscScalar *PETSC_RESTRICT mat = pointMat[f] + pointMatOffsets[f][p];
4478: for (r = 0; r < nRows; r++) {
4479: for (c = 0; c < newNumIndices; c++) {
4480: for (k = 0; k < dof; k++) {
4481: newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4482: }
4483: }
4484: }
4485: }
4486: else {
4487: /* copy this row as is */
4488: for (r = 0; r < dof; r++) {
4489: for (c = 0; c < newNumIndices; c++) {
4490: newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4491: }
4492: }
4493: }
4494: oldOff += dof;
4495: }
4496: }
4497: }
4498: else {
4499: PetscInt oldOff = 0;
4501: for (p = 0; p < numPoints; p++) {
4502: PetscInt rStart = newPointOffsets[0][p];
4503: PetscInt b = points[2 * p];
4504: PetscInt c, r, k;
4505: PetscInt dof;
4507: PetscSectionGetDof(section,b,&dof);
4508: if (pointMatOffsets[0][p] < pointMatOffsets[0][p + 1]) {
4509: PetscInt nRows = newPointOffsets[0][p+1]-rStart;
4510: const PetscScalar *PETSC_RESTRICT mat = pointMat[0] + pointMatOffsets[0][p];
4512: for (r = 0; r < nRows; r++) {
4513: for (c = 0; c < newNumIndices; c++) {
4514: for (k = 0; k < dof; k++) {
4515: newValues[(rStart + r) * newNumIndices + c] += mat[k * nRows + r] * tmpValues[(oldOff + k) * newNumIndices + c];
4516: }
4517: }
4518: }
4519: }
4520: else {
4521: /* copy this row as is */
4522: for (r = 0; r < dof; r++) {
4523: for (c = 0; c < newNumIndices; c++) {
4524: newValues[(rStart + r) * newNumIndices + c] = tmpValues[(oldOff + r) * newNumIndices + c];
4525: }
4526: }
4527: }
4528: oldOff += dof;
4529: }
4530: }
4532: /* clean up */
4533: DMRestoreWorkArray(dm,maxDof,PETSC_INT,&indices);
4534: DMRestoreWorkArray(dm,maxAnchor*maxDof,PETSC_INT,&newIndices);
4535: if (numFields) {
4536: for (f = 0; f < numFields; f++) {
4537: DMRestoreWorkArray(dm,pointMatOffsets[f][numPoints],PETSC_SCALAR,&pointMat[f]);
4538: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[f]);
4539: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[f]);
4540: }
4541: }
4542: else {
4543: DMRestoreWorkArray(dm,pointMatOffsets[0][numPoints],PETSC_SCALAR,&pointMat[0]);
4544: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&pointMatOffsets[0]);
4545: DMRestoreWorkArray(dm,numPoints+1,PETSC_INT,&newPointOffsets[0]);
4546: }
4547: DMRestoreWorkArray(dm,newNumIndices*numIndices,PETSC_SCALAR,&tmpValues);
4548: ISRestoreIndices(aIS,&anchors);
4550: /* output */
4551: *outNumPoints = newNumPoints;
4552: *outNumIndices = newNumIndices;
4553: *outPoints = newPoints;
4554: *outValues = newValues;
4555: for (f = 0; f < numFields; f++) {
4556: offsets[f] = newOffsets[f];
4557: }
4558: return(0);
4559: }
4563: /*@C
4564: DMPlexMatSetClosure - Set an array of the values on the closure of 'point'
4566: Not collective
4568: Input Parameters:
4569: + dm - The DM
4570: . section - The section describing the layout in v, or NULL to use the default section
4571: . globalSection - The section describing the layout in v, or NULL to use the default global section
4572: . A - The matrix
4573: . point - The sieve point in the DM
4574: . values - The array of values
4575: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
4577: Fortran Notes:
4578: This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
4580: Level: intermediate
4582: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
4583: @*/
4584: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4585: {
4586: DM_Plex *mesh = (DM_Plex*) dm->data;
4587: PetscSection clSection;
4588: IS clPoints;
4589: PetscInt *points = NULL, *newPoints;
4590: const PetscInt *clp;
4591: PetscInt *indices;
4592: PetscInt offsets[32];
4593: PetscInt numFields, numPoints, newNumPoints, numIndices, newNumIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
4594: PetscScalar *newValues;
4595: PetscErrorCode ierr;
4599: if (!section) {DMGetDefaultSection(dm, §ion);}
4601: if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
4604: PetscSectionGetNumFields(section, &numFields);
4605: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4606: PetscMemzero(offsets, 32 * sizeof(PetscInt));
4607: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4608: if (!clPoints) {
4609: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4610: /* Compress out points not in the section */
4611: PetscSectionGetChart(section, &pStart, &pEnd);
4612: for (p = 0, q = 0; p < numPoints*2; p += 2) {
4613: if ((points[p] >= pStart) && (points[p] < pEnd)) {
4614: points[q*2] = points[p];
4615: points[q*2+1] = points[p+1];
4616: ++q;
4617: }
4618: }
4619: numPoints = q;
4620: } else {
4621: PetscInt dof, off;
4623: PetscSectionGetDof(clSection, point, &dof);
4624: numPoints = dof/2;
4625: PetscSectionGetOffset(clSection, point, &off);
4626: ISGetIndices(clPoints, &clp);
4627: points = (PetscInt *) &clp[off];
4628: }
4629: for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
4630: PetscInt fdof;
4632: PetscSectionGetDof(section, points[p], &dof);
4633: for (f = 0; f < numFields; ++f) {
4634: PetscSectionGetFieldDof(section, points[p], f, &fdof);
4635: offsets[f+1] += fdof;
4636: }
4637: numIndices += dof;
4638: }
4639: for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];
4641: if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
4642: DMPlexAnchorsModifyMat(dm,section,numPoints,numIndices,points,values,&newNumPoints,&newNumIndices,&newPoints,&newValues,offsets);
4643: if (newNumPoints) {
4644: if (!clPoints) {
4645: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4646: } else {
4647: ISRestoreIndices(clPoints, &clp);
4648: }
4649: numPoints = newNumPoints;
4650: numIndices = newNumIndices;
4651: points = newPoints;
4652: values = newValues;
4653: }
4654: DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
4655: if (numFields) {
4656: for (p = 0; p < numPoints*2; p += 2) {
4657: PetscInt o = points[p+1];
4658: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4659: indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
4660: }
4661: } else {
4662: for (p = 0, off = 0; p < numPoints*2; p += 2) {
4663: PetscInt o = points[p+1];
4664: PetscSectionGetOffset(globalSection, points[p], &globalOff);
4665: indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
4666: }
4667: }
4668: if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
4669: MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
4670: if (mesh->printFEM > 1) {
4671: PetscInt i;
4672: PetscPrintf(PETSC_COMM_SELF, " Indices:");
4673: for (i = 0; i < numIndices; ++i) {PetscPrintf(PETSC_COMM_SELF, " %d", indices[i]);}
4674: PetscPrintf(PETSC_COMM_SELF, "\n");
4675: }
4676: if (ierr) {
4677: PetscMPIInt rank;
4678: PetscErrorCode ierr2;
4680: ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4681: ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4682: ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
4683: ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
4684:
4685: }
4686: if (newNumPoints) {
4687: DMRestoreWorkArray(dm,newNumIndices*newNumIndices,PETSC_SCALAR,&newValues);
4688: DMRestoreWorkArray(dm,2*newNumPoints,PETSC_INT,&newPoints);
4689: }
4690: else {
4691: if (!clPoints) {
4692: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4693: } else {
4694: ISRestoreIndices(clPoints, &clp);
4695: }
4696: }
4697: DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
4698: return(0);
4699: }
4703: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
4704: {
4705: DM_Plex *mesh = (DM_Plex*) dmf->data;
4706: PetscInt *fpoints = NULL, *ftotpoints = NULL;
4707: PetscInt *cpoints = NULL;
4708: PetscInt *findices, *cindices;
4709: PetscInt foffsets[32], coffsets[32];
4710: CellRefiner cellRefiner;
4711: PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4712: PetscErrorCode ierr;
4717: if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4719: if (!csection) {DMGetDefaultSection(dmc, &csection);}
4721: if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4723: if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4726: PetscSectionGetNumFields(fsection, &numFields);
4727: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4728: PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4729: PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4730: /* Column indices */
4731: DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4732: maxFPoints = numCPoints;
4733: /* Compress out points not in the section */
4734: /* TODO: Squeeze out points with 0 dof as well */
4735: PetscSectionGetChart(csection, &pStart, &pEnd);
4736: for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4737: if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4738: cpoints[q*2] = cpoints[p];
4739: cpoints[q*2+1] = cpoints[p+1];
4740: ++q;
4741: }
4742: }
4743: numCPoints = q;
4744: for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4745: PetscInt fdof;
4747: PetscSectionGetDof(csection, cpoints[p], &dof);
4748: if (!dof) continue;
4749: for (f = 0; f < numFields; ++f) {
4750: PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4751: coffsets[f+1] += fdof;
4752: }
4753: numCIndices += dof;
4754: }
4755: for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4756: /* Row indices */
4757: DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4758: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4759: DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4760: for (r = 0, q = 0; r < numSubcells; ++r) {
4761: /* TODO Map from coarse to fine cells */
4762: DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4763: /* Compress out points not in the section */
4764: PetscSectionGetChart(fsection, &pStart, &pEnd);
4765: for (p = 0; p < numFPoints*2; p += 2) {
4766: if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4767: PetscSectionGetDof(fsection, fpoints[p], &dof);
4768: if (!dof) continue;
4769: for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4770: if (s < q) continue;
4771: ftotpoints[q*2] = fpoints[p];
4772: ftotpoints[q*2+1] = fpoints[p+1];
4773: ++q;
4774: }
4775: }
4776: DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4777: }
4778: numFPoints = q;
4779: for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4780: PetscInt fdof;
4782: PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4783: if (!dof) continue;
4784: for (f = 0; f < numFields; ++f) {
4785: PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4786: foffsets[f+1] += fdof;
4787: }
4788: numFIndices += dof;
4789: }
4790: for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];
4792: if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4793: if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4794: DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4795: DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4796: if (numFields) {
4797: for (p = 0; p < numFPoints*2; p += 2) {
4798: PetscInt o = ftotpoints[p+1];
4799: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4800: indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4801: }
4802: for (p = 0; p < numCPoints*2; p += 2) {
4803: PetscInt o = cpoints[p+1];
4804: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4805: indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4806: }
4807: } else {
4808: for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4809: PetscInt o = ftotpoints[p+1];
4810: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4811: indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4812: }
4813: for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4814: PetscInt o = cpoints[p+1];
4815: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4816: indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4817: }
4818: }
4819: if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
4820: MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
4821: if (ierr) {
4822: PetscMPIInt rank;
4823: PetscErrorCode ierr2;
4825: ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
4826: ierr2 = (*PetscErrorPrintf)("[%d]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
4827: ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
4828: ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
4829: ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
4830:
4831: }
4832: DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4833: DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4834: DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
4835: DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
4836: return(0);
4837: }
4841: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
4842: {
4843: PetscInt *fpoints = NULL, *ftotpoints = NULL;
4844: PetscInt *cpoints = NULL;
4845: PetscInt foffsets[32], coffsets[32];
4846: CellRefiner cellRefiner;
4847: PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
4853: if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
4855: if (!csection) {DMGetDefaultSection(dmc, &csection);}
4857: if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
4859: if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
4861: PetscSectionGetNumFields(fsection, &numFields);
4862: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
4863: PetscMemzero(foffsets, 32 * sizeof(PetscInt));
4864: PetscMemzero(coffsets, 32 * sizeof(PetscInt));
4865: /* Column indices */
4866: DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4867: maxFPoints = numCPoints;
4868: /* Compress out points not in the section */
4869: /* TODO: Squeeze out points with 0 dof as well */
4870: PetscSectionGetChart(csection, &pStart, &pEnd);
4871: for (p = 0, q = 0; p < numCPoints*2; p += 2) {
4872: if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
4873: cpoints[q*2] = cpoints[p];
4874: cpoints[q*2+1] = cpoints[p+1];
4875: ++q;
4876: }
4877: }
4878: numCPoints = q;
4879: for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
4880: PetscInt fdof;
4882: PetscSectionGetDof(csection, cpoints[p], &dof);
4883: if (!dof) continue;
4884: for (f = 0; f < numFields; ++f) {
4885: PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
4886: coffsets[f+1] += fdof;
4887: }
4888: numCIndices += dof;
4889: }
4890: for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
4891: /* Row indices */
4892: DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
4893: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
4894: DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
4895: for (r = 0, q = 0; r < numSubcells; ++r) {
4896: /* TODO Map from coarse to fine cells */
4897: DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
4898: /* Compress out points not in the section */
4899: PetscSectionGetChart(fsection, &pStart, &pEnd);
4900: for (p = 0; p < numFPoints*2; p += 2) {
4901: if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
4902: PetscSectionGetDof(fsection, fpoints[p], &dof);
4903: if (!dof) continue;
4904: for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
4905: if (s < q) continue;
4906: ftotpoints[q*2] = fpoints[p];
4907: ftotpoints[q*2+1] = fpoints[p+1];
4908: ++q;
4909: }
4910: }
4911: DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
4912: }
4913: numFPoints = q;
4914: for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
4915: PetscInt fdof;
4917: PetscSectionGetDof(fsection, ftotpoints[p], &dof);
4918: if (!dof) continue;
4919: for (f = 0; f < numFields; ++f) {
4920: PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
4921: foffsets[f+1] += fdof;
4922: }
4923: numFIndices += dof;
4924: }
4925: for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];
4927: if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
4928: if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
4929: if (numFields) {
4930: for (p = 0; p < numFPoints*2; p += 2) {
4931: PetscInt o = ftotpoints[p+1];
4932: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4933: indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
4934: }
4935: for (p = 0; p < numCPoints*2; p += 2) {
4936: PetscInt o = cpoints[p+1];
4937: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4938: indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
4939: }
4940: } else {
4941: for (p = 0, off = 0; p < numFPoints*2; p += 2) {
4942: PetscInt o = ftotpoints[p+1];
4943: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
4944: indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
4945: }
4946: for (p = 0, off = 0; p < numCPoints*2; p += 2) {
4947: PetscInt o = cpoints[p+1];
4948: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
4949: indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
4950: }
4951: }
4952: DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
4953: DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
4954: return(0);
4955: }
4959: /*@
4960: DMPlexGetHybridBounds - Get the first mesh point of each dimension which is a hybrid
4962: Input Parameter:
4963: . dm - The DMPlex object
4965: Output Parameters:
4966: + cMax - The first hybrid cell
4967: . fMax - The first hybrid face
4968: . eMax - The first hybrid edge
4969: - vMax - The first hybrid vertex
4971: Level: developer
4973: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
4974: @*/
4975: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
4976: {
4977: DM_Plex *mesh = (DM_Plex*) dm->data;
4978: PetscInt dim;
4983: DMGetDimension(dm, &dim);
4984: if (cMax) *cMax = mesh->hybridPointMax[dim];
4985: if (fMax) *fMax = mesh->hybridPointMax[dim-1];
4986: if (eMax) *eMax = mesh->hybridPointMax[1];
4987: if (vMax) *vMax = mesh->hybridPointMax[0];
4988: return(0);
4989: }
4993: /*@
4994: DMPlexSetHybridBounds - Set the first mesh point of each dimension which is a hybrid
4996: Input Parameters:
4997: . dm - The DMPlex object
4998: . cMax - The first hybrid cell
4999: . fMax - The first hybrid face
5000: . eMax - The first hybrid edge
5001: - vMax - The first hybrid vertex
5003: Level: developer
5005: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5006: @*/
5007: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5008: {
5009: DM_Plex *mesh = (DM_Plex*) dm->data;
5010: PetscInt dim;
5015: DMGetDimension(dm, &dim);
5016: if (cMax >= 0) mesh->hybridPointMax[dim] = cMax;
5017: if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5018: if (eMax >= 0) mesh->hybridPointMax[1] = eMax;
5019: if (vMax >= 0) mesh->hybridPointMax[0] = vMax;
5020: return(0);
5021: }
5025: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5026: {
5027: DM_Plex *mesh = (DM_Plex*) dm->data;
5032: *cellHeight = mesh->vtkCellHeight;
5033: return(0);
5034: }
5038: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5039: {
5040: DM_Plex *mesh = (DM_Plex*) dm->data;
5044: mesh->vtkCellHeight = cellHeight;
5045: return(0);
5046: }
5050: /* We can easily have a form that takes an IS instead */
5051: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5052: {
5053: PetscSection section, globalSection;
5054: PetscInt *numbers, p;
5058: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
5059: PetscSectionSetChart(section, pStart, pEnd);
5060: for (p = pStart; p < pEnd; ++p) {
5061: PetscSectionSetDof(section, p, 1);
5062: }
5063: PetscSectionSetUp(section);
5064: PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, &globalSection);
5065: PetscMalloc1(pEnd - pStart, &numbers);
5066: for (p = pStart; p < pEnd; ++p) {
5067: PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5068: if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5069: else numbers[p-pStart] += shift;
5070: }
5071: ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5072: if (globalSize) {
5073: PetscLayout layout;
5074: PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5075: PetscLayoutGetSize(layout, globalSize);
5076: PetscLayoutDestroy(&layout);
5077: }
5078: PetscSectionDestroy(§ion);
5079: PetscSectionDestroy(&globalSection);
5080: return(0);
5081: }
5085: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5086: {
5087: DM_Plex *mesh = (DM_Plex*) dm->data;
5088: PetscInt cellHeight, cStart, cEnd, cMax;
5093: if (!mesh->globalCellNumbers) {
5094: DMPlexGetVTKCellHeight(dm, &cellHeight);
5095: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5096: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5097: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
5098: DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, &mesh->globalCellNumbers);
5099: }
5100: *globalCellNumbers = mesh->globalCellNumbers;
5101: return(0);
5102: }
5106: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5107: {
5108: DM_Plex *mesh = (DM_Plex*) dm->data;
5109: PetscInt vStart, vEnd, vMax;
5114: if (!mesh->globalVertexNumbers) {
5115: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5116: DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5117: if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
5118: DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, &mesh->globalVertexNumbers);
5119: }
5120: *globalVertexNumbers = mesh->globalVertexNumbers;
5121: return(0);
5122: }
5126: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5127: {
5128: IS nums[4];
5129: PetscInt depths[4];
5130: PetscInt depth, d, shift = 0;
5135: DMPlexGetDepth(dm, &depth);
5136: /* For unstratified meshes use dim instead of depth */
5137: if (depth < 0) {DMGetDimension(dm, &depth);}
5138: depths[0] = depth; depths[1] = 0;
5139: for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5140: for (d = 0; d <= depth; ++d) {
5141: PetscInt pStart, pEnd, gsize;
5143: DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5144: DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5145: shift += gsize;
5146: }
5147: ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5148: for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5149: return(0);
5150: }
5155: /*@C
5156: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
5157: the local section and an SF describing the section point overlap.
5159: Input Parameters:
5160: + s - The PetscSection for the local field layout
5161: . sf - The SF describing parallel layout of the section points
5162: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
5163: . label - The label specifying the points
5164: - labelValue - The label stratum specifying the points
5166: Output Parameter:
5167: . gsection - The PetscSection for the global field layout
5169: Note: This gives negative sizes and offsets to points not owned by this process
5171: Level: developer
5173: .seealso: PetscSectionCreate()
5174: @*/
5175: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
5176: {
5177: PetscInt *neg = NULL, *tmpOff = NULL;
5178: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
5182: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
5183: PetscSectionGetChart(s, &pStart, &pEnd);
5184: PetscSectionSetChart(*gsection, pStart, pEnd);
5185: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
5186: if (nroots >= 0) {
5187: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
5188: PetscCalloc1(nroots, &neg);
5189: if (nroots > pEnd-pStart) {
5190: PetscCalloc1(nroots, &tmpOff);
5191: } else {
5192: tmpOff = &(*gsection)->atlasDof[-pStart];
5193: }
5194: }
5195: /* Mark ghost points with negative dof */
5196: for (p = pStart; p < pEnd; ++p) {
5197: PetscInt value;
5199: DMLabelGetValue(label, p, &value);
5200: if (value != labelValue) continue;
5201: PetscSectionGetDof(s, p, &dof);
5202: PetscSectionSetDof(*gsection, p, dof);
5203: PetscSectionGetConstraintDof(s, p, &cdof);
5204: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
5205: if (neg) neg[p] = -(dof+1);
5206: }
5207: PetscSectionSetUpBC(*gsection);
5208: if (nroots >= 0) {
5209: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5210: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5211: if (nroots > pEnd-pStart) {
5212: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
5213: }
5214: }
5215: /* Calculate new sizes, get proccess offset, and calculate point offsets */
5216: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5217: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
5218: (*gsection)->atlasOff[p] = off;
5219: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
5220: }
5221: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
5222: globalOff -= off;
5223: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5224: (*gsection)->atlasOff[p] += globalOff;
5225: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
5226: }
5227: /* Put in negative offsets for ghost points */
5228: if (nroots >= 0) {
5229: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5230: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5231: if (nroots > pEnd-pStart) {
5232: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
5233: }
5234: }
5235: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
5236: PetscFree(neg);
5237: return(0);
5238: }
5242: /*@
5243: DMPlexCheckSymmetry - Check that the adjacency information in the mesh is symmetric.
5245: Input Parameters:
5246: + dm - The DMPlex object
5248: Note: This is a useful diagnostic when creating meshes programmatically.
5250: Level: developer
5252: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5253: @*/
5254: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5255: {
5256: PetscSection coneSection, supportSection;
5257: const PetscInt *cone, *support;
5258: PetscInt coneSize, c, supportSize, s;
5259: PetscInt pStart, pEnd, p, csize, ssize;
5260: PetscErrorCode ierr;
5264: DMPlexGetConeSection(dm, &coneSection);
5265: DMPlexGetSupportSection(dm, &supportSection);
5266: /* Check that point p is found in the support of its cone points, and vice versa */
5267: DMPlexGetChart(dm, &pStart, &pEnd);
5268: for (p = pStart; p < pEnd; ++p) {
5269: DMPlexGetConeSize(dm, p, &coneSize);
5270: DMPlexGetCone(dm, p, &cone);
5271: for (c = 0; c < coneSize; ++c) {
5272: PetscBool dup = PETSC_FALSE;
5273: PetscInt d;
5274: for (d = c-1; d >= 0; --d) {
5275: if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5276: }
5277: DMPlexGetSupportSize(dm, cone[c], &supportSize);
5278: DMPlexGetSupport(dm, cone[c], &support);
5279: for (s = 0; s < supportSize; ++s) {
5280: if (support[s] == p) break;
5281: }
5282: if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5283: PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5284: for (s = 0; s < coneSize; ++s) {
5285: PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
5286: }
5287: PetscPrintf(PETSC_COMM_SELF, "\n");
5288: PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
5289: for (s = 0; s < supportSize; ++s) {
5290: PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
5291: }
5292: PetscPrintf(PETSC_COMM_SELF, "\n");
5293: if (dup) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
5294: else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
5295: }
5296: }
5297: DMPlexGetSupportSize(dm, p, &supportSize);
5298: DMPlexGetSupport(dm, p, &support);
5299: for (s = 0; s < supportSize; ++s) {
5300: DMPlexGetConeSize(dm, support[s], &coneSize);
5301: DMPlexGetCone(dm, support[s], &cone);
5302: for (c = 0; c < coneSize; ++c) {
5303: if (cone[c] == p) break;
5304: }
5305: if (c >= coneSize) {
5306: PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
5307: for (c = 0; c < supportSize; ++c) {
5308: PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
5309: }
5310: PetscPrintf(PETSC_COMM_SELF, "\n");
5311: PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
5312: for (c = 0; c < coneSize; ++c) {
5313: PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
5314: }
5315: PetscPrintf(PETSC_COMM_SELF, "\n");
5316: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
5317: }
5318: }
5319: }
5320: PetscSectionGetStorageSize(coneSection, &csize);
5321: PetscSectionGetStorageSize(supportSection, &ssize);
5322: if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
5323: return(0);
5324: }
5328: /*@
5329: DMPlexCheckSkeleton - Check that each cell has the correct number of vertices
5331: Input Parameters:
5332: + dm - The DMPlex object
5333: . isSimplex - Are the cells simplices or tensor products
5334: - cellHeight - Normally 0
5336: Note: This is a useful diagnostic when creating meshes programmatically.
5338: Level: developer
5340: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
5341: @*/
5342: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5343: {
5344: PetscInt dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;
5349: DMGetDimension(dm, &dim);
5350: switch (dim) {
5351: case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
5352: case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
5353: case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
5354: default:
5355: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
5356: }
5357: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5358: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5359: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5360: cMax = cMax >= 0 ? cMax : cEnd;
5361: for (c = cStart; c < cMax; ++c) {
5362: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
5364: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5365: for (cl = 0; cl < closureSize*2; cl += 2) {
5366: const PetscInt p = closure[cl];
5367: if ((p >= vStart) && (p < vEnd)) ++coneSize;
5368: }
5369: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5370: if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d vertices != %d", c, coneSize, numCorners);
5371: }
5372: for (c = cMax; c < cEnd; ++c) {
5373: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
5375: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5376: for (cl = 0; cl < closureSize*2; cl += 2) {
5377: const PetscInt p = closure[cl];
5378: if ((p >= vStart) && (p < vEnd)) ++coneSize;
5379: }
5380: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5381: if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has %d vertices > %d", c, coneSize, numHybridCorners);
5382: }
5383: return(0);
5384: }
5388: /*@
5389: DMPlexCheckFaces - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type
5391: Input Parameters:
5392: + dm - The DMPlex object
5393: . isSimplex - Are the cells simplices or tensor products
5394: - cellHeight - Normally 0
5396: Note: This is a useful diagnostic when creating meshes programmatically.
5398: Level: developer
5400: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
5401: @*/
5402: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
5403: {
5404: PetscInt pMax[4];
5405: PetscInt dim, vStart, vEnd, cStart, cEnd, c, h;
5410: DMGetDimension(dm, &dim);
5411: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5412: DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
5413: for (h = cellHeight; h < dim; ++h) {
5414: DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
5415: for (c = cStart; c < cEnd; ++c) {
5416: const PetscInt *cone, *ornt, *faces;
5417: PetscInt numFaces, faceSize, coneSize,f;
5418: PetscInt *closure = NULL, closureSize, cl, numCorners = 0;
5420: if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
5421: DMPlexGetConeSize(dm, c, &coneSize);
5422: DMPlexGetCone(dm, c, &cone);
5423: DMPlexGetConeOrientation(dm, c, &ornt);
5424: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5425: for (cl = 0; cl < closureSize*2; cl += 2) {
5426: const PetscInt p = closure[cl];
5427: if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
5428: }
5429: DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
5430: if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
5431: for (f = 0; f < numFaces; ++f) {
5432: PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;
5434: DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
5435: for (cl = 0; cl < fclosureSize*2; cl += 2) {
5436: const PetscInt p = fclosure[cl];
5437: if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
5438: }
5439: if (fnumCorners != faceSize) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d (%d) of cell %d has %d vertices but should have %d", cone[f], f, c, fnumCorners, faceSize);
5440: for (v = 0; v < fnumCorners; ++v) {
5441: if (fclosure[v] != faces[f*faceSize+v]) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d (%d) of cell %d vertex %d, %d != %d", cone[f], f, c, v, fclosure[v], faces[f*faceSize+v]);
5442: }
5443: DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
5444: }
5445: DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
5446: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
5447: }
5448: }
5449: return(0);
5450: }
5454: /* Pointwise interpolation
5455: Just code FEM for now
5456: u^f = I u^c
5457: sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
5458: u^f_i = sum_j psi^f_i I phi^c_j u^c_j
5459: I_{ij} = psi^f_i phi^c_j
5460: */
5461: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
5462: {
5463: PetscSection gsc, gsf;
5464: PetscInt m, n;
5465: void *ctx;
5469: /*
5470: Loop over coarse cells
5471: Loop over coarse basis functions
5472: Loop over fine cells in coarse cell
5473: Loop over fine dual basis functions
5474: Evaluate coarse basis on fine dual basis quad points
5475: Sum
5476: Update local element matrix
5477: Accumulate to interpolation matrix
5479: Can extend PetscFEIntegrateJacobian_Basic() to do a specialized cell loop
5480: */
5481: DMGetDefaultGlobalSection(dmFine, &gsf);
5482: PetscSectionGetConstrainedStorageSize(gsf, &m);
5483: DMGetDefaultGlobalSection(dmCoarse, &gsc);
5484: PetscSectionGetConstrainedStorageSize(gsc, &n);
5485: /* We need to preallocate properly */
5486: MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
5487: MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
5488: MatSetType(*interpolation, dmCoarse->mattype);
5489: DMGetApplicationContext(dmFine, &ctx);
5490: DMPlexComputeInterpolatorFEM(dmCoarse, dmFine, *interpolation, ctx);
5491: /* Use naive scaling */
5492: DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
5493: return(0);
5494: }
5498: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat)
5499: {
5501: VecScatter ctx;
5504: DMPlexComputeInjectorFEM(dmCoarse, dmFine, &ctx, NULL);
5505: MatCreateScatter(PetscObjectComm((PetscObject)ctx), ctx, mat);
5506: VecScatterDestroy(&ctx);
5507: return(0);
5508: }
5512: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
5513: {
5514: PetscSection section;
5515: IS *bcPoints, *bcComps;
5516: PetscBool *isFE;
5517: PetscInt *bcFields, *numComp, *numDof;
5518: PetscInt depth, dim, numBd, numBC = 0, numFields, bd, bc = 0, f;
5519: PetscInt cStart, cEnd, cEndInterior;
5523: DMGetNumFields(dm, &numFields);
5524: /* FE and FV boundary conditions are handled slightly differently */
5525: PetscMalloc1(numFields, &isFE);
5526: for (f = 0; f < numFields; ++f) {
5527: PetscObject obj;
5528: PetscClassId id;
5530: DMGetField(dm, f, &obj);
5531: PetscObjectGetClassId(obj, &id);
5532: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
5533: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
5534: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
5535: }
5536: /* Allocate boundary point storage for FEM boundaries */
5537: DMPlexGetDepth(dm, &depth);
5538: DMGetDimension(dm, &dim);
5539: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
5540: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
5541: DMPlexGetNumBoundary(dm, &numBd);
5542: for (bd = 0; bd < numBd; ++bd) {
5543: PetscInt field;
5544: PetscBool isEssential;
5546: DMPlexGetBoundary(dm, bd, &isEssential, NULL, NULL, &field, NULL, NULL, NULL, NULL, NULL, NULL);
5547: if (isFE[field] && isEssential) ++numBC;
5548: }
5549: /* Add ghost cell boundaries for FVM */
5550: for (f = 0; f < numFields; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC;
5551: PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);
5552: /* Constrain ghost cells for FV */
5553: for (f = 0; f < numFields; ++f) {
5554: PetscInt *newidx, c;
5556: if (isFE[f] || cEndInterior < 0) continue;
5557: PetscMalloc1(cEnd-cEndInterior,&newidx);
5558: for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c;
5559: bcFields[bc] = f;
5560: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5561: }
5562: /* Handle FEM Dirichlet boundaries */
5563: for (bd = 0; bd < numBd; ++bd) {
5564: const char *bdLabel;
5565: DMLabel label;
5566: const PetscInt *comps;
5567: const PetscInt *values;
5568: PetscInt bd2, field, numComps, numValues;
5569: PetscBool isEssential, duplicate = PETSC_FALSE;
5571: DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);
5572: if (!isFE[field]) continue;
5573: DMPlexGetLabel(dm, bdLabel, &label);
5574: /* Only want to modify label once */
5575: for (bd2 = 0; bd2 < bd; ++bd2) {
5576: const char *bdname;
5577: DMPlexGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
5578: PetscStrcmp(bdname, bdLabel, &duplicate);
5579: if (duplicate) break;
5580: }
5581: if (!duplicate && (isFE[field])) {
5582: DMPlexLabelComplete(dm, label);
5583: DMPlexLabelAddCells(dm, label);
5584: }
5585: /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */
5586: if (isEssential) {
5587: PetscInt *newidx;
5588: PetscInt n, newn = 0, p, v;
5590: bcFields[bc] = field;
5591: if (numComps) {ISCreateGeneral(PetscObjectComm((PetscObject) dm), numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);}
5592: for (v = 0; v < numValues; ++v) {
5593: IS tmp;
5594: const PetscInt *idx;
5596: DMPlexGetStratumIS(dm, bdLabel, values[v], &tmp);
5597: if (!tmp) continue;
5598: ISGetLocalSize(tmp, &n);
5599: ISGetIndices(tmp, &idx);
5600: if (isFE[field]) {
5601: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
5602: } else {
5603: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn;
5604: }
5605: ISRestoreIndices(tmp, &idx);
5606: ISDestroy(&tmp);
5607: }
5608: PetscMalloc1(newn,&newidx);
5609: newn = 0;
5610: for (v = 0; v < numValues; ++v) {
5611: IS tmp;
5612: const PetscInt *idx;
5614: DMPlexGetStratumIS(dm, bdLabel, values[v], &tmp);
5615: if (!tmp) continue;
5616: ISGetLocalSize(tmp, &n);
5617: ISGetIndices(tmp, &idx);
5618: if (isFE[field]) {
5619: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
5620: } else {
5621: for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p];
5622: }
5623: ISRestoreIndices(tmp, &idx);
5624: ISDestroy(&tmp);
5625: }
5626: ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
5627: }
5628: }
5629: /* Handle discretization */
5630: PetscCalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
5631: for (f = 0; f < numFields; ++f) {
5632: PetscObject obj;
5634: DMGetField(dm, f, &obj);
5635: if (isFE[f]) {
5636: PetscFE fe = (PetscFE) obj;
5637: const PetscInt *numFieldDof;
5638: PetscInt d;
5640: PetscFEGetNumComponents(fe, &numComp[f]);
5641: PetscFEGetNumDof(fe, &numFieldDof);
5642: for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
5643: } else {
5644: PetscFV fv = (PetscFV) obj;
5646: PetscFVGetNumComponents(fv, &numComp[f]);
5647: numDof[f*(dim+1)+dim] = numComp[f];
5648: }
5649: }
5650: for (f = 0; f < numFields; ++f) {
5651: PetscInt d;
5652: for (d = 1; d < dim; ++d) {
5653: if ((numDof[f*(dim+1)+d] > 0) && (depth < dim)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
5654: }
5655: }
5656: DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, §ion);
5657: for (f = 0; f < numFields; ++f) {
5658: PetscFE fe;
5659: const char *name;
5661: DMGetField(dm, f, (PetscObject *) &fe);
5662: PetscObjectGetName((PetscObject) fe, &name);
5663: PetscSectionSetFieldName(section, f, name);
5664: }
5665: DMSetDefaultSection(dm, section);
5666: PetscSectionDestroy(§ion);
5667: for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);ISDestroy(&bcComps[bc]);}
5668: PetscFree3(bcFields,bcPoints,bcComps);
5669: PetscFree2(numComp,numDof);
5670: PetscFree(isFE);
5671: return(0);
5672: }
5676: /*@
5677: DMPlexGetCoarseDM - Get the coarse mesh from which this was obtained by refinement
5679: Input Parameter:
5680: . dm - The DMPlex object
5682: Output Parameter:
5683: . cdm - The coarse DM
5685: Level: intermediate
5687: .seealso: DMPlexSetCoarseDM()
5688: @*/
5689: PetscErrorCode DMPlexGetCoarseDM(DM dm, DM *cdm)
5690: {
5694: *cdm = ((DM_Plex *) dm->data)->coarseMesh;
5695: return(0);
5696: }
5700: /*@
5701: DMPlexSetCoarseDM - Set the coarse mesh from which this was obtained by refinement
5703: Input Parameters:
5704: + dm - The DMPlex object
5705: - cdm - The coarse DM
5707: Level: intermediate
5709: .seealso: DMPlexGetCoarseDM()
5710: @*/
5711: PetscErrorCode DMPlexSetCoarseDM(DM dm, DM cdm)
5712: {
5713: DM_Plex *mesh;
5719: mesh = (DM_Plex *) dm->data;
5720: DMDestroy(&mesh->coarseMesh);
5721: mesh->coarseMesh = cdm;
5722: PetscObjectReference((PetscObject) mesh->coarseMesh);
5723: return(0);
5724: }
5726: /* anchors */
5729: /*@
5730: DMPlexGetAnchors - Get the layout of the anchor (point-to-point) constraints. Typically, the user will not have to
5731: call DMPlexGetAnchors() directly: if there are anchors, then DMPlexGetAnchors() is called during DMGetConstraints().
5733: not collective
5735: Input Parameters:
5736: . dm - The DMPlex object
5738: Output Parameters:
5739: + anchorSection - If not NULL, set to the section describing which points anchor the constrained points.
5740: - anchorIS - If not NULL, set to the list of anchors indexed by anchorSection
5743: Level: intermediate
5745: .seealso: DMPlexSetAnchors(), DMGetConstraints(), DMSetConstraints()
5746: @*/
5747: PetscErrorCode DMPlexGetAnchors(DM dm, PetscSection *anchorSection, IS *anchorIS)
5748: {
5749: DM_Plex *plex = (DM_Plex *)dm->data;
5754: if (!plex->anchorSection && !plex->anchorIS && plex->createanchors) {(*plex->createanchors)(dm);}
5755: if (anchorSection) *anchorSection = plex->anchorSection;
5756: if (anchorIS) *anchorIS = plex->anchorIS;
5757: return(0);
5758: }
5762: /*@
5763: DMPlexSetAnchors - Set the layout of the local anchor (point-to-point) constraints. Unlike boundary conditions,
5764: when a point's degrees of freedom in a section are constrained to an outside value, the anchor constraints set a
5765: point's degrees of freedom to be a linear combination of other points' degrees of freedom.
5767: After specifying the layout of constraints with DMPlexSetAnchors(), one specifies the constraints by calling
5768: DMGetConstraints() and filling in the entries in the constraint matrix.
5770: collective on dm
5772: Input Parameters:
5773: + dm - The DMPlex object
5774: . anchorSection - The section that describes the mapping from constrained points to the anchor points listed in anchorIS. Must have a local communicator (PETSC_COMM_SELF or derivative).
5775: - anchorIS - The list of all anchor points. Must have a local communicator (PETSC_COMM_SELF or derivative).
5777: The reference counts of anchorSection and anchorIS are incremented.
5779: Level: intermediate
5781: .seealso: DMPlexGetAnchors(), DMGetConstraints(), DMSetConstraints()
5782: @*/
5783: PetscErrorCode DMPlexSetAnchors(DM dm, PetscSection anchorSection, IS anchorIS)
5784: {
5785: DM_Plex *plex = (DM_Plex *)dm->data;
5786: PetscMPIInt result;
5791: if (anchorSection) {
5793: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorSection),&result);
5794: if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor section must have local communicator");
5795: }
5796: if (anchorIS) {
5798: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)anchorIS),&result);
5799: if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"anchor IS must have local communicator");
5800: }
5802: PetscObjectReference((PetscObject)anchorSection);
5803: PetscSectionDestroy(&plex->anchorSection);
5804: plex->anchorSection = anchorSection;
5806: PetscObjectReference((PetscObject)anchorIS);
5807: ISDestroy(&plex->anchorIS);
5808: plex->anchorIS = anchorIS;
5810: #if defined(PETSC_USE_DEBUG)
5811: if (anchorIS && anchorSection) {
5812: PetscInt size, a, pStart, pEnd;
5813: const PetscInt *anchors;
5815: PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5816: ISGetLocalSize(anchorIS,&size);
5817: ISGetIndices(anchorIS,&anchors);
5818: for (a = 0; a < size; a++) {
5819: PetscInt p;
5821: p = anchors[a];
5822: if (p >= pStart && p < pEnd) {
5823: PetscInt dof;
5825: PetscSectionGetDof(anchorSection,p,&dof);
5826: if (dof) {
5827: PetscErrorCode ierr2;
5829: ierr2 = ISRestoreIndices(anchorIS,&anchors);CHKERRQ(ierr2);
5830: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point %d cannot be constrained and an anchor",p);
5831: }
5832: }
5833: }
5834: ISRestoreIndices(anchorIS,&anchors);
5835: }
5836: #endif
5837: /* reset the generic constraints */
5838: DMSetDefaultConstraints(dm,NULL,NULL);
5839: return(0);
5840: }
5844: static PetscErrorCode DMPlexCreateConstraintSection_Anchors(DM dm, PetscSection section, PetscSection *cSec)
5845: {
5846: PetscSection anchorSection;
5847: PetscInt pStart, pEnd, p, dof, numFields, f;
5852: DMPlexGetAnchors(dm,&anchorSection,NULL);
5853: PetscSectionCreate(PETSC_COMM_SELF,cSec);
5854: PetscSectionGetNumFields(section,&numFields);
5855: PetscSectionSetNumFields(*cSec,numFields);
5856: PetscSectionGetChart(anchorSection,&pStart,&pEnd);
5857: PetscSectionSetChart(*cSec,pStart,pEnd);
5858: for (p = pStart; p < pEnd; p++) {
5859: PetscSectionGetDof(anchorSection,p,&dof);
5860: if (dof) {
5861: PetscSectionGetDof(section,p,&dof);
5862: PetscSectionSetDof(*cSec,p,dof);
5863: for (f = 0; f < numFields; f++) {
5864: PetscSectionGetFieldDof(section,p,f,&dof);
5865: PetscSectionSetFieldDof(*cSec,p,f,dof);
5866: }
5867: }
5868: }
5869: PetscSectionSetUp(*cSec);
5870: return(0);
5871: }
5875: static PetscErrorCode DMPlexCreateConstraintMatrix_Anchors(DM dm, PetscSection section, PetscSection cSec, Mat *cMat)
5876: {
5877: PetscSection aSec;
5878: PetscInt pStart, pEnd, p, dof, aDof, aOff, off, nnz, annz, m, n, q, a, offset, *i, *j;
5879: const PetscInt *anchors;
5880: PetscInt numFields, f;
5881: IS aIS;
5886: PetscSectionGetStorageSize(cSec, &m);
5887: PetscSectionGetStorageSize(section, &n);
5888: MatCreate(PETSC_COMM_SELF,cMat);
5889: MatSetSizes(*cMat,m,n,m,n);
5890: MatSetType(*cMat,MATSEQAIJ);
5891: DMPlexGetAnchors(dm,&aSec,&aIS);
5892: ISGetIndices(aIS,&anchors);
5893: PetscSectionGetChart(aSec,&pStart,&pEnd);
5894: PetscMalloc1(m+1,&i);
5895: i[0] = 0;
5896: PetscSectionGetNumFields(section,&numFields);
5897: for (p = pStart; p < pEnd; p++) {
5898: PetscSectionGetDof(aSec,p,&dof);
5899: if (!dof) continue;
5900: PetscSectionGetOffset(aSec,p,&off);
5901: if (numFields) {
5902: for (f = 0; f < numFields; f++) {
5903: annz = 0;
5904: for (q = 0; q < dof; q++) {
5905: a = anchors[off + q];
5906: PetscSectionGetFieldDof(section,a,f,&aDof);
5907: annz += aDof;
5908: }
5909: PetscSectionGetFieldDof(cSec,p,f,&dof);
5910: PetscSectionGetFieldOffset(cSec,p,f,&off);
5911: for (q = 0; q < dof; q++) {
5912: i[off + q + 1] = i[off + q] + annz;
5913: }
5914: }
5915: }
5916: else {
5917: annz = 0;
5918: for (q = 0; q < dof; q++) {
5919: a = anchors[off + q];
5920: PetscSectionGetDof(section,a,&aDof);
5921: annz += aDof;
5922: }
5923: PetscSectionGetDof(cSec,p,&dof);
5924: PetscSectionGetOffset(cSec,p,&off);
5925: for (q = 0; q < dof; q++) {
5926: i[off + q + 1] = i[off + q] + annz;
5927: }
5928: }
5929: }
5930: nnz = i[m];
5931: PetscMalloc1(nnz,&j);
5932: offset = 0;
5933: for (p = pStart; p < pEnd; p++) {
5934: if (numFields) {
5935: for (f = 0; f < numFields; f++) {
5936: PetscSectionGetFieldDof(cSec,p,f,&dof);
5937: for (q = 0; q < dof; q++) {
5938: PetscInt rDof, rOff, r;
5939: PetscSectionGetDof(aSec,p,&rDof);
5940: PetscSectionGetOffset(aSec,p,&rOff);
5941: for (r = 0; r < rDof; r++) {
5942: PetscInt s;
5944: a = anchors[rOff + r];
5946: PetscSectionGetFieldDof(section,a,f,&aDof);
5947: PetscSectionGetFieldOffset(section,a,f,&aOff);
5948: for (s = 0; s < aDof; s++) {
5949: j[offset++] = aOff + s;
5950: }
5951: }
5952: }
5953: }
5954: }
5955: else {
5956: PetscSectionGetDof(cSec,p,&dof);
5957: for (q = 0; q < dof; q++) {
5958: PetscInt rDof, rOff, r;
5959: PetscSectionGetDof(aSec,p,&rDof);
5960: PetscSectionGetOffset(aSec,p,&rOff);
5961: for (r = 0; r < rDof; r++) {
5962: PetscInt s;
5964: a = anchors[rOff + r];
5966: PetscSectionGetDof(section,a,&aDof);
5967: PetscSectionGetOffset(section,a,&aOff);
5968: for (s = 0; s < aDof; s++) {
5969: j[offset++] = aOff + s;
5970: }
5971: }
5972: }
5973: }
5974: }
5975: MatSeqAIJSetPreallocationCSR(*cMat,i,j,NULL);
5976: PetscFree(i);
5977: PetscFree(j);
5978: ISRestoreIndices(aIS,&anchors);
5979: return(0);
5980: }
5984: PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm)
5985: {
5986: DM_Plex *plex = (DM_Plex *)dm->data;
5987: PetscSection anchorSection, section, cSec;
5988: Mat cMat;
5993: DMPlexGetAnchors(dm,&anchorSection,NULL);
5994: if (anchorSection) {
5995: PetscDS ds;
5996: PetscInt nf;
5998: DMGetDefaultSection(dm,§ion);
5999: DMPlexCreateConstraintSection_Anchors(dm,section,&cSec);
6000: DMPlexCreateConstraintMatrix_Anchors(dm,section,cSec,&cMat);
6001: DMGetDS(dm,&ds);
6002: PetscDSGetNumFields(ds,&nf);
6003: if (nf && plex->computeanchormatrix) {(*plex->computeanchormatrix)(dm,section,cSec,cMat);}
6004: DMSetDefaultConstraints(dm,cSec,cMat);
6005: PetscSectionDestroy(&cSec);
6006: MatDestroy(&cMat);
6007: }
6008: return(0);
6009: }