Actual source code: plex.c
petsc-3.5.4 2015-05-23
1: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <../src/sys/utils/hash.h>
3: #include <petsc-private/isimpl.h>
4: #include <petscsf.h>
6: /* Logging support */
7: PetscLogEvent DMPLEX_Interpolate, DMPLEX_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeField, DMPLEX_DistributeData, DMPLEX_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM, DMPLEX_InterpolatorFEM, DMPLEX_InjectorFEM;
9: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
10: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
11: PETSC_EXTERN PetscErrorCode VecLoad_Default(Vec, PetscViewer);
15: PetscErrorCode DMPlexGetFieldType_Internal(DM dm, PetscSection section, PetscInt field, PetscInt *sStart, PetscInt *sEnd, PetscViewerVTKFieldType *ft)
16: {
17: PetscInt dim, pStart, pEnd, vStart, vEnd, cStart, cEnd, vdof = 0, cdof = 0;
21: *ft = PETSC_VTK_POINT_FIELD;
22: DMPlexGetDimension(dm, &dim);
23: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
24: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
25: PetscSectionGetChart(section, &pStart, &pEnd);
26: if (field >= 0) {
27: if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetFieldDof(section, vStart, field, &vdof);}
28: if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetFieldDof(section, cStart, field, &cdof);}
29: } else {
30: if ((vStart >= pStart) && (vStart < pEnd)) {PetscSectionGetDof(section, vStart, &vdof);}
31: if ((cStart >= pStart) && (cStart < pEnd)) {PetscSectionGetDof(section, cStart, &cdof);}
32: }
33: if (vdof) {
34: *sStart = vStart;
35: *sEnd = vEnd;
36: if (vdof == dim) *ft = PETSC_VTK_POINT_VECTOR_FIELD;
37: else *ft = PETSC_VTK_POINT_FIELD;
38: } else if (cdof) {
39: *sStart = cStart;
40: *sEnd = cEnd;
41: if (cdof == dim) *ft = PETSC_VTK_CELL_VECTOR_FIELD;
42: else *ft = PETSC_VTK_CELL_FIELD;
43: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Could not classify input Vec for VTK");
44: return(0);
45: }
49: PetscErrorCode VecView_Plex_Local(Vec v, PetscViewer viewer)
50: {
51: DM dm;
52: PetscBool isvtk, ishdf5, isseq;
56: VecGetDM(v, &dm);
57: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
58: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
59: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
60: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
61: if (isvtk || ishdf5) {
62: PetscInt numFields;
63: PetscBool fem = PETSC_FALSE;
65: DMGetNumFields(dm, &numFields);
66: if (numFields) {
67: PetscObject fe;
69: DMGetField(dm, 0, &fe);
70: if (fe->classid == PETSCFE_CLASSID) fem = PETSC_TRUE;
71: }
72: if (fem) {DMPlexInsertBoundaryValuesFEM(dm, v);}
73: }
74: if (isvtk) {
75: PetscSection section;
76: PetscViewerVTKFieldType ft;
77: PetscInt pStart, pEnd;
79: DMGetDefaultSection(dm, §ion);
80: DMPlexGetFieldType_Internal(dm, section, PETSC_DETERMINE, &pStart, &pEnd, &ft);
81: PetscObjectReference((PetscObject) dm); /* viewer drops reference */
82: PetscObjectReference((PetscObject) v); /* viewer drops reference */
83: PetscViewerVTKAddField(viewer, (PetscObject) dm, DMPlexVTKWriteAll, ft, (PetscObject) v);
84: } else if (ishdf5) {
85: #if defined(PETSC_HAVE_HDF5)
86: VecView_Plex_Local_HDF5(v, viewer);
87: #else
88: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
89: #endif
90: } else {
91: if (isseq) {VecView_Seq(v, viewer);}
92: else {VecView_MPI(v, viewer);}
93: }
94: return(0);
95: }
99: PetscErrorCode VecView_Plex(Vec v, PetscViewer viewer)
100: {
101: DM dm;
102: PetscBool isvtk, ishdf5, isseq;
106: VecGetDM(v, &dm);
107: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
108: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
109: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
110: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
111: if (isvtk) {
112: Vec locv;
113: const char *name;
115: DMGetLocalVector(dm, &locv);
116: PetscObjectGetName((PetscObject) v, &name);
117: PetscObjectSetName((PetscObject) locv, name);
118: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
119: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
120: VecView_Plex_Local(locv, viewer);
121: DMRestoreLocalVector(dm, &locv);
122: } else if (ishdf5) {
123: #if defined(PETSC_HAVE_HDF5)
124: VecView_Plex_HDF5(v, viewer);
125: #else
126: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
127: #endif
128: } else {
129: if (isseq) {VecView_Seq(v, viewer);}
130: else {VecView_MPI(v, viewer);}
131: }
132: return(0);
133: }
137: PetscErrorCode VecLoad_Plex_Local(Vec v, PetscViewer viewer)
138: {
139: DM dm;
140: PetscBool ishdf5;
144: VecGetDM(v, &dm);
145: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
146: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
147: if (ishdf5) {
148: DM dmBC;
149: Vec gv;
150: const char *name;
152: DMGetOutputDM(dm, &dmBC);
153: DMGetGlobalVector(dmBC, &gv);
154: PetscObjectGetName((PetscObject) v, &name);
155: PetscObjectSetName((PetscObject) gv, name);
156: VecLoad_Default(gv, viewer);
157: DMGlobalToLocalBegin(dmBC, gv, INSERT_VALUES, v);
158: DMGlobalToLocalEnd(dmBC, gv, INSERT_VALUES, v);
159: DMRestoreGlobalVector(dmBC, &gv);
160: } else {
161: VecLoad_Default(v, viewer);
162: }
163: return(0);
164: }
168: PetscErrorCode VecLoad_Plex(Vec v, PetscViewer viewer)
169: {
170: DM dm;
171: PetscBool ishdf5;
175: VecGetDM(v, &dm);
176: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
177: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
178: if (ishdf5) {
179: #if defined(PETSC_HAVE_HDF5)
180: VecLoad_Plex_HDF5(v, viewer);
181: #else
182: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
183: #endif
184: } else {
185: VecLoad_Default(v, viewer);
186: }
187: return(0);
188: }
192: PetscErrorCode DMPlexView_Ascii(DM dm, PetscViewer viewer)
193: {
194: DM_Plex *mesh = (DM_Plex*) dm->data;
195: DM cdm;
196: DMLabel markers;
197: PetscSection coordSection;
198: Vec coordinates;
199: PetscViewerFormat format;
200: PetscErrorCode ierr;
203: DMGetCoordinateDM(dm, &cdm);
204: DMGetDefaultSection(cdm, &coordSection);
205: DMGetCoordinatesLocal(dm, &coordinates);
206: PetscViewerGetFormat(viewer, &format);
207: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
208: const char *name;
209: PetscInt maxConeSize, maxSupportSize;
210: PetscInt pStart, pEnd, p;
211: PetscMPIInt rank, size;
213: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
214: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
215: PetscObjectGetName((PetscObject) dm, &name);
216: DMPlexGetChart(dm, &pStart, &pEnd);
217: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
218: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
219: PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
220: PetscViewerASCIISynchronizedPrintf(viewer, "Max sizes cone: %D support: %D\n", maxConeSize, maxSupportSize);
221: PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
222: PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
223: for (p = pStart; p < pEnd; ++p) {
224: PetscInt dof, off, s;
226: PetscSectionGetDof(mesh->supportSection, p, &dof);
227: PetscSectionGetOffset(mesh->supportSection, p, &off);
228: for (s = off; s < off+dof; ++s) {
229: PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D ----> %D\n", rank, p, mesh->supports[s]);
230: }
231: }
232: PetscViewerFlush(viewer);
233: PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
234: for (p = pStart; p < pEnd; ++p) {
235: PetscInt dof, off, c;
237: PetscSectionGetDof(mesh->coneSection, p, &dof);
238: PetscSectionGetOffset(mesh->coneSection, p, &off);
239: for (c = off; c < off+dof; ++c) {
240: PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D <---- %D (%D)\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]);
241: }
242: }
243: PetscViewerFlush(viewer);
244: PetscSectionGetChart(coordSection, &pStart, NULL);
245: if (pStart >= 0) {PetscSectionVecView(coordSection, coordinates, viewer);}
246: DMPlexGetLabel(dm, "marker", &markers);
247: DMLabelView(markers,viewer);
248: if (size > 1) {
249: PetscSF sf;
251: DMGetPointSF(dm, &sf);
252: PetscSFView(sf, viewer);
253: }
254: PetscViewerFlush(viewer);
255: } else if (format == PETSC_VIEWER_ASCII_LATEX) {
256: const char *name;
257: const char *colors[3] = {"red", "blue", "green"};
258: const int numColors = 3;
259: PetscReal scale = 2.0;
260: PetscScalar *coords;
261: PetscInt depth, cStart, cEnd, c, vStart, vEnd, v, eStart = 0, eEnd = 0, e, p;
262: PetscMPIInt rank, size;
264: DMPlexGetDepth(dm, &depth);
265: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
266: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
267: PetscObjectGetName((PetscObject) dm, &name);
268: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
269: PetscViewerASCIIPrintf(viewer, "\
270: \\documentclass[crop,multi=false]{standalone}\n\n\
271: \\usepackage{tikz}\n\
272: \\usepackage{pgflibraryshapes}\n\
273: \\usetikzlibrary{backgrounds}\n\
274: \\usetikzlibrary{arrows}\n\
275: \\begin{document}\n\
276: \\section{%s}\n\
277: \\begin{center}\n", name, 8.0/scale);
278: PetscViewerASCIIPrintf(viewer, "Mesh for process ");
279: for (p = 0; p < size; ++p) {
280: if (p > 0 && p == size-1) {
281: PetscViewerASCIIPrintf(viewer, ", and ", colors[p%numColors], p);
282: } else if (p > 0) {
283: PetscViewerASCIIPrintf(viewer, ", ", colors[p%numColors], p);
284: }
285: PetscViewerASCIIPrintf(viewer, "{\\textcolor{%s}%D}", colors[p%numColors], p);
286: }
287: PetscViewerASCIIPrintf(viewer, ".\n\n\n\
288: \\begin{tikzpicture}[scale = %g,font=\\fontsize{8}{8}\\selectfont]\n");
289: /* Plot vertices */
290: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
291: VecGetArray(coordinates, &coords);
292: PetscViewerASCIIPrintf(viewer, "\\path\n");
293: for (v = vStart; v < vEnd; ++v) {
294: PetscInt off, dof, d;
296: PetscSectionGetDof(coordSection, v, &dof);
297: PetscSectionGetOffset(coordSection, v, &off);
298: PetscViewerASCIISynchronizedPrintf(viewer, "(");
299: for (d = 0; d < dof; ++d) {
300: if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
301: PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)(scale*PetscRealPart(coords[off+d])));
302: }
303: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [draw,shape=circle,color=%s] {%D} --\n", v, rank, colors[rank%numColors], v);
304: }
305: VecRestoreArray(coordinates, &coords);
306: PetscViewerFlush(viewer);
307: PetscViewerASCIIPrintf(viewer, "(0,0);\n");
308: /* Plot edges */
309: VecGetArray(coordinates, &coords);
310: PetscViewerASCIIPrintf(viewer, "\\path\n");
311: if (depth > 1) {DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);}
312: for (e = eStart; e < eEnd; ++e) {
313: const PetscInt *cone;
314: PetscInt coneSize, offA, offB, dof, d;
316: DMPlexGetConeSize(dm, e, &coneSize);
317: if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %d cone should have two vertices, not %d", e, coneSize);
318: DMPlexGetCone(dm, e, &cone);
319: PetscSectionGetDof(coordSection, cone[0], &dof);
320: PetscSectionGetOffset(coordSection, cone[0], &offA);
321: PetscSectionGetOffset(coordSection, cone[1], &offB);
322: PetscViewerASCIISynchronizedPrintf(viewer, "(");
323: for (d = 0; d < dof; ++d) {
324: if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
325: PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)(scale*0.5*PetscRealPart(coords[offA+d]+coords[offB+d])));
326: }
327: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [draw,shape=circle,color=%s] {%D} --\n", e, rank, colors[rank%numColors], e);
328: }
329: VecRestoreArray(coordinates, &coords);
330: PetscViewerFlush(viewer);
331: PetscViewerASCIIPrintf(viewer, "(0,0);\n");
332: /* Plot cells */
333: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
334: for (c = cStart; c < cEnd; ++c) {
335: PetscInt *closure = NULL;
336: PetscInt closureSize, firstPoint = -1;
338: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
339: PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] ", colors[rank%numColors]);
340: for (p = 0; p < closureSize*2; p += 2) {
341: const PetscInt point = closure[p];
343: if ((point < vStart) || (point >= vEnd)) continue;
344: if (firstPoint >= 0) {PetscViewerASCIISynchronizedPrintf(viewer, " -- ");}
345: PetscViewerASCIISynchronizedPrintf(viewer, "(%D_%D)", point, rank);
346: if (firstPoint < 0) firstPoint = point;
347: }
348: /* Why doesn't this work? PetscViewerASCIISynchronizedPrintf(viewer, " -- cycle;\n"); */
349: PetscViewerASCIISynchronizedPrintf(viewer, " -- (%D_%D);\n", firstPoint, rank);
350: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
351: }
352: PetscViewerFlush(viewer);
353: PetscViewerASCIIPrintf(viewer, "\\end{tikzpicture}\n\\end{center}\n");
354: PetscViewerASCIIPrintf(viewer, "\\end{document}\n", name);
355: } else {
356: MPI_Comm comm;
357: PetscInt *sizes, *hybsizes;
358: PetscInt locDepth, depth, dim, d, pMax[4];
359: PetscInt pStart, pEnd, p;
360: PetscInt numLabels, l;
361: const char *name;
362: PetscMPIInt size;
364: PetscObjectGetComm((PetscObject)dm,&comm);
365: MPI_Comm_size(comm, &size);
366: DMPlexGetDimension(dm, &dim);
367: PetscObjectGetName((PetscObject) dm, &name);
368: if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dim);}
369: else {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dim);}
370: DMPlexGetDepth(dm, &locDepth);
371: MPI_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
372: DMPlexGetHybridBounds(dm, &pMax[depth], depth > 0 ? &pMax[depth-1] : NULL, &pMax[1], &pMax[0]);
373: PetscMalloc2(size,&sizes,size,&hybsizes);
374: if (depth == 1) {
375: DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd);
376: pEnd = pEnd - pStart;
377: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
378: PetscViewerASCIIPrintf(viewer, " %D-cells:", 0);
379: for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
380: PetscViewerASCIIPrintf(viewer, "\n");
381: DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd);
382: pEnd = pEnd - pStart;
383: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
384: PetscViewerASCIIPrintf(viewer, " %D-cells:", dim);
385: for (p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
386: PetscViewerASCIIPrintf(viewer, "\n");
387: } else {
388: for (d = 0; d <= dim; d++) {
389: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
390: pEnd -= pStart;
391: pMax[d] -= pStart;
392: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
393: MPI_Gather(&pMax[d], 1, MPIU_INT, hybsizes, 1, MPIU_INT, 0, comm);
394: PetscViewerASCIIPrintf(viewer, " %D-cells:", d);
395: for (p = 0; p < size; ++p) {
396: if (hybsizes[p] >= 0) {PetscViewerASCIIPrintf(viewer, " %D (%D)", sizes[p], sizes[p] - hybsizes[p]);}
397: else {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
398: }
399: PetscViewerASCIIPrintf(viewer, "\n");
400: }
401: }
402: PetscFree2(sizes,hybsizes);
403: DMPlexGetNumLabels(dm, &numLabels);
404: if (numLabels) {PetscViewerASCIIPrintf(viewer, "Labels:\n");}
405: for (l = 0; l < numLabels; ++l) {
406: DMLabel label;
407: const char *name;
408: IS valueIS;
409: const PetscInt *values;
410: PetscInt numValues, v;
412: DMPlexGetLabelName(dm, l, &name);
413: DMPlexGetLabel(dm, name, &label);
414: DMLabelGetNumValues(label, &numValues);
415: PetscViewerASCIIPrintf(viewer, " %s: %d strata of sizes (", name, numValues);
416: DMLabelGetValueIS(label, &valueIS);
417: ISGetIndices(valueIS, &values);
418: for (v = 0; v < numValues; ++v) {
419: PetscInt size;
421: DMLabelGetStratumSize(label, values[v], &size);
422: if (v > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
423: PetscViewerASCIIPrintf(viewer, "%d", size);
424: }
425: PetscViewerASCIIPrintf(viewer, ")\n");
426: ISRestoreIndices(valueIS, &values);
427: ISDestroy(&valueIS);
428: }
429: }
430: return(0);
431: }
435: PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
436: {
437: PetscBool iascii, ishdf5;
443: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
444: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
445: if (iascii) {
446: DMPlexView_Ascii(dm, viewer);
447: } else if (ishdf5) {
448: #if defined(PETSC_HAVE_HDF5)
449: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
450: DMPlexView_HDF5(dm, viewer);
451: PetscViewerPopFormat(viewer);
452: #else
453: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
454: #endif
455: }
456: return(0);
457: }
461: PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
462: {
463: PetscBool isbinary, ishdf5;
469: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
470: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
471: if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");}
472: else if (ishdf5) {
473: #if defined(PETSC_HAVE_HDF5)
474: DMPlexLoad_HDF5(dm, viewer);
475: #else
476: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
477: #endif
478: }
479: return(0);
480: }
484: static PetscErrorCode BoundaryDestroy(DMBoundary *boundary)
485: {
486: DMBoundary b, next;
490: if (!boundary) return(0);
491: b = *boundary;
492: *boundary = NULL;
493: for (; b; b = next) {
494: next = b->next;
495: PetscFree(b->ids);
496: PetscFree(b->name);
497: PetscFree(b->labelname);
498: PetscFree(b);
499: }
500: return(0);
501: }
505: PetscErrorCode DMDestroy_Plex(DM dm)
506: {
507: DM_Plex *mesh = (DM_Plex*) dm->data;
508: DMLabel next = mesh->labels;
512: if (--mesh->refct > 0) return(0);
513: PetscSectionDestroy(&mesh->coneSection);
514: PetscFree(mesh->cones);
515: PetscFree(mesh->coneOrientations);
516: PetscSectionDestroy(&mesh->supportSection);
517: PetscFree(mesh->supports);
518: PetscFree(mesh->facesTmp);
519: while (next) {
520: DMLabel tmp = next->next;
522: DMLabelDestroy(&next);
523: next = tmp;
524: }
525: DMDestroy(&mesh->coarseMesh);
526: DMLabelDestroy(&mesh->subpointMap);
527: ISDestroy(&mesh->globalVertexNumbers);
528: ISDestroy(&mesh->globalCellNumbers);
529: BoundaryDestroy(&mesh->boundary);
530: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
531: PetscFree(mesh);
532: return(0);
533: }
537: PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
538: {
539: PetscSection section, sectionGlobal;
540: PetscInt bs = -1;
541: PetscInt localSize;
542: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock;
544: MatType mtype;
547: MatInitializePackage();
548: mtype = dm->mattype;
549: DMGetDefaultSection(dm, §ion);
550: DMGetDefaultGlobalSection(dm, §ionGlobal);
551: /* PetscSectionGetStorageSize(sectionGlobal, &localSize); */
552: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
553: MatCreate(PetscObjectComm((PetscObject)dm), J);
554: MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
555: MatSetType(*J, mtype);
556: MatSetFromOptions(*J);
557: PetscStrcmp(mtype, MATSHELL, &isShell);
558: PetscStrcmp(mtype, MATBAIJ, &isBlock);
559: PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
560: PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
561: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
562: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
563: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
564: if (!isShell) {
565: PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
566: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin;
568: if (bs < 0) {
569: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
570: PetscInt pStart, pEnd, p, dof, cdof;
572: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
573: for (p = pStart; p < pEnd; ++p) {
574: PetscSectionGetDof(sectionGlobal, p, &dof);
575: PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
576: if (dof-cdof) {
577: if (bs < 0) {
578: bs = dof-cdof;
579: } else if (bs != dof-cdof) {
580: /* Layout does not admit a pointwise block size */
581: bs = 1;
582: break;
583: }
584: }
585: }
586: /* Must have same blocksize on all procs (some might have no points) */
587: bsLocal = bs;
588: MPI_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
589: bsLocal = bs < 0 ? bsMax : bs;
590: MPI_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
591: if (bsMin != bsMax) {
592: bs = 1;
593: } else {
594: bs = bsMax;
595: }
596: } else {
597: bs = 1;
598: }
599: }
600: PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);
601: DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix);
602: PetscFree4(dnz, onz, dnzu, onzu);
603: }
604: return(0);
605: }
609: /*@
610: DMPlexGetDimension - Return the topological mesh dimension
612: Not collective
614: Input Parameter:
615: . mesh - The DMPlex
617: Output Parameter:
618: . dim - The topological mesh dimension
620: Level: beginner
622: .seealso: DMPlexCreate()
623: @*/
624: PetscErrorCode DMPlexGetDimension(DM dm, PetscInt *dim)
625: {
626: DM_Plex *mesh = (DM_Plex*) dm->data;
631: *dim = mesh->dim;
632: return(0);
633: }
637: /*@
638: DMPlexSetDimension - Set the topological mesh dimension
640: Collective on mesh
642: Input Parameters:
643: + mesh - The DMPlex
644: - dim - The topological mesh dimension
646: Level: beginner
648: .seealso: DMPlexCreate()
649: @*/
650: PetscErrorCode DMPlexSetDimension(DM dm, PetscInt dim)
651: {
652: DM_Plex *mesh = (DM_Plex*) dm->data;
657: mesh->dim = dim;
658: return(0);
659: }
663: /*@
664: DMPlexGetChart - Return the interval for all mesh points [pStart, pEnd)
666: Not collective
668: Input Parameter:
669: . mesh - The DMPlex
671: Output Parameters:
672: + pStart - The first mesh point
673: - pEnd - The upper bound for mesh points
675: Level: beginner
677: .seealso: DMPlexCreate(), DMPlexSetChart()
678: @*/
679: PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
680: {
681: DM_Plex *mesh = (DM_Plex*) dm->data;
686: PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
687: return(0);
688: }
692: /*@
693: DMPlexSetChart - Set the interval for all mesh points [pStart, pEnd)
695: Not collective
697: Input Parameters:
698: + mesh - The DMPlex
699: . pStart - The first mesh point
700: - pEnd - The upper bound for mesh points
702: Output Parameters:
704: Level: beginner
706: .seealso: DMPlexCreate(), DMPlexGetChart()
707: @*/
708: PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
709: {
710: DM_Plex *mesh = (DM_Plex*) dm->data;
715: PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
716: PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
717: return(0);
718: }
722: /*@
723: DMPlexGetConeSize - Return the number of in-edges for this point in the Sieve DAG
725: Not collective
727: Input Parameters:
728: + mesh - The DMPlex
729: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
731: Output Parameter:
732: . size - The cone size for point p
734: Level: beginner
736: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
737: @*/
738: PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size)
739: {
740: DM_Plex *mesh = (DM_Plex*) dm->data;
746: PetscSectionGetDof(mesh->coneSection, p, size);
747: return(0);
748: }
752: /*@
753: DMPlexSetConeSize - Set the number of in-edges for this point in the Sieve DAG
755: Not collective
757: Input Parameters:
758: + mesh - The DMPlex
759: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
760: - size - The cone size for point p
762: Output Parameter:
764: Note:
765: This should be called after DMPlexSetChart().
767: Level: beginner
769: .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart()
770: @*/
771: PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size)
772: {
773: DM_Plex *mesh = (DM_Plex*) dm->data;
778: PetscSectionSetDof(mesh->coneSection, p, size);
780: mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
781: return(0);
782: }
786: /*@
787: DMPlexAddConeSize - Add the given number of in-edges to this point in the Sieve DAG
789: Not collective
791: Input Parameters:
792: + mesh - The DMPlex
793: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
794: - size - The additional cone size for point p
796: Output Parameter:
798: Note:
799: This should be called after DMPlexSetChart().
801: Level: beginner
803: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart()
804: @*/
805: PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size)
806: {
807: DM_Plex *mesh = (DM_Plex*) dm->data;
808: PetscInt csize;
813: PetscSectionAddDof(mesh->coneSection, p, size);
814: PetscSectionGetDof(mesh->coneSection, p, &csize);
816: mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize);
817: return(0);
818: }
822: /*@C
823: DMPlexGetCone - Return the points on the in-edges for this point in the Sieve DAG
825: Not collective
827: Input Parameters:
828: + mesh - The DMPlex
829: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
831: Output Parameter:
832: . cone - An array of points which are on the in-edges for point p
834: Level: beginner
836: Fortran Notes:
837: Since it returns an array, this routine is only available in Fortran 90, and you must
838: include petsc.h90 in your code.
840: You must also call DMPlexRestoreCone() after you finish using the returned array.
842: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart()
843: @*/
844: PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
845: {
846: DM_Plex *mesh = (DM_Plex*) dm->data;
847: PetscInt off;
853: PetscSectionGetOffset(mesh->coneSection, p, &off);
854: *cone = &mesh->cones[off];
855: return(0);
856: }
860: /*@
861: DMPlexSetCone - Set the points on the in-edges for this point in the Sieve DAG
863: Not collective
865: Input Parameters:
866: + mesh - The DMPlex
867: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
868: - cone - An array of points which are on the in-edges for point p
870: Output Parameter:
872: Note:
873: This should be called after all calls to DMPlexSetConeSize() and DMSetUp().
875: Level: beginner
877: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
878: @*/
879: PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[])
880: {
881: DM_Plex *mesh = (DM_Plex*) dm->data;
882: PetscInt pStart, pEnd;
883: PetscInt dof, off, c;
888: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
889: PetscSectionGetDof(mesh->coneSection, p, &dof);
891: PetscSectionGetOffset(mesh->coneSection, p, &off);
892: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
893: for (c = 0; c < dof; ++c) {
894: if ((cone[c] < pStart) || (cone[c] >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D, %D)", cone[c], pStart, pEnd);
895: mesh->cones[off+c] = cone[c];
896: }
897: return(0);
898: }
902: /*@C
903: DMPlexGetConeOrientation - Return the orientations on the in-edges for this point in the Sieve DAG
905: Not collective
907: Input Parameters:
908: + mesh - The DMPlex
909: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
911: Output Parameter:
912: . coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
913: integer giving the prescription for cone traversal. If it is negative, the cone is
914: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
915: the index of the cone point on which to start.
917: Level: beginner
919: Fortran Notes:
920: Since it returns an array, this routine is only available in Fortran 90, and you must
921: include petsc.h90 in your code.
923: You must also call DMPlexRestoreConeOrientation() after you finish using the returned array.
925: .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
926: @*/
927: PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
928: {
929: DM_Plex *mesh = (DM_Plex*) dm->data;
930: PetscInt off;
935: #if defined(PETSC_USE_DEBUG)
936: {
937: PetscInt dof;
938: PetscSectionGetDof(mesh->coneSection, p, &dof);
940: }
941: #endif
942: PetscSectionGetOffset(mesh->coneSection, p, &off);
944: *coneOrientation = &mesh->coneOrientations[off];
945: return(0);
946: }
950: /*@
951: DMPlexSetConeOrientation - Set the orientations on the in-edges for this point in the Sieve DAG
953: Not collective
955: Input Parameters:
956: + mesh - The DMPlex
957: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
958: - coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
959: integer giving the prescription for cone traversal. If it is negative, the cone is
960: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
961: the index of the cone point on which to start.
963: Output Parameter:
965: Note:
966: This should be called after all calls to DMPlexSetConeSize() and DMSetUp().
968: Level: beginner
970: .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
971: @*/
972: PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
973: {
974: DM_Plex *mesh = (DM_Plex*) dm->data;
975: PetscInt pStart, pEnd;
976: PetscInt dof, off, c;
981: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
982: PetscSectionGetDof(mesh->coneSection, p, &dof);
984: PetscSectionGetOffset(mesh->coneSection, p, &off);
985: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
986: for (c = 0; c < dof; ++c) {
987: PetscInt cdof, o = coneOrientation[c];
989: PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
990: if (o && ((o < -(cdof+1)) || (o >= cdof))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone orientation %D is not in the valid range [%D. %D)", o, -(cdof+1), cdof);
991: mesh->coneOrientations[off+c] = o;
992: }
993: return(0);
994: }
998: PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
999: {
1000: DM_Plex *mesh = (DM_Plex*) dm->data;
1001: PetscInt pStart, pEnd;
1002: PetscInt dof, off;
1007: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1008: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1009: if ((conePoint < pStart) || (conePoint >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D, %D)", conePoint, pStart, pEnd);
1010: PetscSectionGetDof(mesh->coneSection, p, &dof);
1011: PetscSectionGetOffset(mesh->coneSection, p, &off);
1012: if ((conePos < 0) || (conePos >= dof)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone position %D of point %D is not in the valid range [0, %D)", conePos, p, dof);
1013: mesh->cones[off+conePos] = conePoint;
1014: return(0);
1015: }
1019: PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation)
1020: {
1021: DM_Plex *mesh = (DM_Plex*) dm->data;
1022: PetscInt pStart, pEnd;
1023: PetscInt dof, off;
1028: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1029: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1030: PetscSectionGetDof(mesh->coneSection, p, &dof);
1031: PetscSectionGetOffset(mesh->coneSection, p, &off);
1032: if ((conePos < 0) || (conePos >= dof)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone position %D of point %D is not in the valid range [0, %D)", conePos, p, dof);
1033: mesh->coneOrientations[off+conePos] = coneOrientation;
1034: return(0);
1035: }
1039: /*@
1040: DMPlexGetSupportSize - Return the number of out-edges for this point in the Sieve DAG
1042: Not collective
1044: Input Parameters:
1045: + mesh - The DMPlex
1046: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1048: Output Parameter:
1049: . size - The support size for point p
1051: Level: beginner
1053: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize()
1054: @*/
1055: PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1056: {
1057: DM_Plex *mesh = (DM_Plex*) dm->data;
1063: PetscSectionGetDof(mesh->supportSection, p, size);
1064: return(0);
1065: }
1069: /*@
1070: DMPlexSetSupportSize - Set the number of out-edges for this point in the Sieve DAG
1072: Not collective
1074: Input Parameters:
1075: + mesh - The DMPlex
1076: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1077: - size - The support size for point p
1079: Output Parameter:
1081: Note:
1082: This should be called after DMPlexSetChart().
1084: Level: beginner
1086: .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart()
1087: @*/
1088: PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size)
1089: {
1090: DM_Plex *mesh = (DM_Plex*) dm->data;
1095: PetscSectionSetDof(mesh->supportSection, p, size);
1097: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size);
1098: return(0);
1099: }
1103: /*@C
1104: DMPlexGetSupport - Return the points on the out-edges for this point in the Sieve DAG
1106: Not collective
1108: Input Parameters:
1109: + mesh - The DMPlex
1110: - p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1112: Output Parameter:
1113: . support - An array of points which are on the out-edges for point p
1115: Level: beginner
1117: Fortran Notes:
1118: Since it returns an array, this routine is only available in Fortran 90, and you must
1119: include petsc.h90 in your code.
1121: You must also call DMPlexRestoreSupport() after you finish using the returned array.
1123: .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1124: @*/
1125: PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1126: {
1127: DM_Plex *mesh = (DM_Plex*) dm->data;
1128: PetscInt off;
1134: PetscSectionGetOffset(mesh->supportSection, p, &off);
1135: *support = &mesh->supports[off];
1136: return(0);
1137: }
1141: /*@
1142: DMPlexSetSupport - Set the points on the out-edges for this point in the Sieve DAG
1144: Not collective
1146: Input Parameters:
1147: + mesh - The DMPlex
1148: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1149: - support - An array of points which are on the in-edges for point p
1151: Output Parameter:
1153: Note:
1154: This should be called after all calls to DMPlexSetSupportSize() and DMSetUp().
1156: Level: beginner
1158: .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp()
1159: @*/
1160: PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[])
1161: {
1162: DM_Plex *mesh = (DM_Plex*) dm->data;
1163: PetscInt pStart, pEnd;
1164: PetscInt dof, off, c;
1169: PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1170: PetscSectionGetDof(mesh->supportSection, p, &dof);
1172: PetscSectionGetOffset(mesh->supportSection, p, &off);
1173: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1174: for (c = 0; c < dof; ++c) {
1175: if ((support[c] < pStart) || (support[c] >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support point %D is not in the valid range [%D, %D)", support[c], pStart, pEnd);
1176: mesh->supports[off+c] = support[c];
1177: }
1178: return(0);
1179: }
1183: PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint)
1184: {
1185: DM_Plex *mesh = (DM_Plex*) dm->data;
1186: PetscInt pStart, pEnd;
1187: PetscInt dof, off;
1192: PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);
1193: PetscSectionGetDof(mesh->supportSection, p, &dof);
1194: PetscSectionGetOffset(mesh->supportSection, p, &off);
1195: if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1196: if ((supportPoint < pStart) || (supportPoint >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support point %D is not in the valid range [%D, %D)", supportPoint, pStart, pEnd);
1197: if (supportPos >= dof) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support position %D of point %D is not in the valid range [0, %D)", supportPos, p, dof);
1198: mesh->supports[off+supportPos] = supportPoint;
1199: return(0);
1200: }
1204: /*@C
1205: DMPlexGetTransitiveClosure - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1207: Not collective
1209: Input Parameters:
1210: + mesh - The DMPlex
1211: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1212: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1213: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
1215: Output Parameters:
1216: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1217: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
1219: Note:
1220: If using internal storage (points is NULL on input), each call overwrites the last output.
1222: Fortran Notes:
1223: Since it returns an array, this routine is only available in Fortran 90, and you must
1224: include petsc.h90 in your code.
1226: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1228: Level: beginner
1230: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1231: @*/
1232: PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1233: {
1234: DM_Plex *mesh = (DM_Plex*) dm->data;
1235: PetscInt *closure, *fifo;
1236: const PetscInt *tmp = NULL, *tmpO = NULL;
1237: PetscInt tmpSize, t;
1238: PetscInt depth = 0, maxSize;
1239: PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0;
1240: PetscErrorCode ierr;
1244: DMPlexGetDepth(dm, &depth);
1245: /* This is only 1-level */
1246: if (useCone) {
1247: DMPlexGetConeSize(dm, p, &tmpSize);
1248: DMPlexGetCone(dm, p, &tmp);
1249: DMPlexGetConeOrientation(dm, p, &tmpO);
1250: } else {
1251: DMPlexGetSupportSize(dm, p, &tmpSize);
1252: DMPlexGetSupport(dm, p, &tmp);
1253: }
1254: if (depth == 1) {
1255: if (*points) {
1256: closure = *points;
1257: } else {
1258: maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1259: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1260: }
1261: closure[0] = p; closure[1] = 0;
1262: for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1263: closure[closureSize] = tmp[t];
1264: closure[closureSize+1] = tmpO ? tmpO[t] : 0;
1265: }
1266: if (numPoints) *numPoints = closureSize/2;
1267: if (points) *points = closure;
1268: return(0);
1269: }
1270: maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth+1);
1271: DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1272: if (*points) {
1273: closure = *points;
1274: } else {
1275: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1276: }
1277: closure[0] = p; closure[1] = 0;
1278: for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1279: const PetscInt cp = tmp[t];
1280: const PetscInt co = tmpO ? tmpO[t] : 0;
1282: closure[closureSize] = cp;
1283: closure[closureSize+1] = co;
1284: fifo[fifoSize] = cp;
1285: fifo[fifoSize+1] = co;
1286: }
1287: /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1288: while (fifoSize - fifoStart) {
1289: const PetscInt q = fifo[fifoStart];
1290: const PetscInt o = fifo[fifoStart+1];
1291: const PetscInt rev = o >= 0 ? 0 : 1;
1292: const PetscInt off = rev ? -(o+1) : o;
1294: if (useCone) {
1295: DMPlexGetConeSize(dm, q, &tmpSize);
1296: DMPlexGetCone(dm, q, &tmp);
1297: DMPlexGetConeOrientation(dm, q, &tmpO);
1298: } else {
1299: DMPlexGetSupportSize(dm, q, &tmpSize);
1300: DMPlexGetSupport(dm, q, &tmp);
1301: tmpO = NULL;
1302: }
1303: for (t = 0; t < tmpSize; ++t) {
1304: const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize;
1305: const PetscInt cp = tmp[i];
1306: /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1307: /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1308: const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1309: PetscInt co = tmpO ? tmpO[i] : 0;
1310: PetscInt c;
1312: if (rev) {
1313: PetscInt childSize, coff;
1314: DMPlexGetConeSize(dm, cp, &childSize);
1315: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1316: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1317: }
1318: /* Check for duplicate */
1319: for (c = 0; c < closureSize; c += 2) {
1320: if (closure[c] == cp) break;
1321: }
1322: if (c == closureSize) {
1323: closure[closureSize] = cp;
1324: closure[closureSize+1] = co;
1325: fifo[fifoSize] = cp;
1326: fifo[fifoSize+1] = co;
1327: closureSize += 2;
1328: fifoSize += 2;
1329: }
1330: }
1331: fifoStart += 2;
1332: }
1333: if (numPoints) *numPoints = closureSize/2;
1334: if (points) *points = closure;
1335: DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1336: return(0);
1337: }
1341: /*@C
1342: DMPlexGetTransitiveClosure_Internal - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG with a specified initial orientation
1344: Not collective
1346: Input Parameters:
1347: + mesh - The DMPlex
1348: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1349: . orientation - The orientation of the point
1350: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1351: - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
1353: Output Parameters:
1354: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1355: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
1357: Note:
1358: If using internal storage (points is NULL on input), each call overwrites the last output.
1360: Fortran Notes:
1361: Since it returns an array, this routine is only available in Fortran 90, and you must
1362: include petsc.h90 in your code.
1364: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1366: Level: beginner
1368: .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1369: @*/
1370: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1371: {
1372: DM_Plex *mesh = (DM_Plex*) dm->data;
1373: PetscInt *closure, *fifo;
1374: const PetscInt *tmp = NULL, *tmpO = NULL;
1375: PetscInt tmpSize, t;
1376: PetscInt depth = 0, maxSize;
1377: PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0;
1378: PetscErrorCode ierr;
1382: DMPlexGetDepth(dm, &depth);
1383: /* This is only 1-level */
1384: if (useCone) {
1385: DMPlexGetConeSize(dm, p, &tmpSize);
1386: DMPlexGetCone(dm, p, &tmp);
1387: DMPlexGetConeOrientation(dm, p, &tmpO);
1388: } else {
1389: DMPlexGetSupportSize(dm, p, &tmpSize);
1390: DMPlexGetSupport(dm, p, &tmp);
1391: }
1392: if (depth == 1) {
1393: if (*points) {
1394: closure = *points;
1395: } else {
1396: maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
1397: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1398: }
1399: closure[0] = p; closure[1] = ornt;
1400: for (t = 0; t < tmpSize; ++t, closureSize += 2) {
1401: const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1402: closure[closureSize] = tmp[i];
1403: closure[closureSize+1] = tmpO ? tmpO[i] : 0;
1404: }
1405: if (numPoints) *numPoints = closureSize/2;
1406: if (points) *points = closure;
1407: return(0);
1408: }
1409: maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth+1);
1410: DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);
1411: if (*points) {
1412: closure = *points;
1413: } else {
1414: DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);
1415: }
1416: closure[0] = p; closure[1] = ornt;
1417: for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1418: const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
1419: const PetscInt cp = tmp[i];
1420: PetscInt co = tmpO ? tmpO[i] : 0;
1422: if (ornt < 0) {
1423: PetscInt childSize, coff;
1424: DMPlexGetConeSize(dm, cp, &childSize);
1425: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1426: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1427: }
1428: closure[closureSize] = cp;
1429: closure[closureSize+1] = co;
1430: fifo[fifoSize] = cp;
1431: fifo[fifoSize+1] = co;
1432: }
1433: /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
1434: while (fifoSize - fifoStart) {
1435: const PetscInt q = fifo[fifoStart];
1436: const PetscInt o = fifo[fifoStart+1];
1437: const PetscInt rev = o >= 0 ? 0 : 1;
1438: const PetscInt off = rev ? -(o+1) : o;
1440: if (useCone) {
1441: DMPlexGetConeSize(dm, q, &tmpSize);
1442: DMPlexGetCone(dm, q, &tmp);
1443: DMPlexGetConeOrientation(dm, q, &tmpO);
1444: } else {
1445: DMPlexGetSupportSize(dm, q, &tmpSize);
1446: DMPlexGetSupport(dm, q, &tmp);
1447: tmpO = NULL;
1448: }
1449: for (t = 0; t < tmpSize; ++t) {
1450: const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize;
1451: const PetscInt cp = tmp[i];
1452: /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
1453: /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
1454: const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
1455: PetscInt co = tmpO ? tmpO[i] : 0;
1456: PetscInt c;
1458: if (rev) {
1459: PetscInt childSize, coff;
1460: DMPlexGetConeSize(dm, cp, &childSize);
1461: coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
1462: co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
1463: }
1464: /* Check for duplicate */
1465: for (c = 0; c < closureSize; c += 2) {
1466: if (closure[c] == cp) break;
1467: }
1468: if (c == closureSize) {
1469: closure[closureSize] = cp;
1470: closure[closureSize+1] = co;
1471: fifo[fifoSize] = cp;
1472: fifo[fifoSize+1] = co;
1473: closureSize += 2;
1474: fifoSize += 2;
1475: }
1476: }
1477: fifoStart += 2;
1478: }
1479: if (numPoints) *numPoints = closureSize/2;
1480: if (points) *points = closure;
1481: DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);
1482: return(0);
1483: }
1487: /*@C
1488: DMPlexRestoreTransitiveClosure - Restore the array of points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1490: Not collective
1492: Input Parameters:
1493: + mesh - The DMPlex
1494: . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
1495: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1496: . numPoints - The number of points in the closure, so points[] is of size 2*numPoints, zeroed on exit
1497: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...], zeroed on exit
1499: Note:
1500: If not using internal storage (points is not NULL on input), this call is unnecessary
1502: Fortran Notes:
1503: Since it returns an array, this routine is only available in Fortran 90, and you must
1504: include petsc.h90 in your code.
1506: The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1508: Level: beginner
1510: .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
1511: @*/
1512: PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1513: {
1520: DMRestoreWorkArray(dm, 0, PETSC_INT, points);
1521: if (numPoints) *numPoints = 0;
1522: return(0);
1523: }
1527: /*@
1528: DMPlexGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG
1530: Not collective
1532: Input Parameter:
1533: . mesh - The DMPlex
1535: Output Parameters:
1536: + maxConeSize - The maximum number of in-edges
1537: - maxSupportSize - The maximum number of out-edges
1539: Level: beginner
1541: .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
1542: @*/
1543: PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1544: {
1545: DM_Plex *mesh = (DM_Plex*) dm->data;
1549: if (maxConeSize) *maxConeSize = mesh->maxConeSize;
1550: if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1551: return(0);
1552: }
1556: PetscErrorCode DMSetUp_Plex(DM dm)
1557: {
1558: DM_Plex *mesh = (DM_Plex*) dm->data;
1559: PetscInt size;
1564: PetscSectionSetUp(mesh->coneSection);
1565: PetscSectionGetStorageSize(mesh->coneSection, &size);
1566: PetscMalloc1(size, &mesh->cones);
1567: PetscCalloc1(size, &mesh->coneOrientations);
1568: if (mesh->maxSupportSize) {
1569: PetscSectionSetUp(mesh->supportSection);
1570: PetscSectionGetStorageSize(mesh->supportSection, &size);
1571: PetscMalloc1(size, &mesh->supports);
1572: }
1573: return(0);
1574: }
1578: PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1579: {
1583: if (subdm) {DMClone(dm, subdm);}
1584: DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);
1585: return(0);
1586: }
1590: /*@
1591: DMPlexSymmetrize - Creates support (out-edge) information from cone (in-edge) inoformation
1593: Not collective
1595: Input Parameter:
1596: . mesh - The DMPlex
1598: Output Parameter:
1600: Note:
1601: This should be called after all calls to DMPlexSetCone()
1603: Level: beginner
1605: .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone()
1606: @*/
1607: PetscErrorCode DMPlexSymmetrize(DM dm)
1608: {
1609: DM_Plex *mesh = (DM_Plex*) dm->data;
1610: PetscInt *offsets;
1611: PetscInt supportSize;
1612: PetscInt pStart, pEnd, p;
1617: if (mesh->supports) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Supports were already setup in this DMPlex");
1618: /* Calculate support sizes */
1619: DMPlexGetChart(dm, &pStart, &pEnd);
1620: for (p = pStart; p < pEnd; ++p) {
1621: PetscInt dof, off, c;
1623: PetscSectionGetDof(mesh->coneSection, p, &dof);
1624: PetscSectionGetOffset(mesh->coneSection, p, &off);
1625: for (c = off; c < off+dof; ++c) {
1626: PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1627: }
1628: }
1629: for (p = pStart; p < pEnd; ++p) {
1630: PetscInt dof;
1632: PetscSectionGetDof(mesh->supportSection, p, &dof);
1634: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1635: }
1636: PetscSectionSetUp(mesh->supportSection);
1637: /* Calculate supports */
1638: PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1639: PetscMalloc1(supportSize, &mesh->supports);
1640: PetscCalloc1(pEnd - pStart, &offsets);
1641: for (p = pStart; p < pEnd; ++p) {
1642: PetscInt dof, off, c;
1644: PetscSectionGetDof(mesh->coneSection, p, &dof);
1645: PetscSectionGetOffset(mesh->coneSection, p, &off);
1646: for (c = off; c < off+dof; ++c) {
1647: const PetscInt q = mesh->cones[c];
1648: PetscInt offS;
1650: PetscSectionGetOffset(mesh->supportSection, q, &offS);
1652: mesh->supports[offS+offsets[q]] = p;
1653: ++offsets[q];
1654: }
1655: }
1656: PetscFree(offsets);
1657: return(0);
1658: }
1662: /*@
1663: DMPlexStratify - The Sieve DAG for most topologies is a graded poset (http://en.wikipedia.org/wiki/Graded_poset), and
1664: can be illustrated by Hasse Diagram (a http://en.wikipedia.org/wiki/Hasse_diagram). The strata group all points of the
1665: same grade, and this function calculates the strata. This grade can be seen as the height (or depth) of the point in
1666: the DAG.
1668: Not collective
1670: Input Parameter:
1671: . mesh - The DMPlex
1673: Output Parameter:
1675: Notes:
1676: Concretely, DMPlexStratify() creates a new label named "depth" containing the dimension of each element: 0 for vertices,
1677: 1 for edges, and so on. The depth label can be accessed through DMPlexGetDepthLabel() or DMPlexGetDepthStratum(), or
1678: manually via DMPlexGetLabel(). The height is defined implicitly by height = maxDimension - depth, and can be accessed
1679: via DMPlexGetHeightStratum(). For example, cells have height 0 and faces have height 1.
1681: DMPlexStratify() should be called after all calls to DMPlexSymmetrize()
1683: Level: beginner
1685: .seealso: DMPlexCreate(), DMPlexSymmetrize()
1686: @*/
1687: PetscErrorCode DMPlexStratify(DM dm)
1688: {
1689: DMLabel label;
1690: PetscInt pStart, pEnd, p;
1691: PetscInt numRoots = 0, numLeaves = 0;
1696: PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);
1697: /* Calculate depth */
1698: DMPlexGetChart(dm, &pStart, &pEnd);
1699: DMPlexCreateLabel(dm, "depth");
1700: DMPlexGetDepthLabel(dm, &label);
1701: /* Initialize roots and count leaves */
1702: for (p = pStart; p < pEnd; ++p) {
1703: PetscInt coneSize, supportSize;
1705: DMPlexGetConeSize(dm, p, &coneSize);
1706: DMPlexGetSupportSize(dm, p, &supportSize);
1707: if (!coneSize && supportSize) {
1708: ++numRoots;
1709: DMLabelSetValue(label, p, 0);
1710: } else if (!supportSize && coneSize) {
1711: ++numLeaves;
1712: } else if (!supportSize && !coneSize) {
1713: /* Isolated points */
1714: DMLabelSetValue(label, p, 0);
1715: }
1716: }
1717: if (numRoots + numLeaves == (pEnd - pStart)) {
1718: for (p = pStart; p < pEnd; ++p) {
1719: PetscInt coneSize, supportSize;
1721: DMPlexGetConeSize(dm, p, &coneSize);
1722: DMPlexGetSupportSize(dm, p, &supportSize);
1723: if (!supportSize && coneSize) {
1724: DMLabelSetValue(label, p, 1);
1725: }
1726: }
1727: } else {
1728: IS pointIS;
1729: PetscInt numPoints = 0, level = 0;
1731: DMLabelGetStratumIS(label, level, &pointIS);
1732: if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1733: while (numPoints) {
1734: const PetscInt *points;
1735: const PetscInt newLevel = level+1;
1737: ISGetIndices(pointIS, &points);
1738: for (p = 0; p < numPoints; ++p) {
1739: const PetscInt point = points[p];
1740: const PetscInt *support;
1741: PetscInt supportSize, s;
1743: DMPlexGetSupportSize(dm, point, &supportSize);
1744: DMPlexGetSupport(dm, point, &support);
1745: for (s = 0; s < supportSize; ++s) {
1746: DMLabelSetValue(label, support[s], newLevel);
1747: }
1748: }
1749: ISRestoreIndices(pointIS, &points);
1750: ++level;
1751: ISDestroy(&pointIS);
1752: DMLabelGetStratumIS(label, level, &pointIS);
1753: if (pointIS) {ISGetLocalSize(pointIS, &numPoints);}
1754: else {numPoints = 0;}
1755: }
1756: ISDestroy(&pointIS);
1757: }
1758: PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);
1759: return(0);
1760: }
1764: /*@C
1765: DMPlexGetJoin - Get an array for the join of the set of points
1767: Not Collective
1769: Input Parameters:
1770: + dm - The DMPlex object
1771: . numPoints - The number of input points for the join
1772: - points - The input points
1774: Output Parameters:
1775: + numCoveredPoints - The number of points in the join
1776: - coveredPoints - The points in the join
1778: Level: intermediate
1780: Note: Currently, this is restricted to a single level join
1782: Fortran Notes:
1783: Since it returns an array, this routine is only available in Fortran 90, and you must
1784: include petsc.h90 in your code.
1786: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1788: .keywords: mesh
1789: .seealso: DMPlexRestoreJoin(), DMPlexGetMeet()
1790: @*/
1791: PetscErrorCode DMPlexGetJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1792: {
1793: DM_Plex *mesh = (DM_Plex*) dm->data;
1794: PetscInt *join[2];
1795: PetscInt joinSize, i = 0;
1796: PetscInt dof, off, p, c, m;
1804: DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[0]);
1805: DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1]);
1806: /* Copy in support of first point */
1807: PetscSectionGetDof(mesh->supportSection, points[0], &dof);
1808: PetscSectionGetOffset(mesh->supportSection, points[0], &off);
1809: for (joinSize = 0; joinSize < dof; ++joinSize) {
1810: join[i][joinSize] = mesh->supports[off+joinSize];
1811: }
1812: /* Check each successive support */
1813: for (p = 1; p < numPoints; ++p) {
1814: PetscInt newJoinSize = 0;
1816: PetscSectionGetDof(mesh->supportSection, points[p], &dof);
1817: PetscSectionGetOffset(mesh->supportSection, points[p], &off);
1818: for (c = 0; c < dof; ++c) {
1819: const PetscInt point = mesh->supports[off+c];
1821: for (m = 0; m < joinSize; ++m) {
1822: if (point == join[i][m]) {
1823: join[1-i][newJoinSize++] = point;
1824: break;
1825: }
1826: }
1827: }
1828: joinSize = newJoinSize;
1829: i = 1-i;
1830: }
1831: *numCoveredPoints = joinSize;
1832: *coveredPoints = join[i];
1833: DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
1834: return(0);
1835: }
1839: /*@C
1840: DMPlexRestoreJoin - Restore an array for the join of the set of points
1842: Not Collective
1844: Input Parameters:
1845: + dm - The DMPlex object
1846: . numPoints - The number of input points for the join
1847: - points - The input points
1849: Output Parameters:
1850: + numCoveredPoints - The number of points in the join
1851: - coveredPoints - The points in the join
1853: Fortran Notes:
1854: Since it returns an array, this routine is only available in Fortran 90, and you must
1855: include petsc.h90 in your code.
1857: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1859: Level: intermediate
1861: .keywords: mesh
1862: .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet()
1863: @*/
1864: PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1865: {
1873: DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
1874: if (numCoveredPoints) *numCoveredPoints = 0;
1875: return(0);
1876: }
1880: /*@C
1881: DMPlexGetFullJoin - Get an array for the join of the set of points
1883: Not Collective
1885: Input Parameters:
1886: + dm - The DMPlex object
1887: . numPoints - The number of input points for the join
1888: - points - The input points
1890: Output Parameters:
1891: + numCoveredPoints - The number of points in the join
1892: - coveredPoints - The points in the join
1894: Fortran Notes:
1895: Since it returns an array, this routine is only available in Fortran 90, and you must
1896: include petsc.h90 in your code.
1898: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
1900: Level: intermediate
1902: .keywords: mesh
1903: .seealso: DMPlexGetJoin(), DMPlexRestoreJoin(), DMPlexGetMeet()
1904: @*/
1905: PetscErrorCode DMPlexGetFullJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
1906: {
1907: DM_Plex *mesh = (DM_Plex*) dm->data;
1908: PetscInt *offsets, **closures;
1909: PetscInt *join[2];
1910: PetscInt depth = 0, maxSize, joinSize = 0, i = 0;
1911: PetscInt p, d, c, m;
1920: DMPlexGetDepth(dm, &depth);
1921: PetscCalloc1(numPoints, &closures);
1922: DMGetWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
1923: maxSize = PetscPowInt(mesh->maxSupportSize,depth+1);
1924: DMGetWorkArray(dm, maxSize, PETSC_INT, &join[0]);
1925: DMGetWorkArray(dm, maxSize, PETSC_INT, &join[1]);
1927: for (p = 0; p < numPoints; ++p) {
1928: PetscInt closureSize;
1930: DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closures[p]);
1932: offsets[p*(depth+2)+0] = 0;
1933: for (d = 0; d < depth+1; ++d) {
1934: PetscInt pStart, pEnd, i;
1936: DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
1937: for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) {
1938: if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
1939: offsets[p*(depth+2)+d+1] = i;
1940: break;
1941: }
1942: }
1943: if (i == closureSize) offsets[p*(depth+2)+d+1] = i;
1944: }
1945: if (offsets[p*(depth+2)+depth+1] != closureSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Total size of closure %D should be %D", offsets[p*(depth+2)+depth+1], closureSize);
1946: }
1947: for (d = 0; d < depth+1; ++d) {
1948: PetscInt dof;
1950: /* Copy in support of first point */
1951: dof = offsets[d+1] - offsets[d];
1952: for (joinSize = 0; joinSize < dof; ++joinSize) {
1953: join[i][joinSize] = closures[0][(offsets[d]+joinSize)*2];
1954: }
1955: /* Check each successive cone */
1956: for (p = 1; p < numPoints && joinSize; ++p) {
1957: PetscInt newJoinSize = 0;
1959: dof = offsets[p*(depth+2)+d+1] - offsets[p*(depth+2)+d];
1960: for (c = 0; c < dof; ++c) {
1961: const PetscInt point = closures[p][(offsets[p*(depth+2)+d]+c)*2];
1963: for (m = 0; m < joinSize; ++m) {
1964: if (point == join[i][m]) {
1965: join[1-i][newJoinSize++] = point;
1966: break;
1967: }
1968: }
1969: }
1970: joinSize = newJoinSize;
1971: i = 1-i;
1972: }
1973: if (joinSize) break;
1974: }
1975: *numCoveredPoints = joinSize;
1976: *coveredPoints = join[i];
1977: for (p = 0; p < numPoints; ++p) {
1978: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, NULL, &closures[p]);
1979: }
1980: PetscFree(closures);
1981: DMRestoreWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);
1982: DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);
1983: return(0);
1984: }
1988: /*@C
1989: DMPlexGetMeet - Get an array for the meet of the set of points
1991: Not Collective
1993: Input Parameters:
1994: + dm - The DMPlex object
1995: . numPoints - The number of input points for the meet
1996: - points - The input points
1998: Output Parameters:
1999: + numCoveredPoints - The number of points in the meet
2000: - coveredPoints - The points in the meet
2002: Level: intermediate
2004: Note: Currently, this is restricted to a single level meet
2006: Fortran Notes:
2007: Since it returns an array, this routine is only available in Fortran 90, and you must
2008: include petsc.h90 in your code.
2010: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2012: .keywords: mesh
2013: .seealso: DMPlexRestoreMeet(), DMPlexGetJoin()
2014: @*/
2015: PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2016: {
2017: DM_Plex *mesh = (DM_Plex*) dm->data;
2018: PetscInt *meet[2];
2019: PetscInt meetSize, i = 0;
2020: PetscInt dof, off, p, c, m;
2028: DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[0]);
2029: DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1]);
2030: /* Copy in cone of first point */
2031: PetscSectionGetDof(mesh->coneSection, points[0], &dof);
2032: PetscSectionGetOffset(mesh->coneSection, points[0], &off);
2033: for (meetSize = 0; meetSize < dof; ++meetSize) {
2034: meet[i][meetSize] = mesh->cones[off+meetSize];
2035: }
2036: /* Check each successive cone */
2037: for (p = 1; p < numPoints; ++p) {
2038: PetscInt newMeetSize = 0;
2040: PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2041: PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2042: for (c = 0; c < dof; ++c) {
2043: const PetscInt point = mesh->cones[off+c];
2045: for (m = 0; m < meetSize; ++m) {
2046: if (point == meet[i][m]) {
2047: meet[1-i][newMeetSize++] = point;
2048: break;
2049: }
2050: }
2051: }
2052: meetSize = newMeetSize;
2053: i = 1-i;
2054: }
2055: *numCoveringPoints = meetSize;
2056: *coveringPoints = meet[i];
2057: DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2058: return(0);
2059: }
2063: /*@C
2064: DMPlexRestoreMeet - Restore an array for the meet of the set of points
2066: Not Collective
2068: Input Parameters:
2069: + dm - The DMPlex object
2070: . numPoints - The number of input points for the meet
2071: - points - The input points
2073: Output Parameters:
2074: + numCoveredPoints - The number of points in the meet
2075: - coveredPoints - The points in the meet
2077: Level: intermediate
2079: Fortran Notes:
2080: Since it returns an array, this routine is only available in Fortran 90, and you must
2081: include petsc.h90 in your code.
2083: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2085: .keywords: mesh
2086: .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin()
2087: @*/
2088: PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2089: {
2097: DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);
2098: if (numCoveredPoints) *numCoveredPoints = 0;
2099: return(0);
2100: }
2104: /*@C
2105: DMPlexGetFullMeet - Get an array for the meet of the set of points
2107: Not Collective
2109: Input Parameters:
2110: + dm - The DMPlex object
2111: . numPoints - The number of input points for the meet
2112: - points - The input points
2114: Output Parameters:
2115: + numCoveredPoints - The number of points in the meet
2116: - coveredPoints - The points in the meet
2118: Level: intermediate
2120: Fortran Notes:
2121: Since it returns an array, this routine is only available in Fortran 90, and you must
2122: include petsc.h90 in your code.
2124: The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array.
2126: .keywords: mesh
2127: .seealso: DMPlexGetMeet(), DMPlexRestoreMeet(), DMPlexGetJoin()
2128: @*/
2129: PetscErrorCode DMPlexGetFullMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2130: {
2131: DM_Plex *mesh = (DM_Plex*) dm->data;
2132: PetscInt *offsets, **closures;
2133: PetscInt *meet[2];
2134: PetscInt height = 0, maxSize, meetSize = 0, i = 0;
2135: PetscInt p, h, c, m;
2144: DMPlexGetDepth(dm, &height);
2145: PetscMalloc1(numPoints, &closures);
2146: DMGetWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2147: maxSize = PetscPowInt(mesh->maxConeSize,height+1);
2148: DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[0]);
2149: DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[1]);
2151: for (p = 0; p < numPoints; ++p) {
2152: PetscInt closureSize;
2154: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closures[p]);
2156: offsets[p*(height+2)+0] = 0;
2157: for (h = 0; h < height+1; ++h) {
2158: PetscInt pStart, pEnd, i;
2160: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
2161: for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) {
2162: if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) {
2163: offsets[p*(height+2)+h+1] = i;
2164: break;
2165: }
2166: }
2167: if (i == closureSize) offsets[p*(height+2)+h+1] = i;
2168: }
2169: if (offsets[p*(height+2)+height+1] != closureSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Total size of closure %D should be %D", offsets[p*(height+2)+height+1], closureSize);
2170: }
2171: for (h = 0; h < height+1; ++h) {
2172: PetscInt dof;
2174: /* Copy in cone of first point */
2175: dof = offsets[h+1] - offsets[h];
2176: for (meetSize = 0; meetSize < dof; ++meetSize) {
2177: meet[i][meetSize] = closures[0][(offsets[h]+meetSize)*2];
2178: }
2179: /* Check each successive cone */
2180: for (p = 1; p < numPoints && meetSize; ++p) {
2181: PetscInt newMeetSize = 0;
2183: dof = offsets[p*(height+2)+h+1] - offsets[p*(height+2)+h];
2184: for (c = 0; c < dof; ++c) {
2185: const PetscInt point = closures[p][(offsets[p*(height+2)+h]+c)*2];
2187: for (m = 0; m < meetSize; ++m) {
2188: if (point == meet[i][m]) {
2189: meet[1-i][newMeetSize++] = point;
2190: break;
2191: }
2192: }
2193: }
2194: meetSize = newMeetSize;
2195: i = 1-i;
2196: }
2197: if (meetSize) break;
2198: }
2199: *numCoveredPoints = meetSize;
2200: *coveredPoints = meet[i];
2201: for (p = 0; p < numPoints; ++p) {
2202: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, NULL, &closures[p]);
2203: }
2204: PetscFree(closures);
2205: DMRestoreWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);
2206: DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);
2207: return(0);
2208: }
2212: /*@C
2213: DMPlexEqual - Determine if two DMs have the same topology
2215: Not Collective
2217: Input Parameters:
2218: + dmA - A DMPlex object
2219: - dmB - A DMPlex object
2221: Output Parameters:
2222: . equal - PETSC_TRUE if the topologies are identical
2224: Level: intermediate
2226: Notes:
2227: We are not solving graph isomorphism, so we do not permutation.
2229: .keywords: mesh
2230: .seealso: DMPlexGetCone()
2231: @*/
2232: PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal)
2233: {
2234: PetscInt depth, depthB, pStart, pEnd, pStartB, pEndB, p;
2242: *equal = PETSC_FALSE;
2243: DMPlexGetDepth(dmA, &depth);
2244: DMPlexGetDepth(dmB, &depthB);
2245: if (depth != depthB) return(0);
2246: DMPlexGetChart(dmA, &pStart, &pEnd);
2247: DMPlexGetChart(dmB, &pStartB, &pEndB);
2248: if ((pStart != pStartB) || (pEnd != pEndB)) return(0);
2249: for (p = pStart; p < pEnd; ++p) {
2250: const PetscInt *cone, *coneB, *ornt, *orntB, *support, *supportB;
2251: PetscInt coneSize, coneSizeB, c, supportSize, supportSizeB, s;
2253: DMPlexGetConeSize(dmA, p, &coneSize);
2254: DMPlexGetCone(dmA, p, &cone);
2255: DMPlexGetConeOrientation(dmA, p, &ornt);
2256: DMPlexGetConeSize(dmB, p, &coneSizeB);
2257: DMPlexGetCone(dmB, p, &coneB);
2258: DMPlexGetConeOrientation(dmB, p, &orntB);
2259: if (coneSize != coneSizeB) return(0);
2260: for (c = 0; c < coneSize; ++c) {
2261: if (cone[c] != coneB[c]) return(0);
2262: if (ornt[c] != orntB[c]) return(0);
2263: }
2264: DMPlexGetSupportSize(dmA, p, &supportSize);
2265: DMPlexGetSupport(dmA, p, &support);
2266: DMPlexGetSupportSize(dmB, p, &supportSizeB);
2267: DMPlexGetSupport(dmB, p, &supportB);
2268: if (supportSize != supportSizeB) return(0);
2269: for (s = 0; s < supportSize; ++s) {
2270: if (support[s] != supportB[s]) return(0);
2271: }
2272: }
2273: *equal = PETSC_TRUE;
2274: return(0);
2275: }
2279: PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices)
2280: {
2281: MPI_Comm comm;
2285: PetscObjectGetComm((PetscObject)dm,&comm);
2287: switch (cellDim) {
2288: case 0:
2289: *numFaceVertices = 0;
2290: break;
2291: case 1:
2292: *numFaceVertices = 1;
2293: break;
2294: case 2:
2295: switch (numCorners) {
2296: case 3: /* triangle */
2297: *numFaceVertices = 2; /* Edge has 2 vertices */
2298: break;
2299: case 4: /* quadrilateral */
2300: *numFaceVertices = 2; /* Edge has 2 vertices */
2301: break;
2302: case 6: /* quadratic triangle, tri and quad cohesive Lagrange cells */
2303: *numFaceVertices = 3; /* Edge has 3 vertices */
2304: break;
2305: case 9: /* quadratic quadrilateral, quadratic quad cohesive Lagrange cells */
2306: *numFaceVertices = 3; /* Edge has 3 vertices */
2307: break;
2308: default:
2309: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim);
2310: }
2311: break;
2312: case 3:
2313: switch (numCorners) {
2314: case 4: /* tetradehdron */
2315: *numFaceVertices = 3; /* Face has 3 vertices */
2316: break;
2317: case 6: /* tet cohesive cells */
2318: *numFaceVertices = 4; /* Face has 4 vertices */
2319: break;
2320: case 8: /* hexahedron */
2321: *numFaceVertices = 4; /* Face has 4 vertices */
2322: break;
2323: case 9: /* tet cohesive Lagrange cells */
2324: *numFaceVertices = 6; /* Face has 6 vertices */
2325: break;
2326: case 10: /* quadratic tetrahedron */
2327: *numFaceVertices = 6; /* Face has 6 vertices */
2328: break;
2329: case 12: /* hex cohesive Lagrange cells */
2330: *numFaceVertices = 6; /* Face has 6 vertices */
2331: break;
2332: case 18: /* quadratic tet cohesive Lagrange cells */
2333: *numFaceVertices = 6; /* Face has 6 vertices */
2334: break;
2335: case 27: /* quadratic hexahedron, quadratic hex cohesive Lagrange cells */
2336: *numFaceVertices = 9; /* Face has 9 vertices */
2337: break;
2338: default:
2339: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim);
2340: }
2341: break;
2342: default:
2343: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid cell dimension %d", cellDim);
2344: }
2345: return(0);
2346: }
2350: /* Trys to give the mesh a consistent orientation */
2351: PetscErrorCode DMPlexOrient(DM dm)
2352: {
2353: PetscBT seenCells, flippedCells, seenFaces;
2354: PetscInt *faceFIFO, fTop, fBottom;
2355: PetscInt dim, h, cStart, cEnd, c, fStart, fEnd, face, maxConeSize, *revcone, *revconeO;
2359: /* Truth Table
2360: mismatch flips do action mismatch flipA ^ flipB action
2361: F 0 flips no F F F
2362: F 1 flip yes F T T
2363: F 2 flips no T F T
2364: T 0 flips yes T T F
2365: T 1 flip no
2366: T 2 flips yes
2367: */
2368: DMPlexGetDimension(dm, &dim);
2369: DMPlexGetVTKCellHeight(dm, &h);
2370: DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
2371: DMPlexGetHeightStratum(dm, h+1, &fStart, &fEnd);
2372: PetscBTCreate(cEnd - cStart, &seenCells);
2373: PetscBTMemzero(cEnd - cStart, seenCells);
2374: PetscBTCreate(cEnd - cStart, &flippedCells);
2375: PetscBTMemzero(cEnd - cStart, flippedCells);
2376: PetscBTCreate(fEnd - fStart, &seenFaces);
2377: PetscBTMemzero(fEnd - fStart, seenFaces);
2378: PetscMalloc1((fEnd - fStart), &faceFIFO);
2379: fTop = fBottom = 0;
2380: /* Initialize FIFO with first cell */
2381: if (cEnd > cStart) {
2382: const PetscInt *cone;
2383: PetscInt coneSize;
2385: DMPlexGetConeSize(dm, cStart, &coneSize);
2386: DMPlexGetCone(dm, cStart, &cone);
2387: for (c = 0; c < coneSize; ++c) {
2388: faceFIFO[fBottom++] = cone[c];
2389: PetscBTSet(seenFaces, cone[c]-fStart);
2390: }
2391: }
2392: /* Consider each face in FIFO */
2393: while (fTop < fBottom) {
2394: const PetscInt *support, *coneA, *coneB, *coneOA, *coneOB;
2395: PetscInt supportSize, coneSizeA, coneSizeB, posA = -1, posB = -1;
2396: PetscInt seenA, flippedA, seenB, flippedB, mismatch;
2398: face = faceFIFO[fTop++];
2399: DMPlexGetSupportSize(dm, face, &supportSize);
2400: DMPlexGetSupport(dm, face, &support);
2401: if (supportSize < 2) continue;
2402: if (supportSize != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %d", supportSize);
2403: seenA = PetscBTLookup(seenCells, support[0]-cStart);
2404: flippedA = PetscBTLookup(flippedCells, support[0]-cStart) ? 1 : 0;
2405: seenB = PetscBTLookup(seenCells, support[1]-cStart);
2406: flippedB = PetscBTLookup(flippedCells, support[1]-cStart) ? 1 : 0;
2408: DMPlexGetConeSize(dm, support[0], &coneSizeA);
2409: DMPlexGetConeSize(dm, support[1], &coneSizeB);
2410: DMPlexGetCone(dm, support[0], &coneA);
2411: DMPlexGetCone(dm, support[1], &coneB);
2412: DMPlexGetConeOrientation(dm, support[0], &coneOA);
2413: DMPlexGetConeOrientation(dm, support[1], &coneOB);
2414: for (c = 0; c < coneSizeA; ++c) {
2415: if (!PetscBTLookup(seenFaces, coneA[c]-fStart)) {
2416: faceFIFO[fBottom++] = coneA[c];
2417: PetscBTSet(seenFaces, coneA[c]-fStart);
2418: }
2419: if (coneA[c] == face) posA = c;
2420: if (fBottom > fEnd-fStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %d was pushed exceeding capacity %d > %d", coneA[c], fBottom, fEnd-fStart);
2421: }
2422: if (posA < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d could not be located in cell %d", face, support[0]);
2423: for (c = 0; c < coneSizeB; ++c) {
2424: if (!PetscBTLookup(seenFaces, coneB[c]-fStart)) {
2425: faceFIFO[fBottom++] = coneB[c];
2426: PetscBTSet(seenFaces, coneB[c]-fStart);
2427: }
2428: if (coneB[c] == face) posB = c;
2429: if (fBottom > fEnd-fStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %d was pushed exceeding capacity %d > %d", coneA[c], fBottom, fEnd-fStart);
2430: }
2431: if (posB < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d could not be located in cell %d", face, support[1]);
2433: if (dim == 1) {
2434: mismatch = posA == posB;
2435: } else {
2436: mismatch = coneOA[posA] == coneOB[posB];
2437: }
2439: if (mismatch ^ (flippedA ^ flippedB)) {
2440: if (seenA && seenB) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen cells %d and %d do not match: Fault mesh is non-orientable", support[0], support[1]);
2441: if (!seenA && !flippedA) {
2442: PetscBTSet(flippedCells, support[0]-cStart);
2443: } else if (!seenB && !flippedB) {
2444: PetscBTSet(flippedCells, support[1]-cStart);
2445: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
2446: } else if (mismatch && flippedA && flippedB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
2447: PetscBTSet(seenCells, support[0]-cStart);
2448: PetscBTSet(seenCells, support[1]-cStart);
2449: }
2450: /* Now all subdomains are oriented, but we need a consistent parallel orientation */
2451: {
2452: /* Find a representative face (edge) separating pairs of procs */
2453: PetscSF sf;
2454: const PetscInt *lpoints;
2455: const PetscSFNode *rpoints;
2456: PetscInt *neighbors, *nranks;
2457: PetscInt numLeaves, numRoots, numNeighbors = 0, l, n;
2459: DMGetPointSF(dm, &sf);
2460: PetscSFGetGraph(sf, &numRoots, &numLeaves, &lpoints, &rpoints);
2461: if (numLeaves >= 0) {
2462: const PetscInt *cone, *ornt, *support;
2463: PetscInt coneSize, supportSize;
2464: int *rornt, *lornt; /* PetscSF cannot handle smaller than int */
2465: PetscBool *match, flipped = PETSC_FALSE;
2467: PetscMalloc1(numLeaves,&neighbors);
2468: /* I know this is p^2 time in general, but for bounded degree its alright */
2469: for (l = 0; l < numLeaves; ++l) {
2470: const PetscInt face = lpoints[l];
2471: if ((face >= fStart) && (face < fEnd)) {
2472: const PetscInt rank = rpoints[l].rank;
2473: for (n = 0; n < numNeighbors; ++n) if (rank == rpoints[neighbors[n]].rank) break;
2474: if (n >= numNeighbors) {
2475: PetscInt supportSize;
2476: DMPlexGetSupportSize(dm, face, &supportSize);
2477: if (supportSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Boundary faces should see one cell, not %d", supportSize);
2478: neighbors[numNeighbors++] = l;
2479: }
2480: }
2481: }
2482: PetscCalloc4(numNeighbors,&match,numNeighbors,&nranks,numRoots,&rornt,numRoots,&lornt);
2483: for (face = fStart; face < fEnd; ++face) {
2484: DMPlexGetSupportSize(dm, face, &supportSize);
2485: if (supportSize != 1) continue;
2486: DMPlexGetSupport(dm, face, &support);
2488: DMPlexGetCone(dm, support[0], &cone);
2489: DMPlexGetConeSize(dm, support[0], &coneSize);
2490: DMPlexGetConeOrientation(dm, support[0], &ornt);
2491: for (c = 0; c < coneSize; ++c) if (cone[c] == face) break;
2492: if (dim == 1) {
2493: /* Use cone position instead, shifted to -1 or 1 */
2494: rornt[face] = c*2-1;
2495: } else {
2496: if (PetscBTLookup(flippedCells, support[0]-cStart)) rornt[face] = ornt[c] < 0 ? -1 : 1;
2497: else rornt[face] = ornt[c] < 0 ? 1 : -1;
2498: }
2499: }
2500: /* Mark each edge with match or nomatch */
2501: PetscSFBcastBegin(sf, MPI_INT, rornt, lornt);
2502: PetscSFBcastEnd(sf, MPI_INT, rornt, lornt);
2503: for (n = 0; n < numNeighbors; ++n) {
2504: const PetscInt face = lpoints[neighbors[n]];
2506: if (rornt[face]*lornt[face] < 0) match[n] = PETSC_TRUE;
2507: else match[n] = PETSC_FALSE;
2508: nranks[n] = rpoints[neighbors[n]].rank;
2509: }
2510: /* Collect the graph on 0 */
2511: {
2512: MPI_Comm comm = PetscObjectComm((PetscObject) sf);
2513: Mat G;
2514: PetscBT seenProcs, flippedProcs;
2515: PetscInt *procFIFO, pTop, pBottom;
2516: PetscInt *adj = NULL;
2517: PetscBool *val = NULL;
2518: PetscMPIInt *recvcounts = NULL, *displs = NULL, p;
2519: PetscMPIInt N = numNeighbors, numProcs = 0, rank;
2520: PetscInt debug = 0;
2522: MPI_Comm_rank(comm, &rank);
2523: if (!rank) {MPI_Comm_size(comm, &numProcs);}
2524: PetscCalloc2(numProcs,&recvcounts,numProcs+1,&displs);
2525: MPI_Gather(&N, 1, MPI_INT, recvcounts, 1, MPI_INT, 0, comm);
2526: for (p = 0; p < numProcs; ++p) {
2527: displs[p+1] = displs[p] + recvcounts[p];
2528: }
2529: if (!rank) {PetscMalloc2(displs[numProcs],&adj,displs[numProcs],&val);}
2530: MPI_Gatherv(nranks, numNeighbors, MPIU_INT, adj, recvcounts, displs, MPIU_INT, 0, comm);
2531: MPI_Gatherv(match, numNeighbors, MPIU_BOOL, val, recvcounts, displs, MPIU_BOOL, 0, comm);
2532: if (debug) {
2533: for (p = 0; p < numProcs; ++p) {
2534: PetscPrintf(comm, "Proc %d:\n", p);
2535: for (n = 0; n < recvcounts[p]; ++n) {
2536: PetscPrintf(comm, " edge %d (%d):\n", adj[displs[p]+n], val[displs[p]+n]);
2537: }
2538: }
2539: }
2540: /* Symmetrize the graph */
2541: MatCreate(PETSC_COMM_SELF, &G);
2542: MatSetSizes(G, numProcs, numProcs, numProcs, numProcs);
2543: MatSetUp(G);
2544: for (p = 0; p < numProcs; ++p) {
2545: for (n = 0; n < recvcounts[p]; ++n) {
2546: const PetscInt r = p;
2547: const PetscInt q = adj[displs[p]+n];
2548: const PetscScalar o = val[displs[p]+n] ? 1.0 : 0.0;
2550: MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES);
2551: MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES);
2552: }
2553: }
2554: MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY);
2555: MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY);
2557: PetscBTCreate(numProcs, &seenProcs);
2558: PetscBTMemzero(numProcs, seenProcs);
2559: PetscBTCreate(numProcs, &flippedProcs);
2560: PetscBTMemzero(numProcs, flippedProcs);
2561: PetscMalloc1(numProcs,&procFIFO);
2562: pTop = pBottom = 0;
2563: for (p = 0; p < numProcs; ++p) {
2564: if (PetscBTLookup(seenProcs, p)) continue;
2565: /* Initialize FIFO with next proc */
2566: procFIFO[pBottom++] = p;
2567: PetscBTSet(seenProcs, p);
2568: /* Consider each proc in FIFO */
2569: while (pTop < pBottom) {
2570: const PetscScalar *ornt;
2571: const PetscInt *neighbors;
2572: PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors;
2574: proc = procFIFO[pTop++];
2575: flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0;
2576: MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt);
2577: /* Loop over neighboring procs */
2578: for (n = 0; n < numNeighbors; ++n) {
2579: nproc = neighbors[n];
2580: mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1;
2581: seen = PetscBTLookup(seenProcs, nproc);
2582: flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0;
2584: if (mismatch ^ (flippedA ^ flippedB)) {
2585: if (seen) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen procs %d and %d do not match: Fault mesh is non-orientable", proc, nproc);
2586: if (!flippedB) {
2587: PetscBTSet(flippedProcs, nproc);
2588: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable");
2589: } else if (mismatch && flippedA && flippedB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable");
2590: if (!seen) {
2591: procFIFO[pBottom++] = nproc;
2592: PetscBTSet(seenProcs, nproc);
2593: }
2594: }
2595: }
2596: }
2597: PetscFree(procFIFO);
2598: MatDestroy(&G);
2600: PetscFree2(recvcounts,displs);
2601: PetscFree2(adj,val);
2602: {
2603: PetscBool *flips;
2605: PetscMalloc1(numProcs,&flips);
2606: for (p = 0; p < numProcs; ++p) {
2607: flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE;
2608: if (debug && flips[p]) {PetscPrintf(comm, "Flipping Proc %d:\n", p);}
2609: }
2610: MPI_Scatter(flips, 1, MPIU_BOOL, &flipped, 1, MPIU_BOOL, 0, comm);
2611: PetscFree(flips);
2612: }
2613: PetscBTDestroy(&seenProcs);
2614: PetscBTDestroy(&flippedProcs);
2615: }
2616: PetscFree4(match,nranks,rornt,lornt);
2617: PetscFree(neighbors);
2618: if (flipped) {for (c = cStart; c < cEnd; ++c) {PetscBTNegate(flippedCells, c-cStart);}}
2619: }
2620: }
2621: /* Reverse flipped cells in the mesh */
2622: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
2623: DMGetWorkArray(dm, maxConeSize, PETSC_INT, &revcone);
2624: DMGetWorkArray(dm, maxConeSize, PETSC_INT, &revconeO);
2625: for (c = cStart; c < cEnd; ++c) {
2626: const PetscInt *cone, *coneO, *support;
2627: PetscInt coneSize, supportSize, faceSize, cp, sp;
2629: if (!PetscBTLookup(flippedCells, c-cStart)) continue;
2630: DMPlexGetConeSize(dm, c, &coneSize);
2631: DMPlexGetCone(dm, c, &cone);
2632: DMPlexGetConeOrientation(dm, c, &coneO);
2633: for (cp = 0; cp < coneSize; ++cp) {
2634: const PetscInt rcp = coneSize-cp-1;
2636: DMPlexGetConeSize(dm, cone[rcp], &faceSize);
2637: revcone[cp] = cone[rcp];
2638: revconeO[cp] = coneO[rcp] >= 0 ? -(faceSize-coneO[rcp]) : faceSize+coneO[rcp];
2639: }
2640: DMPlexSetCone(dm, c, revcone);
2641: DMPlexSetConeOrientation(dm, c, revconeO);
2642: /* Reverse orientations of support */
2643: faceSize = coneSize;
2644: DMPlexGetSupportSize(dm, c, &supportSize);
2645: DMPlexGetSupport(dm, c, &support);
2646: for (sp = 0; sp < supportSize; ++sp) {
2647: DMPlexGetConeSize(dm, support[sp], &coneSize);
2648: DMPlexGetCone(dm, support[sp], &cone);
2649: DMPlexGetConeOrientation(dm, support[sp], &coneO);
2650: for (cp = 0; cp < coneSize; ++cp) {
2651: if (cone[cp] != c) continue;
2652: DMPlexInsertConeOrientation(dm, support[sp], cp, coneO[cp] >= 0 ? -(faceSize-coneO[cp]) : faceSize+coneO[cp]);
2653: }
2654: }
2655: }
2656: DMRestoreWorkArray(dm, maxConeSize, PETSC_INT, &revcone);
2657: DMRestoreWorkArray(dm, maxConeSize, PETSC_INT, &revconeO);
2658: PetscBTDestroy(&seenCells);
2659: PetscBTDestroy(&flippedCells);
2660: PetscBTDestroy(&seenFaces);
2661: PetscFree(faceFIFO);
2662: return(0);
2663: }
2667: PetscErrorCode DMPlexInvertCell_Internal(PetscInt dim, PetscInt numCorners, PetscInt cone[])
2668: {
2669: int tmpc;
2672: if (dim != 3) return(0);
2673: switch (numCorners) {
2674: case 4:
2675: tmpc = cone[0];
2676: cone[0] = cone[1];
2677: cone[1] = tmpc;
2678: break;
2679: case 8:
2680: tmpc = cone[1];
2681: cone[1] = cone[3];
2682: cone[3] = tmpc;
2683: break;
2684: default: break;
2685: }
2686: return(0);
2687: }
2691: /*@C
2692: DMPlexInvertCell - This flips tetrahedron and hexahedron orientation since Plex stores them internally with outward normals. Other cells are left untouched.
2694: Input Parameters:
2695: + numCorners - The number of vertices in a cell
2696: - cone - The incoming cone
2698: Output Parameter:
2699: . cone - The inverted cone (in-place)
2701: Level: developer
2703: .seealso: DMPlexGenerate()
2704: @*/
2705: PetscErrorCode DMPlexInvertCell(PetscInt dim, PetscInt numCorners, int cone[])
2706: {
2707: int tmpc;
2710: if (dim != 3) return(0);
2711: switch (numCorners) {
2712: case 4:
2713: tmpc = cone[0];
2714: cone[0] = cone[1];
2715: cone[1] = tmpc;
2716: break;
2717: case 8:
2718: tmpc = cone[1];
2719: cone[1] = cone[3];
2720: cone[3] = tmpc;
2721: break;
2722: default: break;
2723: }
2724: return(0);
2725: }
2729: /* This is to fix the tetrahedron orientation from TetGen */
2730: PETSC_UNUSED static PetscErrorCode DMPlexInvertCells_Internal(PetscInt dim, PetscInt numCells, PetscInt numCorners, int cells[])
2731: {
2732: PetscInt bound = numCells*numCorners, coff;
2736: for (coff = 0; coff < bound; coff += numCorners) {
2737: DMPlexInvertCell(dim, numCorners, &cells[coff]);
2738: }
2739: return(0);
2740: }
2742: #if defined(PETSC_HAVE_TRIANGLE)
2743: #include <triangle.h>
2747: PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx)
2748: {
2750: inputCtx->numberofpoints = 0;
2751: inputCtx->numberofpointattributes = 0;
2752: inputCtx->pointlist = NULL;
2753: inputCtx->pointattributelist = NULL;
2754: inputCtx->pointmarkerlist = NULL;
2755: inputCtx->numberofsegments = 0;
2756: inputCtx->segmentlist = NULL;
2757: inputCtx->segmentmarkerlist = NULL;
2758: inputCtx->numberoftriangleattributes = 0;
2759: inputCtx->trianglelist = NULL;
2760: inputCtx->numberofholes = 0;
2761: inputCtx->holelist = NULL;
2762: inputCtx->numberofregions = 0;
2763: inputCtx->regionlist = NULL;
2764: return(0);
2765: }
2769: PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
2770: {
2772: outputCtx->numberofpoints = 0;
2773: outputCtx->pointlist = NULL;
2774: outputCtx->pointattributelist = NULL;
2775: outputCtx->pointmarkerlist = NULL;
2776: outputCtx->numberoftriangles = 0;
2777: outputCtx->trianglelist = NULL;
2778: outputCtx->triangleattributelist = NULL;
2779: outputCtx->neighborlist = NULL;
2780: outputCtx->segmentlist = NULL;
2781: outputCtx->segmentmarkerlist = NULL;
2782: outputCtx->numberofedges = 0;
2783: outputCtx->edgelist = NULL;
2784: outputCtx->edgemarkerlist = NULL;
2785: return(0);
2786: }
2790: PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
2791: {
2793: free(outputCtx->pointlist);
2794: free(outputCtx->pointmarkerlist);
2795: free(outputCtx->segmentlist);
2796: free(outputCtx->segmentmarkerlist);
2797: free(outputCtx->edgelist);
2798: free(outputCtx->edgemarkerlist);
2799: free(outputCtx->trianglelist);
2800: free(outputCtx->neighborlist);
2801: return(0);
2802: }
2806: PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
2807: {
2808: MPI_Comm comm;
2809: PetscInt dim = 2;
2810: const PetscBool createConvexHull = PETSC_FALSE;
2811: const PetscBool constrained = PETSC_FALSE;
2812: struct triangulateio in;
2813: struct triangulateio out;
2814: PetscInt vStart, vEnd, v, eStart, eEnd, e;
2815: PetscMPIInt rank;
2816: PetscErrorCode ierr;
2819: PetscObjectGetComm((PetscObject)boundary,&comm);
2820: MPI_Comm_rank(comm, &rank);
2821: InitInput_Triangle(&in);
2822: InitOutput_Triangle(&out);
2823: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
2825: in.numberofpoints = vEnd - vStart;
2826: if (in.numberofpoints > 0) {
2827: PetscSection coordSection;
2828: Vec coordinates;
2829: PetscScalar *array;
2831: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
2832: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
2833: DMGetCoordinatesLocal(boundary, &coordinates);
2834: DMGetCoordinateSection(boundary, &coordSection);
2835: VecGetArray(coordinates, &array);
2836: for (v = vStart; v < vEnd; ++v) {
2837: const PetscInt idx = v - vStart;
2838: PetscInt off, d;
2840: PetscSectionGetOffset(coordSection, v, &off);
2841: for (d = 0; d < dim; ++d) {
2842: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
2843: }
2844: DMPlexGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);
2845: }
2846: VecRestoreArray(coordinates, &array);
2847: }
2848: DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
2849: in.numberofsegments = eEnd - eStart;
2850: if (in.numberofsegments > 0) {
2851: PetscMalloc1(in.numberofsegments*2, &in.segmentlist);
2852: PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);
2853: for (e = eStart; e < eEnd; ++e) {
2854: const PetscInt idx = e - eStart;
2855: const PetscInt *cone;
2857: DMPlexGetCone(boundary, e, &cone);
2859: in.segmentlist[idx*2+0] = cone[0] - vStart;
2860: in.segmentlist[idx*2+1] = cone[1] - vStart;
2862: DMPlexGetLabelValue(boundary, "marker", e, &in.segmentmarkerlist[idx]);
2863: }
2864: }
2865: #if 0 /* Do not currently support holes */
2866: PetscReal *holeCoords;
2867: PetscInt h, d;
2869: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
2870: if (in.numberofholes > 0) {
2871: PetscMalloc1(in.numberofholes*dim, &in.holelist);
2872: for (h = 0; h < in.numberofholes; ++h) {
2873: for (d = 0; d < dim; ++d) {
2874: in.holelist[h*dim+d] = holeCoords[h*dim+d];
2875: }
2876: }
2877: }
2878: #endif
2879: if (!rank) {
2880: char args[32];
2882: /* Take away 'Q' for verbose output */
2883: PetscStrcpy(args, "pqezQ");
2884: if (createConvexHull) {
2885: PetscStrcat(args, "c");
2886: }
2887: if (constrained) {
2888: PetscStrcpy(args, "zepDQ");
2889: }
2890: triangulate(args, &in, &out, NULL);
2891: }
2892: PetscFree(in.pointlist);
2893: PetscFree(in.pointmarkerlist);
2894: PetscFree(in.segmentlist);
2895: PetscFree(in.segmentmarkerlist);
2896: PetscFree(in.holelist);
2898: {
2899: const PetscInt numCorners = 3;
2900: const PetscInt numCells = out.numberoftriangles;
2901: const PetscInt numVertices = out.numberofpoints;
2902: const int *cells = out.trianglelist;
2903: const double *meshCoords = out.pointlist;
2905: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
2906: /* Set labels */
2907: for (v = 0; v < numVertices; ++v) {
2908: if (out.pointmarkerlist[v]) {
2909: DMPlexSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);
2910: }
2911: }
2912: if (interpolate) {
2913: for (e = 0; e < out.numberofedges; e++) {
2914: if (out.edgemarkerlist[e]) {
2915: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
2916: const PetscInt *edges;
2917: PetscInt numEdges;
2919: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
2920: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
2921: DMPlexSetLabelValue(*dm, "marker", edges[0], out.edgemarkerlist[e]);
2922: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
2923: }
2924: }
2925: }
2926: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
2927: }
2928: #if 0 /* Do not currently support holes */
2929: DMPlexCopyHoles(*dm, boundary);
2930: #endif
2931: FiniOutput_Triangle(&out);
2932: return(0);
2933: }
2937: PetscErrorCode DMPlexRefine_Triangle(DM dm, double *maxVolumes, DM *dmRefined)
2938: {
2939: MPI_Comm comm;
2940: PetscInt dim = 2;
2941: struct triangulateio in;
2942: struct triangulateio out;
2943: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
2944: PetscMPIInt rank;
2945: PetscErrorCode ierr;
2948: PetscObjectGetComm((PetscObject)dm,&comm);
2949: MPI_Comm_rank(comm, &rank);
2950: InitInput_Triangle(&in);
2951: InitOutput_Triangle(&out);
2952: DMPlexGetDepth(dm, &depth);
2953: MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
2954: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
2956: in.numberofpoints = vEnd - vStart;
2957: if (in.numberofpoints > 0) {
2958: PetscSection coordSection;
2959: Vec coordinates;
2960: PetscScalar *array;
2962: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
2963: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
2964: DMGetCoordinatesLocal(dm, &coordinates);
2965: DMGetCoordinateSection(dm, &coordSection);
2966: VecGetArray(coordinates, &array);
2967: for (v = vStart; v < vEnd; ++v) {
2968: const PetscInt idx = v - vStart;
2969: PetscInt off, d;
2971: PetscSectionGetOffset(coordSection, v, &off);
2972: for (d = 0; d < dim; ++d) {
2973: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
2974: }
2975: DMPlexGetLabelValue(dm, "marker", v, &in.pointmarkerlist[idx]);
2976: }
2977: VecRestoreArray(coordinates, &array);
2978: }
2979: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2981: in.numberofcorners = 3;
2982: in.numberoftriangles = cEnd - cStart;
2984: in.trianglearealist = (double*) maxVolumes;
2985: if (in.numberoftriangles > 0) {
2986: PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
2987: for (c = cStart; c < cEnd; ++c) {
2988: const PetscInt idx = c - cStart;
2989: PetscInt *closure = NULL;
2990: PetscInt closureSize;
2992: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
2993: if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
2994: for (v = 0; v < 3; ++v) {
2995: in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
2996: }
2997: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
2998: }
2999: }
3000: /* TODO: Segment markers are missing on input */
3001: #if 0 /* Do not currently support holes */
3002: PetscReal *holeCoords;
3003: PetscInt h, d;
3005: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
3006: if (in.numberofholes > 0) {
3007: PetscMalloc1(in.numberofholes*dim, &in.holelist);
3008: for (h = 0; h < in.numberofholes; ++h) {
3009: for (d = 0; d < dim; ++d) {
3010: in.holelist[h*dim+d] = holeCoords[h*dim+d];
3011: }
3012: }
3013: }
3014: #endif
3015: if (!rank) {
3016: char args[32];
3018: /* Take away 'Q' for verbose output */
3019: PetscStrcpy(args, "pqezQra");
3020: triangulate(args, &in, &out, NULL);
3021: }
3022: PetscFree(in.pointlist);
3023: PetscFree(in.pointmarkerlist);
3024: PetscFree(in.segmentlist);
3025: PetscFree(in.segmentmarkerlist);
3026: PetscFree(in.trianglelist);
3028: {
3029: const PetscInt numCorners = 3;
3030: const PetscInt numCells = out.numberoftriangles;
3031: const PetscInt numVertices = out.numberofpoints;
3032: const int *cells = out.trianglelist;
3033: const double *meshCoords = out.pointlist;
3034: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
3036: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
3037: /* Set labels */
3038: for (v = 0; v < numVertices; ++v) {
3039: if (out.pointmarkerlist[v]) {
3040: DMPlexSetLabelValue(*dmRefined, "marker", v+numCells, out.pointmarkerlist[v]);
3041: }
3042: }
3043: if (interpolate) {
3044: PetscInt e;
3046: for (e = 0; e < out.numberofedges; e++) {
3047: if (out.edgemarkerlist[e]) {
3048: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3049: const PetscInt *edges;
3050: PetscInt numEdges;
3052: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3053: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3054: DMPlexSetLabelValue(*dmRefined, "marker", edges[0], out.edgemarkerlist[e]);
3055: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3056: }
3057: }
3058: }
3059: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
3060: }
3061: #if 0 /* Do not currently support holes */
3062: DMPlexCopyHoles(*dm, boundary);
3063: #endif
3064: FiniOutput_Triangle(&out);
3065: return(0);
3066: }
3067: #endif
3069: #if defined(PETSC_HAVE_TETGEN)
3070: #include <tetgen.h>
3073: PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
3074: {
3075: MPI_Comm comm;
3076: const PetscInt dim = 3;
3077: ::tetgenio in;
3078: ::tetgenio out;
3079: PetscInt vStart, vEnd, v, fStart, fEnd, f;
3080: PetscMPIInt rank;
3084: PetscObjectGetComm((PetscObject)boundary,&comm);
3085: MPI_Comm_rank(comm, &rank);
3086: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
3087: in.numberofpoints = vEnd - vStart;
3088: if (in.numberofpoints > 0) {
3089: PetscSection coordSection;
3090: Vec coordinates;
3091: PetscScalar *array;
3093: in.pointlist = new double[in.numberofpoints*dim];
3094: in.pointmarkerlist = new int[in.numberofpoints];
3096: DMGetCoordinatesLocal(boundary, &coordinates);
3097: DMGetCoordinateSection(boundary, &coordSection);
3098: VecGetArray(coordinates, &array);
3099: for (v = vStart; v < vEnd; ++v) {
3100: const PetscInt idx = v - vStart;
3101: PetscInt off, d;
3103: PetscSectionGetOffset(coordSection, v, &off);
3104: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3105: DMPlexGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);
3106: }
3107: VecRestoreArray(coordinates, &array);
3108: }
3109: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
3111: in.numberoffacets = fEnd - fStart;
3112: if (in.numberoffacets > 0) {
3113: in.facetlist = new tetgenio::facet[in.numberoffacets];
3114: in.facetmarkerlist = new int[in.numberoffacets];
3115: for (f = fStart; f < fEnd; ++f) {
3116: const PetscInt idx = f - fStart;
3117: PetscInt *points = NULL, numPoints, p, numVertices = 0, v;
3119: in.facetlist[idx].numberofpolygons = 1;
3120: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
3121: in.facetlist[idx].numberofholes = 0;
3122: in.facetlist[idx].holelist = NULL;
3124: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3125: for (p = 0; p < numPoints*2; p += 2) {
3126: const PetscInt point = points[p];
3127: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
3128: }
3130: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
3131: poly->numberofvertices = numVertices;
3132: poly->vertexlist = new int[poly->numberofvertices];
3133: for (v = 0; v < numVertices; ++v) {
3134: const PetscInt vIdx = points[v] - vStart;
3135: poly->vertexlist[v] = vIdx;
3136: }
3137: DMPlexGetLabelValue(boundary, "marker", f, &in.facetmarkerlist[idx]);
3138: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3139: }
3140: }
3141: if (!rank) {
3142: char args[32];
3144: /* Take away 'Q' for verbose output */
3145: PetscStrcpy(args, "pqezQ");
3146: ::tetrahedralize(args, &in, &out);
3147: }
3148: {
3149: const PetscInt numCorners = 4;
3150: const PetscInt numCells = out.numberoftetrahedra;
3151: const PetscInt numVertices = out.numberofpoints;
3152: const double *meshCoords = out.pointlist;
3153: int *cells = out.tetrahedronlist;
3155: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
3156: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
3157: /* Set labels */
3158: for (v = 0; v < numVertices; ++v) {
3159: if (out.pointmarkerlist[v]) {
3160: DMPlexSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);
3161: }
3162: }
3163: if (interpolate) {
3164: PetscInt e;
3166: for (e = 0; e < out.numberofedges; e++) {
3167: if (out.edgemarkerlist[e]) {
3168: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3169: const PetscInt *edges;
3170: PetscInt numEdges;
3172: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
3173: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3174: DMPlexSetLabelValue(*dm, "marker", edges[0], out.edgemarkerlist[e]);
3175: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
3176: }
3177: }
3178: for (f = 0; f < out.numberoftrifaces; f++) {
3179: if (out.trifacemarkerlist[f]) {
3180: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
3181: const PetscInt *faces;
3182: PetscInt numFaces;
3184: DMPlexGetJoin(*dm, 3, vertices, &numFaces, &faces);
3185: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3186: DMPlexSetLabelValue(*dm, "marker", faces[0], out.trifacemarkerlist[f]);
3187: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
3188: }
3189: }
3190: }
3191: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
3192: }
3193: return(0);
3194: }
3198: PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
3199: {
3200: MPI_Comm comm;
3201: const PetscInt dim = 3;
3202: ::tetgenio in;
3203: ::tetgenio out;
3204: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
3205: PetscMPIInt rank;
3209: PetscObjectGetComm((PetscObject)dm,&comm);
3210: MPI_Comm_rank(comm, &rank);
3211: DMPlexGetDepth(dm, &depth);
3212: MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
3213: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3215: in.numberofpoints = vEnd - vStart;
3216: if (in.numberofpoints > 0) {
3217: PetscSection coordSection;
3218: Vec coordinates;
3219: PetscScalar *array;
3221: in.pointlist = new double[in.numberofpoints*dim];
3222: in.pointmarkerlist = new int[in.numberofpoints];
3224: DMGetCoordinatesLocal(dm, &coordinates);
3225: DMGetCoordinateSection(dm, &coordSection);
3226: VecGetArray(coordinates, &array);
3227: for (v = vStart; v < vEnd; ++v) {
3228: const PetscInt idx = v - vStart;
3229: PetscInt off, d;
3231: PetscSectionGetOffset(coordSection, v, &off);
3232: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3233: DMPlexGetLabelValue(dm, "marker", v, &in.pointmarkerlist[idx]);
3234: }
3235: VecRestoreArray(coordinates, &array);
3236: }
3237: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
3239: in.numberofcorners = 4;
3240: in.numberoftetrahedra = cEnd - cStart;
3241: in.tetrahedronvolumelist = (double*) maxVolumes;
3242: if (in.numberoftetrahedra > 0) {
3243: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
3244: for (c = cStart; c < cEnd; ++c) {
3245: const PetscInt idx = c - cStart;
3246: PetscInt *closure = NULL;
3247: PetscInt closureSize;
3249: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3250: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
3251: for (v = 0; v < 4; ++v) {
3252: in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
3253: }
3254: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3255: }
3256: }
3257: /* TODO: Put in boundary faces with markers */
3258: if (!rank) {
3259: char args[32];
3261: /* Take away 'Q' for verbose output */
3262: /*PetscStrcpy(args, "qezQra"); */
3263: PetscStrcpy(args, "qezraVVVV");
3264: ::tetrahedralize(args, &in, &out);
3265: }
3266: in.tetrahedronvolumelist = NULL;
3268: {
3269: const PetscInt numCorners = 4;
3270: const PetscInt numCells = out.numberoftetrahedra;
3271: const PetscInt numVertices = out.numberofpoints;
3272: const double *meshCoords = out.pointlist;
3273: int *cells = out.tetrahedronlist;
3275: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
3277: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
3278: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
3279: /* Set labels */
3280: for (v = 0; v < numVertices; ++v) {
3281: if (out.pointmarkerlist[v]) {
3282: DMPlexSetLabelValue(*dmRefined, "marker", v+numCells, out.pointmarkerlist[v]);
3283: }
3284: }
3285: if (interpolate) {
3286: PetscInt e, f;
3288: for (e = 0; e < out.numberofedges; e++) {
3289: if (out.edgemarkerlist[e]) {
3290: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3291: const PetscInt *edges;
3292: PetscInt numEdges;
3294: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3295: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3296: DMPlexSetLabelValue(*dmRefined, "marker", edges[0], out.edgemarkerlist[e]);
3297: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3298: }
3299: }
3300: for (f = 0; f < out.numberoftrifaces; f++) {
3301: if (out.trifacemarkerlist[f]) {
3302: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
3303: const PetscInt *faces;
3304: PetscInt numFaces;
3306: DMPlexGetJoin(*dmRefined, 3, vertices, &numFaces, &faces);
3307: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3308: DMPlexSetLabelValue(*dmRefined, "marker", faces[0], out.trifacemarkerlist[f]);
3309: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
3310: }
3311: }
3312: }
3313: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
3314: }
3315: return(0);
3316: }
3317: #endif
3319: #if defined(PETSC_HAVE_CTETGEN)
3320: #include <ctetgen.h>
3324: PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
3325: {
3326: MPI_Comm comm;
3327: const PetscInt dim = 3;
3328: PLC *in, *out;
3329: PetscInt verbose = 0, vStart, vEnd, v, fStart, fEnd, f;
3330: PetscMPIInt rank;
3334: PetscObjectGetComm((PetscObject)boundary,&comm);
3335: PetscOptionsGetInt(((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
3336: MPI_Comm_rank(comm, &rank);
3337: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
3338: PLCCreate(&in);
3339: PLCCreate(&out);
3341: in->numberofpoints = vEnd - vStart;
3342: if (in->numberofpoints > 0) {
3343: PetscSection coordSection;
3344: Vec coordinates;
3345: PetscScalar *array;
3347: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
3348: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
3349: DMGetCoordinatesLocal(boundary, &coordinates);
3350: DMGetCoordinateSection(boundary, &coordSection);
3351: VecGetArray(coordinates, &array);
3352: for (v = vStart; v < vEnd; ++v) {
3353: const PetscInt idx = v - vStart;
3354: PetscInt off, d, m;
3356: PetscSectionGetOffset(coordSection, v, &off);
3357: for (d = 0; d < dim; ++d) {
3358: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3359: }
3360: DMPlexGetLabelValue(boundary, "marker", v, &m);
3362: in->pointmarkerlist[idx] = (int) m;
3363: }
3364: VecRestoreArray(coordinates, &array);
3365: }
3366: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
3368: in->numberoffacets = fEnd - fStart;
3369: if (in->numberoffacets > 0) {
3370: PetscMalloc1(in->numberoffacets, &in->facetlist);
3371: PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);
3372: for (f = fStart; f < fEnd; ++f) {
3373: const PetscInt idx = f - fStart;
3374: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, m;
3375: polygon *poly;
3377: in->facetlist[idx].numberofpolygons = 1;
3379: PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);
3381: in->facetlist[idx].numberofholes = 0;
3382: in->facetlist[idx].holelist = NULL;
3384: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3385: for (p = 0; p < numPoints*2; p += 2) {
3386: const PetscInt point = points[p];
3387: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
3388: }
3390: poly = in->facetlist[idx].polygonlist;
3391: poly->numberofvertices = numVertices;
3392: PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
3393: for (v = 0; v < numVertices; ++v) {
3394: const PetscInt vIdx = points[v] - vStart;
3395: poly->vertexlist[v] = vIdx;
3396: }
3397: DMPlexGetLabelValue(boundary, "marker", f, &m);
3398: in->facetmarkerlist[idx] = (int) m;
3399: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3400: }
3401: }
3402: if (!rank) {
3403: TetGenOpts t;
3405: TetGenOptsInitialize(&t);
3406: t.in = boundary; /* Should go away */
3407: t.plc = 1;
3408: t.quality = 1;
3409: t.edgesout = 1;
3410: t.zeroindex = 1;
3411: t.quiet = 1;
3412: t.verbose = verbose;
3413: TetGenCheckOpts(&t);
3414: TetGenTetrahedralize(&t, in, out);
3415: }
3416: {
3417: const PetscInt numCorners = 4;
3418: const PetscInt numCells = out->numberoftetrahedra;
3419: const PetscInt numVertices = out->numberofpoints;
3420: const double *meshCoords = out->pointlist;
3421: int *cells = out->tetrahedronlist;
3423: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
3424: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
3425: /* Set labels */
3426: for (v = 0; v < numVertices; ++v) {
3427: if (out->pointmarkerlist[v]) {
3428: DMPlexSetLabelValue(*dm, "marker", v+numCells, out->pointmarkerlist[v]);
3429: }
3430: }
3431: if (interpolate) {
3432: PetscInt e;
3434: for (e = 0; e < out->numberofedges; e++) {
3435: if (out->edgemarkerlist[e]) {
3436: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
3437: const PetscInt *edges;
3438: PetscInt numEdges;
3440: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
3441: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3442: DMPlexSetLabelValue(*dm, "marker", edges[0], out->edgemarkerlist[e]);
3443: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
3444: }
3445: }
3446: for (f = 0; f < out->numberoftrifaces; f++) {
3447: if (out->trifacemarkerlist[f]) {
3448: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
3449: const PetscInt *faces;
3450: PetscInt numFaces;
3452: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
3453: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3454: DMPlexSetLabelValue(*dm, "marker", faces[0], out->trifacemarkerlist[f]);
3455: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
3456: }
3457: }
3458: }
3459: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
3460: }
3462: PLCDestroy(&in);
3463: PLCDestroy(&out);
3464: return(0);
3465: }
3469: PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
3470: {
3471: MPI_Comm comm;
3472: const PetscInt dim = 3;
3473: PLC *in, *out;
3474: PetscInt verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
3475: PetscMPIInt rank;
3479: PetscObjectGetComm((PetscObject)dm,&comm);
3480: PetscOptionsGetInt(((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
3481: MPI_Comm_rank(comm, &rank);
3482: DMPlexGetDepth(dm, &depth);
3483: MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
3484: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3485: PLCCreate(&in);
3486: PLCCreate(&out);
3488: in->numberofpoints = vEnd - vStart;
3489: if (in->numberofpoints > 0) {
3490: PetscSection coordSection;
3491: Vec coordinates;
3492: PetscScalar *array;
3494: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
3495: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
3496: DMGetCoordinatesLocal(dm, &coordinates);
3497: DMGetCoordinateSection(dm, &coordSection);
3498: VecGetArray(coordinates, &array);
3499: for (v = vStart; v < vEnd; ++v) {
3500: const PetscInt idx = v - vStart;
3501: PetscInt off, d, m;
3503: PetscSectionGetOffset(coordSection, v, &off);
3504: for (d = 0; d < dim; ++d) {
3505: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3506: }
3507: DMPlexGetLabelValue(dm, "marker", v, &m);
3509: in->pointmarkerlist[idx] = (int) m;
3510: }
3511: VecRestoreArray(coordinates, &array);
3512: }
3513: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
3515: in->numberofcorners = 4;
3516: in->numberoftetrahedra = cEnd - cStart;
3517: in->tetrahedronvolumelist = maxVolumes;
3518: if (in->numberoftetrahedra > 0) {
3519: PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
3520: for (c = cStart; c < cEnd; ++c) {
3521: const PetscInt idx = c - cStart;
3522: PetscInt *closure = NULL;
3523: PetscInt closureSize;
3525: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3526: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
3527: for (v = 0; v < 4; ++v) {
3528: in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
3529: }
3530: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3531: }
3532: }
3533: if (!rank) {
3534: TetGenOpts t;
3536: TetGenOptsInitialize(&t);
3538: t.in = dm; /* Should go away */
3539: t.refine = 1;
3540: t.varvolume = 1;
3541: t.quality = 1;
3542: t.edgesout = 1;
3543: t.zeroindex = 1;
3544: t.quiet = 1;
3545: t.verbose = verbose; /* Change this */
3547: TetGenCheckOpts(&t);
3548: TetGenTetrahedralize(&t, in, out);
3549: }
3550: {
3551: const PetscInt numCorners = 4;
3552: const PetscInt numCells = out->numberoftetrahedra;
3553: const PetscInt numVertices = out->numberofpoints;
3554: const double *meshCoords = out->pointlist;
3555: int *cells = out->tetrahedronlist;
3556: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
3558: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
3559: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
3560: /* Set labels */
3561: for (v = 0; v < numVertices; ++v) {
3562: if (out->pointmarkerlist[v]) {
3563: DMPlexSetLabelValue(*dmRefined, "marker", v+numCells, out->pointmarkerlist[v]);
3564: }
3565: }
3566: if (interpolate) {
3567: PetscInt e, f;
3569: for (e = 0; e < out->numberofedges; e++) {
3570: if (out->edgemarkerlist[e]) {
3571: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
3572: const PetscInt *edges;
3573: PetscInt numEdges;
3575: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3576: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3577: DMPlexSetLabelValue(*dmRefined, "marker", edges[0], out->edgemarkerlist[e]);
3578: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
3579: }
3580: }
3581: for (f = 0; f < out->numberoftrifaces; f++) {
3582: if (out->trifacemarkerlist[f]) {
3583: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
3584: const PetscInt *faces;
3585: PetscInt numFaces;
3587: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
3588: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3589: DMPlexSetLabelValue(*dmRefined, "marker", faces[0], out->trifacemarkerlist[f]);
3590: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
3591: }
3592: }
3593: }
3594: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
3595: }
3596: PLCDestroy(&in);
3597: PLCDestroy(&out);
3598: return(0);
3599: }
3600: #endif
3604: /*@C
3605: DMPlexGenerate - Generates a mesh.
3607: Not Collective
3609: Input Parameters:
3610: + boundary - The DMPlex boundary object
3611: . name - The mesh generation package name
3612: - interpolate - Flag to create intermediate mesh elements
3614: Output Parameter:
3615: . mesh - The DMPlex object
3617: Level: intermediate
3619: .keywords: mesh, elements
3620: .seealso: DMPlexCreate(), DMRefine()
3621: @*/
3622: PetscErrorCode DMPlexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
3623: {
3624: PetscInt dim;
3625: char genname[1024];
3626: PetscBool isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
3632: DMPlexGetDimension(boundary, &dim);
3633: PetscOptionsGetString(((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);
3634: if (flg) name = genname;
3635: if (name) {
3636: PetscStrcmp(name, "triangle", &isTriangle);
3637: PetscStrcmp(name, "tetgen", &isTetgen);
3638: PetscStrcmp(name, "ctetgen", &isCTetgen);
3639: }
3640: switch (dim) {
3641: case 1:
3642: if (!name || isTriangle) {
3643: #if defined(PETSC_HAVE_TRIANGLE)
3644: DMPlexGenerate_Triangle(boundary, interpolate, mesh);
3645: #else
3646: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
3647: #endif
3648: } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
3649: break;
3650: case 2:
3651: if (!name || isCTetgen) {
3652: #if defined(PETSC_HAVE_CTETGEN)
3653: DMPlexGenerate_CTetgen(boundary, interpolate, mesh);
3654: #else
3655: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
3656: #endif
3657: } else if (isTetgen) {
3658: #if defined(PETSC_HAVE_TETGEN)
3659: DMPlexGenerate_Tetgen(boundary, interpolate, mesh);
3660: #else
3661: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
3662: #endif
3663: } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
3664: break;
3665: default:
3666: SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
3667: }
3668: return(0);
3669: }
3673: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
3674: {
3675: PetscReal refinementLimit;
3676: PetscInt dim, cStart, cEnd;
3677: char genname[1024], *name = NULL;
3678: PetscBool isUniform, isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
3682: DMPlexGetRefinementUniform(dm, &isUniform);
3683: if (isUniform) {
3684: CellRefiner cellRefiner;
3686: DMPlexGetCellRefiner_Internal(dm, &cellRefiner);
3687: DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);
3688: DMPlexCopyBoundary(dm, *dmRefined);
3689: return(0);
3690: }
3691: DMPlexGetRefinementLimit(dm, &refinementLimit);
3692: if (refinementLimit == 0.0) return(0);
3693: DMPlexGetDimension(dm, &dim);
3694: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
3695: PetscOptionsGetString(((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);
3696: if (flg) name = genname;
3697: if (name) {
3698: PetscStrcmp(name, "triangle", &isTriangle);
3699: PetscStrcmp(name, "tetgen", &isTetgen);
3700: PetscStrcmp(name, "ctetgen", &isCTetgen);
3701: }
3702: switch (dim) {
3703: case 2:
3704: if (!name || isTriangle) {
3705: #if defined(PETSC_HAVE_TRIANGLE)
3706: double *maxVolumes;
3707: PetscInt c;
3709: PetscMalloc1((cEnd - cStart), &maxVolumes);
3710: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
3711: DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);
3712: PetscFree(maxVolumes);
3713: #else
3714: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
3715: #endif
3716: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
3717: break;
3718: case 3:
3719: if (!name || isCTetgen) {
3720: #if defined(PETSC_HAVE_CTETGEN)
3721: PetscReal *maxVolumes;
3722: PetscInt c;
3724: PetscMalloc1((cEnd - cStart), &maxVolumes);
3725: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
3726: DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);
3727: #else
3728: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
3729: #endif
3730: } else if (isTetgen) {
3731: #if defined(PETSC_HAVE_TETGEN)
3732: double *maxVolumes;
3733: PetscInt c;
3735: PetscMalloc1((cEnd - cStart), &maxVolumes);
3736: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
3737: DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);
3738: #else
3739: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
3740: #endif
3741: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
3742: break;
3743: default:
3744: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
3745: }
3746: DMPlexCopyBoundary(dm, *dmRefined);
3747: return(0);
3748: }
3752: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
3753: {
3754: DM cdm = dm;
3755: PetscInt r;
3756: PetscBool isUniform;
3760: DMPlexGetRefinementUniform(dm, &isUniform);
3761: if (!isUniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Non-uniform refinement is incompatible with the hierarchy");
3762: for (r = 0; r < nlevels; ++r) {
3763: CellRefiner cellRefiner;
3765: DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);
3766: DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);
3767: DMPlexSetCoarseDM(dmRefined[r], cdm);
3768: cdm = dmRefined[r];
3769: }
3770: return(0);
3771: }
3775: PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened)
3776: {
3777: DM_Plex *mesh = (DM_Plex*) dm->data;
3781: PetscObjectReference((PetscObject) mesh->coarseMesh);
3782: *dmCoarsened = mesh->coarseMesh;
3783: return(0);
3784: }
3788: PetscErrorCode DMPlexLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
3789: {
3790: PetscInt d;
3793: if (!dm->maxCell) {
3794: for (d = 0; d < dim; ++d) out[d] = in[d];
3795: } else {
3796: for (d = 0; d < dim; ++d) {
3797: if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
3798: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
3799: } else {
3800: out[d] = in[d];
3801: }
3802: }
3803: }
3804: return(0);
3805: }
3809: PetscErrorCode DMPlexLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
3810: {
3811: PetscInt d;
3814: if (!dm->maxCell) {
3815: for (d = 0; d < dim; ++d) out[d] += in[d];
3816: } else {
3817: for (d = 0; d < dim; ++d) {
3818: if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
3819: out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
3820: } else {
3821: out[d] += in[d];
3822: }
3823: }
3824: }
3825: return(0);
3826: }
3830: PetscErrorCode DMPlexLocalizeCoordinates(DM dm)
3831: {
3832: PetscSection coordSection, cSection;
3833: Vec coordinates, cVec;
3834: PetscScalar *coords, *coords2, *anchor;
3835: PetscInt Nc, cStart, cEnd, c, vStart, vEnd, v, dof, d, off, off2, bs, coordSize;
3840: if (!dm->maxCell) return(0);
3841: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
3842: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3843: DMGetCoordinatesLocal(dm, &coordinates);
3844: DMGetCoordinateSection(dm, &coordSection);
3845: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
3846: PetscSectionSetNumFields(cSection, 1);
3847: PetscSectionGetFieldComponents(coordSection, 0, &Nc);
3848: PetscSectionSetFieldComponents(cSection, 0, Nc);
3849: PetscSectionSetChart(cSection, cStart, vEnd);
3850: for (v = vStart; v < vEnd; ++v) {
3851: PetscSectionGetDof(coordSection, v, &dof);
3852: PetscSectionSetDof(cSection, v, dof);
3853: PetscSectionSetFieldDof(cSection, v, 0, dof);
3854: }
3855: for (c = cStart; c < cEnd; ++c) {
3856: DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, NULL);
3857: PetscSectionSetDof(cSection, c, dof);
3858: PetscSectionSetFieldDof(cSection, c, 0, dof);
3859: }
3860: PetscSectionSetUp(cSection);
3861: PetscSectionGetStorageSize(cSection, &coordSize);
3862: VecCreate(PetscObjectComm((PetscObject) dm), &cVec);
3863: VecGetBlockSize(coordinates, &bs);
3864: VecSetBlockSize(cVec, bs);
3865: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
3866: VecSetType(cVec,VECSTANDARD);
3867: VecGetArray(coordinates, &coords);
3868: VecGetArray(cVec, &coords2);
3869: for (v = vStart; v < vEnd; ++v) {
3870: PetscSectionGetDof(coordSection, v, &dof);
3871: PetscSectionGetOffset(coordSection, v, &off);
3872: PetscSectionGetOffset(cSection, v, &off2);
3873: for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
3874: }
3875: DMGetWorkArray(dm, 3, PETSC_SCALAR, &anchor);
3876: for (c = cStart; c < cEnd; ++c) {
3877: PetscScalar *cellCoords = NULL;
3878: PetscInt b;
3880: DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
3881: PetscSectionGetOffset(cSection, c, &off2);
3882: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
3883: for (d = 0; d < dof/bs; ++d) {DMPlexLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
3884: DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);
3885: }
3886: DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);
3887: VecRestoreArray(coordinates, &coords);
3888: VecRestoreArray(cVec, &coords2);
3889: DMSetCoordinateSection(dm, cSection);
3890: DMSetCoordinatesLocal(dm, cVec);
3891: VecDestroy(&cVec);
3892: PetscSectionDestroy(&cSection);
3893: return(0);
3894: }
3898: /*@
3899: DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point
3901: Not Collective
3903: Input Parameter:
3904: . dm - The DMPlex object
3906: Output Parameter:
3907: . depthLabel - The DMLabel recording point depth
3909: Level: developer
3911: .keywords: mesh, points
3912: .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
3913: @*/
3914: PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel)
3915: {
3916: DM_Plex *mesh = (DM_Plex*) dm->data;
3922: if (!mesh->depthLabel) {DMPlexGetLabel(dm, "depth", &mesh->depthLabel);}
3923: *depthLabel = mesh->depthLabel;
3924: return(0);
3925: }
3929: /*@
3930: DMPlexGetDepth - Get the depth of the DAG representing this mesh
3932: Not Collective
3934: Input Parameter:
3935: . dm - The DMPlex object
3937: Output Parameter:
3938: . depth - The number of strata (breadth first levels) in the DAG
3940: Level: developer
3942: .keywords: mesh, points
3943: .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum()
3944: @*/
3945: PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth)
3946: {
3947: DMLabel label;
3948: PetscInt d = 0;
3954: DMPlexGetDepthLabel(dm, &label);
3955: if (label) {DMLabelGetNumValues(label, &d);}
3956: *depth = d-1;
3957: return(0);
3958: }
3962: /*@
3963: DMPlexGetDepthStratum - Get the bounds [start, end) for all points at a certain depth.
3965: Not Collective
3967: Input Parameters:
3968: + dm - The DMPlex object
3969: - stratumValue - The requested depth
3971: Output Parameters:
3972: + start - The first point at this depth
3973: - end - One beyond the last point at this depth
3975: Level: developer
3977: .keywords: mesh, points
3978: .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth()
3979: @*/
3980: PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
3981: {
3982: DMLabel label;
3983: PetscInt pStart, pEnd;
3990: DMPlexGetChart(dm, &pStart, &pEnd);
3991: if (pStart == pEnd) return(0);
3992: if (stratumValue < 0) {
3993: if (start) *start = pStart;
3994: if (end) *end = pEnd;
3995: return(0);
3996: }
3997: DMPlexGetDepthLabel(dm, &label);
3998: if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
3999: DMLabelGetStratumBounds(label, stratumValue, start, end);
4000: return(0);
4001: }
4005: /*@
4006: DMPlexGetHeightStratum - Get the bounds [start, end) for all points at a certain height.
4008: Not Collective
4010: Input Parameters:
4011: + dm - The DMPlex object
4012: - stratumValue - The requested height
4014: Output Parameters:
4015: + start - The first point at this height
4016: - end - One beyond the last point at this height
4018: Level: developer
4020: .keywords: mesh, points
4021: .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth()
4022: @*/
4023: PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
4024: {
4025: DMLabel label;
4026: PetscInt depth, pStart, pEnd;
4033: DMPlexGetChart(dm, &pStart, &pEnd);
4034: if (pStart == pEnd) return(0);
4035: if (stratumValue < 0) {
4036: if (start) *start = pStart;
4037: if (end) *end = pEnd;
4038: return(0);
4039: }
4040: DMPlexGetDepthLabel(dm, &label);
4041: if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");
4042: DMLabelGetNumValues(label, &depth);
4043: DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);
4044: return(0);
4045: }
4049: /* Set the number of dof on each point and separate by fields */
4050: PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section)
4051: {
4052: PetscInt *numDofTot;
4053: PetscInt depth, pStart = 0, pEnd = 0;
4054: PetscInt p, d, dep, f;
4058: PetscMalloc1((dim+1), &numDofTot);
4059: for (d = 0; d <= dim; ++d) {
4060: numDofTot[d] = 0;
4061: for (f = 0; f < numFields; ++f) numDofTot[d] += numDof[f*(dim+1)+d];
4062: }
4063: PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);
4064: if (numFields > 0) {
4065: PetscSectionSetNumFields(*section, numFields);
4066: if (numComp) {
4067: for (f = 0; f < numFields; ++f) {
4068: PetscSectionSetFieldComponents(*section, f, numComp[f]);
4069: }
4070: }
4071: }
4072: DMPlexGetChart(dm, &pStart, &pEnd);
4073: PetscSectionSetChart(*section, pStart, pEnd);
4074: DMPlexGetDepth(dm, &depth);
4075: for (dep = 0; dep <= depth; ++dep) {
4076: d = dim == depth ? dep : (!dep ? 0 : dim);
4077: DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);
4078: for (p = pStart; p < pEnd; ++p) {
4079: for (f = 0; f < numFields; ++f) {
4080: PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
4081: }
4082: PetscSectionSetDof(*section, p, numDofTot[d]);
4083: }
4084: }
4085: PetscFree(numDofTot);
4086: return(0);
4087: }
4091: /* Set the number of dof on each point and separate by fields
4092: If constDof is PETSC_DETERMINE, constrain every dof on the point
4093: */
4094: PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC,const PetscInt bcField[],const IS bcPoints[], PetscInt constDof, PetscSection section)
4095: {
4096: PetscInt numFields;
4097: PetscInt bc;
4101: PetscSectionGetNumFields(section, &numFields);
4102: for (bc = 0; bc < numBC; ++bc) {
4103: PetscInt field = 0;
4104: const PetscInt *idx;
4105: PetscInt n, i;
4107: if (numFields) field = bcField[bc];
4108: ISGetLocalSize(bcPoints[bc], &n);
4109: ISGetIndices(bcPoints[bc], &idx);
4110: for (i = 0; i < n; ++i) {
4111: const PetscInt p = idx[i];
4112: PetscInt numConst = constDof;
4114: /* Constrain every dof on the point */
4115: if (numConst < 0) {
4116: if (numFields) {
4117: PetscSectionGetFieldDof(section, p, field, &numConst);
4118: } else {
4119: PetscSectionGetDof(section, p, &numConst);
4120: }
4121: }
4122: if (numFields) {
4123: PetscSectionAddFieldConstraintDof(section, p, field, numConst);
4124: }
4125: PetscSectionAddConstraintDof(section, p, numConst);
4126: }
4127: ISRestoreIndices(bcPoints[bc], &idx);
4128: }
4129: return(0);
4130: }
4134: /* Set the constrained indices on each point and separate by fields */
4135: PetscErrorCode DMPlexCreateSectionBCIndicesAll(DM dm, PetscSection section)
4136: {
4137: PetscInt *maxConstraints;
4138: PetscInt numFields, f, pStart = 0, pEnd = 0, p;
4142: PetscSectionGetNumFields(section, &numFields);
4143: PetscSectionGetChart(section, &pStart, &pEnd);
4144: PetscMalloc1((numFields+1), &maxConstraints);
4145: for (f = 0; f <= numFields; ++f) maxConstraints[f] = 0;
4146: for (p = pStart; p < pEnd; ++p) {
4147: PetscInt cdof;
4149: if (numFields) {
4150: for (f = 0; f < numFields; ++f) {
4151: PetscSectionGetFieldConstraintDof(section, p, f, &cdof);
4152: maxConstraints[f] = PetscMax(maxConstraints[f], cdof);
4153: }
4154: } else {
4155: PetscSectionGetConstraintDof(section, p, &cdof);
4156: maxConstraints[0] = PetscMax(maxConstraints[0], cdof);
4157: }
4158: }
4159: for (f = 0; f < numFields; ++f) {
4160: maxConstraints[numFields] += maxConstraints[f];
4161: }
4162: if (maxConstraints[numFields]) {
4163: PetscInt *indices;
4165: PetscMalloc1(maxConstraints[numFields], &indices);
4166: for (p = pStart; p < pEnd; ++p) {
4167: PetscInt cdof, d;
4169: PetscSectionGetConstraintDof(section, p, &cdof);
4170: if (cdof) {
4171: if (cdof > maxConstraints[numFields]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, "Likely memory corruption, point %D cDof %D > maxConstraints %D", p, cdof, maxConstraints[numFields]);
4172: if (numFields) {
4173: PetscInt numConst = 0, foff = 0;
4175: for (f = 0; f < numFields; ++f) {
4176: PetscInt cfdof, fdof;
4178: PetscSectionGetFieldDof(section, p, f, &fdof);
4179: PetscSectionGetFieldConstraintDof(section, p, f, &cfdof);
4180: /* Change constraint numbering from absolute local dof number to field relative local dof number */
4181: for (d = 0; d < cfdof; ++d) indices[numConst+d] = d;
4182: PetscSectionSetFieldConstraintIndices(section, p, f, &indices[numConst]);
4183: for (d = 0; d < cfdof; ++d) indices[numConst+d] += foff;
4184: numConst += cfdof;
4185: foff += fdof;
4186: }
4187: if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
4188: } else {
4189: for (d = 0; d < cdof; ++d) indices[d] = d;
4190: }
4191: PetscSectionSetConstraintIndices(section, p, indices);
4192: }
4193: }
4194: PetscFree(indices);
4195: }
4196: PetscFree(maxConstraints);
4197: return(0);
4198: }
4202: /* Set the constrained field indices on each point */
4203: PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt field, IS bcPoints, IS constraintIndices, PetscSection section)
4204: {
4205: const PetscInt *points, *indices;
4206: PetscInt numFields, maxDof, numPoints, p, numConstraints;
4207: PetscErrorCode ierr;
4210: PetscSectionGetNumFields(section, &numFields);
4211: if ((field < 0) || (field >= numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, numFields);
4213: ISGetLocalSize(bcPoints, &numPoints);
4214: ISGetIndices(bcPoints, &points);
4215: if (!constraintIndices) {
4216: PetscInt *idx, i;
4218: PetscSectionGetMaxDof(section, &maxDof);
4219: PetscMalloc1(maxDof, &idx);
4220: for (i = 0; i < maxDof; ++i) idx[i] = i;
4221: for (p = 0; p < numPoints; ++p) {
4222: PetscSectionSetFieldConstraintIndices(section, points[p], field, idx);
4223: }
4224: PetscFree(idx);
4225: } else {
4226: ISGetLocalSize(constraintIndices, &numConstraints);
4227: ISGetIndices(constraintIndices, &indices);
4228: for (p = 0; p < numPoints; ++p) {
4229: PetscInt fcdof;
4231: PetscSectionGetFieldConstraintDof(section, points[p], field, &fcdof);
4232: if (fcdof != numConstraints) SETERRQ4(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Section point %d field %d has %d constraints, but yo ugave %d indices", p, field, fcdof, numConstraints);
4233: PetscSectionSetFieldConstraintIndices(section, points[p], field, indices);
4234: }
4235: ISRestoreIndices(constraintIndices, &indices);
4236: }
4237: ISRestoreIndices(bcPoints, &points);
4238: return(0);
4239: }
4243: /* Set the constrained indices on each point and separate by fields */
4244: PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section)
4245: {
4246: PetscInt *indices;
4247: PetscInt numFields, maxDof, f, pStart = 0, pEnd = 0, p;
4251: PetscSectionGetMaxDof(section, &maxDof);
4252: PetscMalloc1(maxDof, &indices);
4253: PetscSectionGetNumFields(section, &numFields);
4254: if (!numFields) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This function only works after users have set field constraint indices.");
4255: PetscSectionGetChart(section, &pStart, &pEnd);
4256: for (p = pStart; p < pEnd; ++p) {
4257: PetscInt cdof, d;
4259: PetscSectionGetConstraintDof(section, p, &cdof);
4260: if (cdof) {
4261: PetscInt numConst = 0, foff = 0;
4263: for (f = 0; f < numFields; ++f) {
4264: const PetscInt *fcind;
4265: PetscInt fdof, fcdof;
4267: PetscSectionGetFieldDof(section, p, f, &fdof);
4268: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
4269: if (fcdof) {PetscSectionGetFieldConstraintIndices(section, p, f, &fcind);}
4270: /* Change constraint numbering from field relative local dof number to absolute local dof number */
4271: for (d = 0; d < fcdof; ++d) indices[numConst+d] = fcind[d]+foff;
4272: foff += fdof;
4273: numConst += fcdof;
4274: }
4275: if (cdof != numConst) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof);
4276: PetscSectionSetConstraintIndices(section, p, indices);
4277: }
4278: }
4279: PetscFree(indices);
4280: return(0);
4281: }
4285: /*@C
4286: DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided.
4288: Not Collective
4290: Input Parameters:
4291: + dm - The DMPlex object
4292: . dim - The spatial dimension of the problem
4293: . numFields - The number of fields in the problem
4294: . numComp - An array of size numFields that holds the number of components for each field
4295: . numDof - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
4296: . numBC - The number of boundary conditions
4297: . bcField - An array of size numBC giving the field number for each boundry condition
4298: . bcPoints - An array of size numBC giving an IS holding the sieve points to which each boundary condition applies
4299: - perm - Optional permutation of the chart, or NULL
4301: Output Parameter:
4302: . section - The PetscSection object
4304: Notes: numDof[f*(dim+1)+d] gives the number of dof for field f on sieve points of dimension d. For instance, numDof[1] is the
4305: number of dof for field 0 on each edge.
4307: The chart permutation is the same one set using PetscSectionSetPermutation()
4309: Level: developer
4311: Fortran Notes:
4312: A Fortran 90 version is available as DMPlexCreateSectionF90()
4314: .keywords: mesh, elements
4315: .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
4316: @*/
4317: PetscErrorCode DMPlexCreateSection(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscInt numBC,const PetscInt bcField[],const IS bcPoints[], IS perm, PetscSection *section)
4318: {
4322: DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);
4323: DMPlexCreateSectionBCDof(dm, numBC, bcField, bcPoints, PETSC_DETERMINE, *section);
4324: if (perm) {PetscSectionSetPermutation(*section, perm);}
4325: PetscSectionSetUp(*section);
4326: if (numBC) {DMPlexCreateSectionBCIndicesAll(dm, *section);}
4327: PetscSectionViewFromOptions(*section,NULL,"-section_view");
4328: return(0);
4329: }
4333: PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm)
4334: {
4335: PetscSection section;
4339: DMClone(dm, cdm);
4340: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
4341: DMSetDefaultSection(*cdm, section);
4342: PetscSectionDestroy(§ion);
4343: return(0);
4344: }
4348: PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section)
4349: {
4350: DM_Plex *mesh = (DM_Plex*) dm->data;
4354: if (section) *section = mesh->coneSection;
4355: return(0);
4356: }
4360: PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section)
4361: {
4362: DM_Plex *mesh = (DM_Plex*) dm->data;
4366: if (section) *section = mesh->supportSection;
4367: return(0);
4368: }
4372: PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[])
4373: {
4374: DM_Plex *mesh = (DM_Plex*) dm->data;
4378: if (cones) *cones = mesh->cones;
4379: return(0);
4380: }
4384: PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[])
4385: {
4386: DM_Plex *mesh = (DM_Plex*) dm->data;
4390: if (coneOrientations) *coneOrientations = mesh->coneOrientations;
4391: return(0);
4392: }
4394: /******************************** FEM Support **********************************/
4398: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
4399: {
4400: PetscScalar *array, *vArray;
4401: const PetscInt *cone, *coneO;
4402: PetscInt pStart, pEnd, p, numPoints, size = 0, offset = 0;
4403: PetscErrorCode ierr;
4406: PetscSectionGetChart(section, &pStart, &pEnd);
4407: DMPlexGetConeSize(dm, point, &numPoints);
4408: DMPlexGetCone(dm, point, &cone);
4409: DMPlexGetConeOrientation(dm, point, &coneO);
4410: if (!values || !*values) {
4411: if ((point >= pStart) && (point < pEnd)) {
4412: PetscInt dof;
4414: PetscSectionGetDof(section, point, &dof);
4415: size += dof;
4416: }
4417: for (p = 0; p < numPoints; ++p) {
4418: const PetscInt cp = cone[p];
4419: PetscInt dof;
4421: if ((cp < pStart) || (cp >= pEnd)) continue;
4422: PetscSectionGetDof(section, cp, &dof);
4423: size += dof;
4424: }
4425: if (!values) {
4426: if (csize) *csize = size;
4427: return(0);
4428: }
4429: DMGetWorkArray(dm, size, PETSC_SCALAR, &array);
4430: } else {
4431: array = *values;
4432: }
4433: size = 0;
4434: VecGetArray(v, &vArray);
4435: if ((point >= pStart) && (point < pEnd)) {
4436: PetscInt dof, off, d;
4437: PetscScalar *varr;
4439: PetscSectionGetDof(section, point, &dof);
4440: PetscSectionGetOffset(section, point, &off);
4441: varr = &vArray[off];
4442: for (d = 0; d < dof; ++d, ++offset) {
4443: array[offset] = varr[d];
4444: }
4445: size += dof;
4446: }
4447: for (p = 0; p < numPoints; ++p) {
4448: const PetscInt cp = cone[p];
4449: PetscInt o = coneO[p];
4450: PetscInt dof, off, d;
4451: PetscScalar *varr;
4453: if ((cp < pStart) || (cp >= pEnd)) continue;
4454: PetscSectionGetDof(section, cp, &dof);
4455: PetscSectionGetOffset(section, cp, &off);
4456: varr = &vArray[off];
4457: if (o >= 0) {
4458: for (d = 0; d < dof; ++d, ++offset) {
4459: array[offset] = varr[d];
4460: }
4461: } else {
4462: for (d = dof-1; d >= 0; --d, ++offset) {
4463: array[offset] = varr[d];
4464: }
4465: }
4466: size += dof;
4467: }
4468: VecRestoreArray(v, &vArray);
4469: if (!*values) {
4470: if (csize) *csize = size;
4471: *values = array;
4472: } else {
4473: if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
4474: *csize = size;
4475: }
4476: return(0);
4477: }
4481: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
4482: {
4483: PetscInt offset = 0, p;
4487: *size = 0;
4488: for (p = 0; p < numPoints*2; p += 2) {
4489: const PetscInt point = points[p];
4490: const PetscInt o = points[p+1];
4491: PetscInt dof, off, d;
4492: const PetscScalar *varr;
4494: PetscSectionGetDof(section, point, &dof);
4495: PetscSectionGetOffset(section, point, &off);
4496: varr = &vArray[off];
4497: if (o >= 0) {
4498: for (d = 0; d < dof; ++d, ++offset) array[offset] = varr[d];
4499: } else {
4500: for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d];
4501: }
4502: }
4503: *size = offset;
4504: return(0);
4505: }
4509: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Fields_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], PetscInt numFields, const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
4510: {
4511: PetscInt offset = 0, f;
4515: *size = 0;
4516: for (f = 0; f < numFields; ++f) {
4517: PetscInt fcomp, p;
4519: PetscSectionGetFieldComponents(section, f, &fcomp);
4520: for (p = 0; p < numPoints*2; p += 2) {
4521: const PetscInt point = points[p];
4522: const PetscInt o = points[p+1];
4523: PetscInt fdof, foff, d, c;
4524: const PetscScalar *varr;
4526: PetscSectionGetFieldDof(section, point, f, &fdof);
4527: PetscSectionGetFieldOffset(section, point, f, &foff);
4528: varr = &vArray[foff];
4529: if (o >= 0) {
4530: for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d];
4531: } else {
4532: for (d = fdof/fcomp-1; d >= 0; --d) {
4533: for (c = 0; c < fcomp; ++c, ++offset) {
4534: array[offset] = varr[d*fcomp+c];
4535: }
4536: }
4537: }
4538: }
4539: }
4540: *size = offset;
4541: return(0);
4542: }
4546: /*@C
4547: DMPlexVecGetClosure - Get an array of the values on the closure of 'point'
4549: Not collective
4551: Input Parameters:
4552: + dm - The DM
4553: . section - The section describing the layout in v, or NULL to use the default section
4554: . v - The local vector
4555: - point - The sieve point in the DM
4557: Output Parameters:
4558: + csize - The number of values in the closure, or NULL
4559: - values - The array of values, which is a borrowed array and should not be freed
4561: Fortran Notes:
4562: Since it returns an array, this routine is only available in Fortran 90, and you must
4563: include petsc.h90 in your code.
4565: The csize argument is not present in the Fortran 90 binding since it is internal to the array.
4567: Level: intermediate
4569: .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
4570: @*/
4571: PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
4572: {
4573: PetscSection clSection;
4574: IS clPoints;
4575: PetscScalar *array, *vArray;
4576: PetscInt *points = NULL;
4577: const PetscInt *clp;
4578: PetscInt depth, numFields, numPoints, size;
4579: PetscErrorCode ierr;
4583: if (!section) {DMGetDefaultSection(dm, §ion);}
4586: DMPlexGetDepth(dm, &depth);
4587: PetscSectionGetNumFields(section, &numFields);
4588: if (depth == 1 && numFields < 2) {
4589: DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);
4590: return(0);
4591: }
4592: /* Get points */
4593: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4594: if (!clPoints) {
4595: PetscInt pStart, pEnd, p, q;
4597: PetscSectionGetChart(section, &pStart, &pEnd);
4598: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4599: /* Compress out points not in the section */
4600: for (p = 0, q = 0; p < numPoints*2; p += 2) {
4601: if ((points[p] >= pStart) && (points[p] < pEnd)) {
4602: points[q*2] = points[p];
4603: points[q*2+1] = points[p+1];
4604: ++q;
4605: }
4606: }
4607: numPoints = q;
4608: } else {
4609: PetscInt dof, off;
4611: PetscSectionGetDof(clSection, point, &dof);
4612: PetscSectionGetOffset(clSection, point, &off);
4613: ISGetIndices(clPoints, &clp);
4614: numPoints = dof/2;
4615: points = (PetscInt *) &clp[off];
4616: }
4617: /* Get array */
4618: if (!values || !*values) {
4619: PetscInt asize = 0, dof, p;
4621: for (p = 0; p < numPoints*2; p += 2) {
4622: PetscSectionGetDof(section, points[p], &dof);
4623: asize += dof;
4624: }
4625: if (!values) {
4626: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4627: else {ISRestoreIndices(clPoints, &clp);}
4628: if (csize) *csize = asize;
4629: return(0);
4630: }
4631: DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);
4632: } else {
4633: array = *values;
4634: }
4635: VecGetArray(v, &vArray);
4636: /* Get values */
4637: if (numFields > 0) {DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);}
4638: else {DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);}
4639: /* Cleanup points */
4640: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
4641: else {ISRestoreIndices(clPoints, &clp);}
4642: /* Cleanup array */
4643: VecRestoreArray(v, &vArray);
4644: if (!*values) {
4645: if (csize) *csize = size;
4646: *values = array;
4647: } else {
4648: if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size);
4649: *csize = size;
4650: }
4651: return(0);
4652: }
4656: /*@C
4657: DMPlexVecRestoreClosure - Restore the array of the values on the closure of 'point'
4659: Not collective
4661: Input Parameters:
4662: + dm - The DM
4663: . section - The section describing the layout in v, or NULL to use the default section
4664: . v - The local vector
4665: . point - The sieve point in the DM
4666: . csize - The number of values in the closure, or NULL
4667: - values - The array of values, which is a borrowed array and should not be freed
4669: Fortran Notes:
4670: Since it returns an array, this routine is only available in Fortran 90, and you must
4671: include petsc.h90 in your code.
4673: The csize argument is not present in the Fortran 90 binding since it is internal to the array.
4675: Level: intermediate
4677: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
4678: @*/
4679: PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
4680: {
4681: PetscInt size = 0;
4685: /* Should work without recalculating size */
4686: DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);
4687: return(0);
4688: }
4690: PETSC_STATIC_INLINE void add (PetscScalar *x, PetscScalar y) {*x += y;}
4691: PETSC_STATIC_INLINE void insert(PetscScalar *x, PetscScalar y) {*x = y;}
4695: PETSC_STATIC_INLINE PetscErrorCode updatePoint_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, PetscInt orientation, const PetscScalar values[], PetscScalar array[])
4696: {
4697: PetscInt cdof; /* The number of constraints on this point */
4698: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
4699: PetscScalar *a;
4700: PetscInt off, cind = 0, k;
4701: PetscErrorCode ierr;
4704: PetscSectionGetConstraintDof(section, point, &cdof);
4705: PetscSectionGetOffset(section, point, &off);
4706: a = &array[off];
4707: if (!cdof || setBC) {
4708: if (orientation >= 0) {
4709: for (k = 0; k < dof; ++k) {
4710: fuse(&a[k], values[k]);
4711: }
4712: } else {
4713: for (k = 0; k < dof; ++k) {
4714: fuse(&a[k], values[dof-k-1]);
4715: }
4716: }
4717: } else {
4718: PetscSectionGetConstraintIndices(section, point, &cdofs);
4719: if (orientation >= 0) {
4720: for (k = 0; k < dof; ++k) {
4721: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4722: fuse(&a[k], values[k]);
4723: }
4724: } else {
4725: for (k = 0; k < dof; ++k) {
4726: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4727: fuse(&a[k], values[dof-k-1]);
4728: }
4729: }
4730: }
4731: return(0);
4732: }
4736: PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[])
4737: {
4738: PetscInt cdof; /* The number of constraints on this point */
4739: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
4740: PetscScalar *a;
4741: PetscInt off, cind = 0, k;
4742: PetscErrorCode ierr;
4745: PetscSectionGetConstraintDof(section, point, &cdof);
4746: PetscSectionGetOffset(section, point, &off);
4747: a = &array[off];
4748: if (cdof) {
4749: PetscSectionGetConstraintIndices(section, point, &cdofs);
4750: if (orientation >= 0) {
4751: for (k = 0; k < dof; ++k) {
4752: if ((cind < cdof) && (k == cdofs[cind])) {
4753: fuse(&a[k], values[k]);
4754: ++cind;
4755: }
4756: }
4757: } else {
4758: for (k = 0; k < dof; ++k) {
4759: if ((cind < cdof) && (k == cdofs[cind])) {
4760: fuse(&a[k], values[dof-k-1]);
4761: ++cind;
4762: }
4763: }
4764: }
4765: }
4766: return(0);
4767: }
4771: PETSC_STATIC_INLINE PetscErrorCode updatePointFields_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, const PetscScalar values[], PetscInt *offset, PetscScalar array[])
4772: {
4773: PetscScalar *a;
4774: PetscInt fdof, foff, fcdof, foffset = *offset;
4775: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4776: PetscInt cind = 0, k, c;
4777: PetscErrorCode ierr;
4780: PetscSectionGetFieldDof(section, point, f, &fdof);
4781: PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
4782: PetscSectionGetFieldOffset(section, point, f, &foff);
4783: a = &array[foff];
4784: if (!fcdof || setBC) {
4785: if (o >= 0) {
4786: for (k = 0; k < fdof; ++k) fuse(&a[k], values[foffset+k]);
4787: } else {
4788: for (k = fdof/fcomp-1; k >= 0; --k) {
4789: for (c = 0; c < fcomp; ++c) {
4790: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
4791: }
4792: }
4793: }
4794: } else {
4795: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
4796: if (o >= 0) {
4797: for (k = 0; k < fdof; ++k) {
4798: if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
4799: fuse(&a[k], values[foffset+k]);
4800: }
4801: } else {
4802: for (k = fdof/fcomp-1; k >= 0; --k) {
4803: for (c = 0; c < fcomp; ++c) {
4804: if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
4805: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
4806: }
4807: }
4808: }
4809: }
4810: *offset += fdof;
4811: return(0);
4812: }
4816: PETSC_STATIC_INLINE PetscErrorCode updatePointFieldsBC_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), const PetscScalar values[], PetscInt *offset, PetscScalar array[])
4817: {
4818: PetscScalar *a;
4819: PetscInt fdof, foff, fcdof, foffset = *offset;
4820: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4821: PetscInt cind = 0, k, c;
4822: PetscErrorCode ierr;
4825: PetscSectionGetFieldDof(section, point, f, &fdof);
4826: PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
4827: PetscSectionGetFieldOffset(section, point, f, &foff);
4828: a = &array[foff];
4829: if (fcdof) {
4830: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
4831: if (o >= 0) {
4832: for (k = 0; k < fdof; ++k) {
4833: if ((cind < fcdof) && (k == fcdofs[cind])) {
4834: fuse(&a[k], values[foffset+k]);
4835: ++cind;
4836: }
4837: }
4838: } else {
4839: for (k = fdof/fcomp-1; k >= 0; --k) {
4840: for (c = 0; c < fcomp; ++c) {
4841: if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {
4842: fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]);
4843: ++cind;
4844: }
4845: }
4846: }
4847: }
4848: }
4849: *offset += fdof;
4850: return(0);
4851: }
4855: PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
4856: {
4857: PetscScalar *array;
4858: const PetscInt *cone, *coneO;
4859: PetscInt pStart, pEnd, p, numPoints, off, dof;
4860: PetscErrorCode ierr;
4863: PetscSectionGetChart(section, &pStart, &pEnd);
4864: DMPlexGetConeSize(dm, point, &numPoints);
4865: DMPlexGetCone(dm, point, &cone);
4866: DMPlexGetConeOrientation(dm, point, &coneO);
4867: VecGetArray(v, &array);
4868: for (p = 0, off = 0; p <= numPoints; ++p, off += dof) {
4869: const PetscInt cp = !p ? point : cone[p-1];
4870: const PetscInt o = !p ? 0 : coneO[p-1];
4872: if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;}
4873: PetscSectionGetDof(section, cp, &dof);
4874: /* ADD_VALUES */
4875: {
4876: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
4877: PetscScalar *a;
4878: PetscInt cdof, coff, cind = 0, k;
4880: PetscSectionGetConstraintDof(section, cp, &cdof);
4881: PetscSectionGetOffset(section, cp, &coff);
4882: a = &array[coff];
4883: if (!cdof) {
4884: if (o >= 0) {
4885: for (k = 0; k < dof; ++k) {
4886: a[k] += values[off+k];
4887: }
4888: } else {
4889: for (k = 0; k < dof; ++k) {
4890: a[k] += values[off+dof-k-1];
4891: }
4892: }
4893: } else {
4894: PetscSectionGetConstraintIndices(section, cp, &cdofs);
4895: if (o >= 0) {
4896: for (k = 0; k < dof; ++k) {
4897: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4898: a[k] += values[off+k];
4899: }
4900: } else {
4901: for (k = 0; k < dof; ++k) {
4902: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4903: a[k] += values[off+dof-k-1];
4904: }
4905: }
4906: }
4907: }
4908: }
4909: VecRestoreArray(v, &array);
4910: return(0);
4911: }
4915: /*@C
4916: DMPlexVecSetClosure - Set an array of the values on the closure of 'point'
4918: Not collective
4920: Input Parameters:
4921: + dm - The DM
4922: . section - The section describing the layout in v, or NULL to use the default section
4923: . v - The local vector
4924: . point - The sieve point in the DM
4925: . values - The array of values
4926: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
4928: Fortran Notes:
4929: This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
4931: Level: intermediate
4933: .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
4934: @*/
4935: PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
4936: {
4937: PetscSection clSection;
4938: IS clPoints;
4939: PetscScalar *array;
4940: PetscInt *points = NULL;
4941: const PetscInt *clp;
4942: PetscInt depth, numFields, numPoints, p;
4943: PetscErrorCode ierr;
4947: if (!section) {DMGetDefaultSection(dm, §ion);}
4950: DMPlexGetDepth(dm, &depth);
4951: PetscSectionGetNumFields(section, &numFields);
4952: if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
4953: DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);
4954: return(0);
4955: }
4956: /* Get points */
4957: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
4958: if (!clPoints) {
4959: PetscInt pStart, pEnd, q;
4961: PetscSectionGetChart(section, &pStart, &pEnd);
4962: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4963: /* Compress out points not in the section */
4964: for (p = 0, q = 0; p < numPoints*2; p += 2) {
4965: if ((points[p] >= pStart) && (points[p] < pEnd)) {
4966: points[q*2] = points[p];
4967: points[q*2+1] = points[p+1];
4968: ++q;
4969: }
4970: }
4971: numPoints = q;
4972: } else {
4973: PetscInt dof, off;
4975: PetscSectionGetDof(clSection, point, &dof);
4976: PetscSectionGetOffset(clSection, point, &off);
4977: ISGetIndices(clPoints, &clp);
4978: numPoints = dof/2;
4979: points = (PetscInt *) &clp[off];
4980: }
4981: /* Get array */
4982: VecGetArray(v, &array);
4983: /* Get values */
4984: if (numFields > 0) {
4985: PetscInt offset = 0, fcomp, f;
4986: for (f = 0; f < numFields; ++f) {
4987: PetscSectionGetFieldComponents(section, f, &fcomp);
4988: switch (mode) {
4989: case INSERT_VALUES:
4990: for (p = 0; p < numPoints*2; p += 2) {
4991: const PetscInt point = points[p];
4992: const PetscInt o = points[p+1];
4993: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
4994: } break;
4995: case INSERT_ALL_VALUES:
4996: for (p = 0; p < numPoints*2; p += 2) {
4997: const PetscInt point = points[p];
4998: const PetscInt o = points[p+1];
4999: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
5000: } break;
5001: case INSERT_BC_VALUES:
5002: for (p = 0; p < numPoints*2; p += 2) {
5003: const PetscInt point = points[p];
5004: const PetscInt o = points[p+1];
5005: updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
5006: } break;
5007: case ADD_VALUES:
5008: for (p = 0; p < numPoints*2; p += 2) {
5009: const PetscInt point = points[p];
5010: const PetscInt o = points[p+1];
5011: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
5012: } break;
5013: case ADD_ALL_VALUES:
5014: for (p = 0; p < numPoints*2; p += 2) {
5015: const PetscInt point = points[p];
5016: const PetscInt o = points[p+1];
5017: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
5018: } break;
5019: default:
5020: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
5021: }
5022: }
5023: } else {
5024: PetscInt dof, off;
5026: switch (mode) {
5027: case INSERT_VALUES:
5028: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5029: PetscInt o = points[p+1];
5030: PetscSectionGetDof(section, points[p], &dof);
5031: updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
5032: } break;
5033: case INSERT_ALL_VALUES:
5034: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5035: PetscInt o = points[p+1];
5036: PetscSectionGetDof(section, points[p], &dof);
5037: updatePoint_private(section, points[p], dof, insert, PETSC_TRUE, o, &values[off], array);
5038: } break;
5039: case INSERT_BC_VALUES:
5040: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5041: PetscInt o = points[p+1];
5042: PetscSectionGetDof(section, points[p], &dof);
5043: updatePointBC_private(section, points[p], dof, insert, o, &values[off], array);
5044: } break;
5045: case ADD_VALUES:
5046: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5047: PetscInt o = points[p+1];
5048: PetscSectionGetDof(section, points[p], &dof);
5049: updatePoint_private(section, points[p], dof, add, PETSC_FALSE, o, &values[off], array);
5050: } break;
5051: case ADD_ALL_VALUES:
5052: for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5053: PetscInt o = points[p+1];
5054: PetscSectionGetDof(section, points[p], &dof);
5055: updatePoint_private(section, points[p], dof, add, PETSC_TRUE, o, &values[off], array);
5056: } break;
5057: default:
5058: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
5059: }
5060: }
5061: /* Cleanup points */
5062: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
5063: else {ISRestoreIndices(clPoints, &clp);}
5064: /* Cleanup array */
5065: VecRestoreArray(v, &array);
5066: return(0);
5067: }
5071: PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode)
5072: {
5073: PetscSection clSection;
5074: IS clPoints;
5075: PetscScalar *array;
5076: PetscInt *points = NULL;
5077: const PetscInt *clp;
5078: PetscInt numFields, numPoints, p;
5079: PetscInt offset = 0, fcomp, f;
5080: PetscErrorCode ierr;
5084: if (!section) {DMGetDefaultSection(dm, §ion);}
5087: PetscSectionGetNumFields(section, &numFields);
5088: /* Get points */
5089: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
5090: if (!clPoints) {
5091: PetscInt pStart, pEnd, q;
5093: PetscSectionGetChart(section, &pStart, &pEnd);
5094: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
5095: /* Compress out points not in the section */
5096: for (p = 0, q = 0; p < numPoints*2; p += 2) {
5097: if ((points[p] >= pStart) && (points[p] < pEnd)) {
5098: points[q*2] = points[p];
5099: points[q*2+1] = points[p+1];
5100: ++q;
5101: }
5102: }
5103: numPoints = q;
5104: } else {
5105: PetscInt dof, off;
5107: PetscSectionGetDof(clSection, point, &dof);
5108: PetscSectionGetOffset(clSection, point, &off);
5109: ISGetIndices(clPoints, &clp);
5110: numPoints = dof/2;
5111: points = (PetscInt *) &clp[off];
5112: }
5113: /* Get array */
5114: VecGetArray(v, &array);
5115: /* Get values */
5116: for (f = 0; f < numFields; ++f) {
5117: PetscSectionGetFieldComponents(section, f, &fcomp);
5118: if (!fieldActive[f]) {
5119: for (p = 0; p < numPoints*2; p += 2) {
5120: PetscInt fdof;
5121: PetscSectionGetFieldDof(section, points[p], f, &fdof);
5122: offset += fdof;
5123: }
5124: continue;
5125: }
5126: switch (mode) {
5127: case INSERT_VALUES:
5128: for (p = 0; p < numPoints*2; p += 2) {
5129: const PetscInt point = points[p];
5130: const PetscInt o = points[p+1];
5131: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array);
5132: } break;
5133: case INSERT_ALL_VALUES:
5134: for (p = 0; p < numPoints*2; p += 2) {
5135: const PetscInt point = points[p];
5136: const PetscInt o = points[p+1];
5137: updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array);
5138: } break;
5139: case INSERT_BC_VALUES:
5140: for (p = 0; p < numPoints*2; p += 2) {
5141: const PetscInt point = points[p];
5142: const PetscInt o = points[p+1];
5143: updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array);
5144: } break;
5145: case ADD_VALUES:
5146: for (p = 0; p < numPoints*2; p += 2) {
5147: const PetscInt point = points[p];
5148: const PetscInt o = points[p+1];
5149: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array);
5150: } break;
5151: case ADD_ALL_VALUES:
5152: for (p = 0; p < numPoints*2; p += 2) {
5153: const PetscInt point = points[p];
5154: const PetscInt o = points[p+1];
5155: updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array);
5156: } break;
5157: default:
5158: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
5159: }
5160: }
5161: /* Cleanup points */
5162: if (!clPoints) {DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);}
5163: else {ISRestoreIndices(clPoints, &clp);}
5164: /* Cleanup array */
5165: VecRestoreArray(v, &array);
5166: return(0);
5167: }
5171: PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[])
5172: {
5173: PetscMPIInt rank;
5174: PetscInt i, j;
5178: MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);
5179: PetscViewerASCIIPrintf(viewer, "[%D]mat for sieve point %D\n", rank, point);
5180: for (i = 0; i < numRIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%D]mat row indices[%D] = %D\n", rank, i, rindices[i]);}
5181: for (i = 0; i < numCIndices; i++) {PetscViewerASCIIPrintf(viewer, "[%D]mat col indices[%D] = %D\n", rank, i, cindices[i]);}
5182: numCIndices = numCIndices ? numCIndices : numRIndices;
5183: for (i = 0; i < numRIndices; i++) {
5184: PetscViewerASCIIPrintf(viewer, "[%D]", rank);
5185: for (j = 0; j < numCIndices; j++) {
5186: #if defined(PETSC_USE_COMPLEX)
5187: PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));
5188: #else
5189: PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);
5190: #endif
5191: }
5192: PetscViewerASCIIPrintf(viewer, "\n");
5193: }
5194: return(0);
5195: }
5199: /* . off - The global offset of this point */
5200: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[])
5201: {
5202: PetscInt dof; /* The number of unknowns on this point */
5203: PetscInt cdof; /* The number of constraints on this point */
5204: const PetscInt *cdofs; /* The indices of the constrained dofs on this point */
5205: PetscInt cind = 0, k;
5206: PetscErrorCode ierr;
5209: PetscSectionGetDof(section, point, &dof);
5210: PetscSectionGetConstraintDof(section, point, &cdof);
5211: if (!cdof || setBC) {
5212: if (orientation >= 0) {
5213: for (k = 0; k < dof; ++k) indices[*loff+k] = off+k;
5214: } else {
5215: for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k;
5216: }
5217: } else {
5218: PetscSectionGetConstraintIndices(section, point, &cdofs);
5219: if (orientation >= 0) {
5220: for (k = 0; k < dof; ++k) {
5221: if ((cind < cdof) && (k == cdofs[cind])) {
5222: /* Insert check for returning constrained indices */
5223: indices[*loff+k] = -(off+k+1);
5224: ++cind;
5225: } else {
5226: indices[*loff+k] = off+k-cind;
5227: }
5228: }
5229: } else {
5230: for (k = 0; k < dof; ++k) {
5231: if ((cind < cdof) && (k == cdofs[cind])) {
5232: /* Insert check for returning constrained indices */
5233: indices[*loff+dof-k-1] = -(off+k+1);
5234: ++cind;
5235: } else {
5236: indices[*loff+dof-k-1] = off+k-cind;
5237: }
5238: }
5239: }
5240: }
5241: *loff += dof;
5242: return(0);
5243: }
5247: /* . off - The global offset of this point */
5248: PetscErrorCode indicesPointFields_private(PetscSection section, PetscInt point, PetscInt off, PetscInt foffs[], PetscBool setBC, PetscInt orientation, PetscInt indices[])
5249: {
5250: PetscInt numFields, foff, f;
5254: PetscSectionGetNumFields(section, &numFields);
5255: for (f = 0, foff = 0; f < numFields; ++f) {
5256: PetscInt fdof, fcomp, cfdof;
5257: const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
5258: PetscInt cind = 0, k, c;
5260: PetscSectionGetFieldComponents(section, f, &fcomp);
5261: PetscSectionGetFieldDof(section, point, f, &fdof);
5262: PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
5263: if (!cfdof || setBC) {
5264: if (orientation >= 0) {
5265: for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k;
5266: } else {
5267: for (k = fdof/fcomp-1; k >= 0; --k) {
5268: for (c = 0; c < fcomp; ++c) {
5269: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
5270: }
5271: }
5272: }
5273: } else {
5274: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
5275: if (orientation >= 0) {
5276: for (k = 0; k < fdof; ++k) {
5277: if ((cind < cfdof) && (k == fcdofs[cind])) {
5278: indices[foffs[f]+k] = -(off+foff+k+1);
5279: ++cind;
5280: } else {
5281: indices[foffs[f]+k] = off+foff+k-cind;
5282: }
5283: }
5284: } else {
5285: for (k = fdof/fcomp-1; k >= 0; --k) {
5286: for (c = 0; c < fcomp; ++c) {
5287: if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
5288: indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
5289: ++cind;
5290: } else {
5291: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind;
5292: }
5293: }
5294: }
5295: }
5296: }
5297: foff += fdof - cfdof;
5298: foffs[f] += fdof;
5299: }
5300: return(0);
5301: }
5305: /*@C
5306: DMPlexMatSetClosure - Set an array of the values on the closure of 'point'
5308: Not collective
5310: Input Parameters:
5311: + dm - The DM
5312: . section - The section describing the layout in v, or NULL to use the default section
5313: . globalSection - The section describing the layout in v, or NULL to use the default global section
5314: . A - The matrix
5315: . point - The sieve point in the DM
5316: . values - The array of values
5317: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
5319: Fortran Notes:
5320: This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
5322: Level: intermediate
5324: .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure()
5325: @*/
5326: PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
5327: {
5328: DM_Plex *mesh = (DM_Plex*) dm->data;
5329: PetscSection clSection;
5330: IS clPoints;
5331: PetscInt *points = NULL;
5332: const PetscInt *clp;
5333: PetscInt *indices;
5334: PetscInt offsets[32];
5335: PetscInt numFields, numPoints, numIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
5336: PetscErrorCode ierr;
5340: if (!section) {DMGetDefaultSection(dm, §ion);}
5342: if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
5345: PetscSectionGetNumFields(section, &numFields);
5346: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
5347: PetscMemzero(offsets, 32 * sizeof(PetscInt));
5348: PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);
5349: if (!clPoints) {
5350: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
5351: /* Compress out points not in the section */
5352: PetscSectionGetChart(section, &pStart, &pEnd);
5353: for (p = 0, q = 0; p < numPoints*2; p += 2) {
5354: if ((points[p] >= pStart) && (points[p] < pEnd)) {
5355: points[q*2] = points[p];
5356: points[q*2+1] = points[p+1];
5357: ++q;
5358: }
5359: }
5360: numPoints = q;
5361: } else {
5362: PetscInt dof, off;
5364: PetscSectionGetDof(clSection, point, &dof);
5365: numPoints = dof/2;
5366: PetscSectionGetOffset(clSection, point, &off);
5367: ISGetIndices(clPoints, &clp);
5368: points = (PetscInt *) &clp[off];
5369: }
5370: for (p = 0, numIndices = 0; p < numPoints*2; p += 2) {
5371: PetscInt fdof;
5373: PetscSectionGetDof(section, points[p], &dof);
5374: for (f = 0; f < numFields; ++f) {
5375: PetscSectionGetFieldDof(section, points[p], f, &fdof);
5376: offsets[f+1] += fdof;
5377: }
5378: numIndices += dof;
5379: }
5380: for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f];
5382: if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices);
5383: DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);
5384: if (numFields) {
5385: for (p = 0; p < numPoints*2; p += 2) {
5386: PetscInt o = points[p+1];
5387: PetscSectionGetOffset(globalSection, points[p], &globalOff);
5388: indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
5389: }
5390: } else {
5391: for (p = 0, off = 0; p < numPoints*2; p += 2) {
5392: PetscInt o = points[p+1];
5393: PetscSectionGetOffset(globalSection, points[p], &globalOff);
5394: indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices);
5395: }
5396: }
5397: if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);}
5398: MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
5399: if (ierr) {
5400: PetscMPIInt rank;
5401: PetscErrorCode ierr2;
5403: ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
5404: ierr2 = (*PetscErrorPrintf)("[%D]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
5405: ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2);
5406: ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2);
5407:
5408: }
5409: if (!clPoints) {
5410: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
5411: } else {
5412: ISRestoreIndices(clPoints, &clp);
5413: }
5414: DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);
5415: return(0);
5416: }
5420: PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode)
5421: {
5422: DM_Plex *mesh = (DM_Plex*) dmf->data;
5423: PetscInt *fpoints = NULL, *ftotpoints = NULL;
5424: PetscInt *cpoints = NULL;
5425: PetscInt *findices, *cindices;
5426: PetscInt foffsets[32], coffsets[32];
5427: CellRefiner cellRefiner;
5428: PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
5429: PetscErrorCode ierr;
5434: if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
5436: if (!csection) {DMGetDefaultSection(dmc, &csection);}
5438: if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
5440: if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
5443: PetscSectionGetNumFields(fsection, &numFields);
5444: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
5445: PetscMemzero(foffsets, 32 * sizeof(PetscInt));
5446: PetscMemzero(coffsets, 32 * sizeof(PetscInt));
5447: /* Column indices */
5448: DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5449: maxFPoints = numCPoints;
5450: /* Compress out points not in the section */
5451: /* TODO: Squeeze out points with 0 dof as well */
5452: PetscSectionGetChart(csection, &pStart, &pEnd);
5453: for (p = 0, q = 0; p < numCPoints*2; p += 2) {
5454: if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
5455: cpoints[q*2] = cpoints[p];
5456: cpoints[q*2+1] = cpoints[p+1];
5457: ++q;
5458: }
5459: }
5460: numCPoints = q;
5461: for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
5462: PetscInt fdof;
5464: PetscSectionGetDof(csection, cpoints[p], &dof);
5465: if (!dof) continue;
5466: for (f = 0; f < numFields; ++f) {
5467: PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
5468: coffsets[f+1] += fdof;
5469: }
5470: numCIndices += dof;
5471: }
5472: for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
5473: /* Row indices */
5474: DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
5475: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
5476: DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
5477: for (r = 0, q = 0; r < numSubcells; ++r) {
5478: /* TODO Map from coarse to fine cells */
5479: DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
5480: /* Compress out points not in the section */
5481: PetscSectionGetChart(fsection, &pStart, &pEnd);
5482: for (p = 0; p < numFPoints*2; p += 2) {
5483: if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
5484: PetscSectionGetDof(fsection, fpoints[p], &dof);
5485: if (!dof) continue;
5486: for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
5487: if (s < q) continue;
5488: ftotpoints[q*2] = fpoints[p];
5489: ftotpoints[q*2+1] = fpoints[p+1];
5490: ++q;
5491: }
5492: }
5493: DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
5494: }
5495: numFPoints = q;
5496: for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
5497: PetscInt fdof;
5499: PetscSectionGetDof(fsection, ftotpoints[p], &dof);
5500: if (!dof) continue;
5501: for (f = 0; f < numFields; ++f) {
5502: PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
5503: foffsets[f+1] += fdof;
5504: }
5505: numFIndices += dof;
5506: }
5507: for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];
5509: if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
5510: if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
5511: DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);
5512: DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
5513: if (numFields) {
5514: for (p = 0; p < numFPoints*2; p += 2) {
5515: PetscInt o = ftotpoints[p+1];
5516: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5517: indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
5518: }
5519: for (p = 0; p < numCPoints*2; p += 2) {
5520: PetscInt o = cpoints[p+1];
5521: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5522: indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5523: }
5524: } else {
5525: for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5526: PetscInt o = ftotpoints[p+1];
5527: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5528: indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5529: }
5530: for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5531: PetscInt o = cpoints[p+1];
5532: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5533: indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5534: }
5535: }
5536: if (mesh->printSetValues) {DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);}
5537: MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode);
5538: if (ierr) {
5539: PetscMPIInt rank;
5540: PetscErrorCode ierr2;
5542: ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2);
5543: ierr2 = (*PetscErrorPrintf)("[%D]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2);
5544: ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2);
5545: ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2);
5546: ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2);
5547:
5548: }
5549: DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5550: DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5551: DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);
5552: DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);
5553: return(0);
5554: }
5558: PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])
5559: {
5560: PetscInt *fpoints = NULL, *ftotpoints = NULL;
5561: PetscInt *cpoints = NULL;
5562: PetscInt foffsets[32], coffsets[32];
5563: CellRefiner cellRefiner;
5564: PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f;
5570: if (!fsection) {DMGetDefaultSection(dmf, &fsection);}
5572: if (!csection) {DMGetDefaultSection(dmc, &csection);}
5574: if (!globalFSection) {DMGetDefaultGlobalSection(dmf, &globalFSection);}
5576: if (!globalCSection) {DMGetDefaultGlobalSection(dmc, &globalCSection);}
5578: PetscSectionGetNumFields(fsection, &numFields);
5579: if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields);
5580: PetscMemzero(foffsets, 32 * sizeof(PetscInt));
5581: PetscMemzero(coffsets, 32 * sizeof(PetscInt));
5582: /* Column indices */
5583: DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5584: maxFPoints = numCPoints;
5585: /* Compress out points not in the section */
5586: /* TODO: Squeeze out points with 0 dof as well */
5587: PetscSectionGetChart(csection, &pStart, &pEnd);
5588: for (p = 0, q = 0; p < numCPoints*2; p += 2) {
5589: if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) {
5590: cpoints[q*2] = cpoints[p];
5591: cpoints[q*2+1] = cpoints[p+1];
5592: ++q;
5593: }
5594: }
5595: numCPoints = q;
5596: for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) {
5597: PetscInt fdof;
5599: PetscSectionGetDof(csection, cpoints[p], &dof);
5600: if (!dof) continue;
5601: for (f = 0; f < numFields; ++f) {
5602: PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);
5603: coffsets[f+1] += fdof;
5604: }
5605: numCIndices += dof;
5606: }
5607: for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f];
5608: /* Row indices */
5609: DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);
5610: CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);
5611: DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);
5612: for (r = 0, q = 0; r < numSubcells; ++r) {
5613: /* TODO Map from coarse to fine cells */
5614: DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);
5615: /* Compress out points not in the section */
5616: PetscSectionGetChart(fsection, &pStart, &pEnd);
5617: for (p = 0; p < numFPoints*2; p += 2) {
5618: if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) {
5619: PetscSectionGetDof(fsection, fpoints[p], &dof);
5620: if (!dof) continue;
5621: for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break;
5622: if (s < q) continue;
5623: ftotpoints[q*2] = fpoints[p];
5624: ftotpoints[q*2+1] = fpoints[p+1];
5625: ++q;
5626: }
5627: }
5628: DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);
5629: }
5630: numFPoints = q;
5631: for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) {
5632: PetscInt fdof;
5634: PetscSectionGetDof(fsection, ftotpoints[p], &dof);
5635: if (!dof) continue;
5636: for (f = 0; f < numFields; ++f) {
5637: PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);
5638: foffsets[f+1] += fdof;
5639: }
5640: numFIndices += dof;
5641: }
5642: for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f];
5644: if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices);
5645: if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices);
5646: if (numFields) {
5647: for (p = 0; p < numFPoints*2; p += 2) {
5648: PetscInt o = ftotpoints[p+1];
5649: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5650: indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices);
5651: }
5652: for (p = 0; p < numCPoints*2; p += 2) {
5653: PetscInt o = cpoints[p+1];
5654: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5655: indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices);
5656: }
5657: } else {
5658: for (p = 0, off = 0; p < numFPoints*2; p += 2) {
5659: PetscInt o = ftotpoints[p+1];
5660: PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);
5661: indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices);
5662: }
5663: for (p = 0, off = 0; p < numCPoints*2; p += 2) {
5664: PetscInt o = cpoints[p+1];
5665: PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);
5666: indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices);
5667: }
5668: }
5669: DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);
5670: DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);
5671: return(0);
5672: }
5676: /*@
5677: DMPlexGetHybridBounds - Get the first mesh point of each dimension which is a hybrid
5679: Input Parameter:
5680: . dm - The DMPlex object
5682: Output Parameters:
5683: + cMax - The first hybrid cell
5684: . fMax - The first hybrid face
5685: . eMax - The first hybrid edge
5686: - vMax - The first hybrid vertex
5688: Level: developer
5690: .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds()
5691: @*/
5692: PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
5693: {
5694: DM_Plex *mesh = (DM_Plex*) dm->data;
5695: PetscInt dim;
5700: DMPlexGetDimension(dm, &dim);
5701: if (cMax) *cMax = mesh->hybridPointMax[dim];
5702: if (fMax) *fMax = mesh->hybridPointMax[dim-1];
5703: if (eMax) *eMax = mesh->hybridPointMax[1];
5704: if (vMax) *vMax = mesh->hybridPointMax[0];
5705: return(0);
5706: }
5710: /*@
5711: DMPlexSetHybridBounds - Set the first mesh point of each dimension which is a hybrid
5713: Input Parameters:
5714: . dm - The DMPlex object
5715: . cMax - The first hybrid cell
5716: . fMax - The first hybrid face
5717: . eMax - The first hybrid edge
5718: - vMax - The first hybrid vertex
5720: Level: developer
5722: .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds()
5723: @*/
5724: PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax)
5725: {
5726: DM_Plex *mesh = (DM_Plex*) dm->data;
5727: PetscInt dim;
5732: DMPlexGetDimension(dm, &dim);
5733: if (cMax >= 0) mesh->hybridPointMax[dim] = cMax;
5734: if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax;
5735: if (eMax >= 0) mesh->hybridPointMax[1] = eMax;
5736: if (vMax >= 0) mesh->hybridPointMax[0] = vMax;
5737: return(0);
5738: }
5742: PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight)
5743: {
5744: DM_Plex *mesh = (DM_Plex*) dm->data;
5749: *cellHeight = mesh->vtkCellHeight;
5750: return(0);
5751: }
5755: PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight)
5756: {
5757: DM_Plex *mesh = (DM_Plex*) dm->data;
5761: mesh->vtkCellHeight = cellHeight;
5762: return(0);
5763: }
5767: /* We can easily have a form that takes an IS instead */
5768: static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
5769: {
5770: PetscSection section, globalSection;
5771: PetscInt *numbers, p;
5775: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
5776: PetscSectionSetChart(section, pStart, pEnd);
5777: for (p = pStart; p < pEnd; ++p) {
5778: PetscSectionSetDof(section, p, 1);
5779: }
5780: PetscSectionSetUp(section);
5781: PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, &globalSection);
5782: PetscMalloc1((pEnd - pStart), &numbers);
5783: for (p = pStart; p < pEnd; ++p) {
5784: PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);
5785: if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
5786: else numbers[p-pStart] += shift;
5787: }
5788: ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);
5789: if (globalSize) {
5790: PetscLayout layout;
5791: PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);
5792: PetscLayoutGetSize(layout, globalSize);
5793: PetscLayoutDestroy(&layout);
5794: }
5795: PetscSectionDestroy(§ion);
5796: PetscSectionDestroy(&globalSection);
5797: return(0);
5798: }
5802: PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers)
5803: {
5804: DM_Plex *mesh = (DM_Plex*) dm->data;
5805: PetscInt cellHeight, cStart, cEnd, cMax;
5810: if (!mesh->globalCellNumbers) {
5811: DMPlexGetVTKCellHeight(dm, &cellHeight);
5812: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
5813: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
5814: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
5815: DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, &mesh->globalCellNumbers);
5816: }
5817: *globalCellNumbers = mesh->globalCellNumbers;
5818: return(0);
5819: }
5823: PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers)
5824: {
5825: DM_Plex *mesh = (DM_Plex*) dm->data;
5826: PetscInt vStart, vEnd, vMax;
5831: if (!mesh->globalVertexNumbers) {
5832: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
5833: DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
5834: if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
5835: DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, &mesh->globalVertexNumbers);
5836: }
5837: *globalVertexNumbers = mesh->globalVertexNumbers;
5838: return(0);
5839: }
5843: PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
5844: {
5845: IS nums[4];
5846: PetscInt depths[4];
5847: PetscInt depth, d, shift = 0;
5852: DMPlexGetDepth(dm, &depth);
5853: depths[0] = depth; depths[1] = 0;
5854: for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
5855: for (d = 0; d <= depth; ++d) {
5856: PetscInt pStart, pEnd, gsize;
5858: DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);
5859: DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);
5860: shift += gsize;
5861: }
5862: ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
5863: for (d = 0; d <= depth; ++d) {ISDestroy(&nums[d]);}
5864: return(0);
5865: }
5870: /*@C
5871: PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using
5872: the local section and an SF describing the section point overlap.
5874: Input Parameters:
5875: + s - The PetscSection for the local field layout
5876: . sf - The SF describing parallel layout of the section points
5877: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
5878: . label - The label specifying the points
5879: - labelValue - The label stratum specifying the points
5881: Output Parameter:
5882: . gsection - The PetscSection for the global field layout
5884: Note: This gives negative sizes and offsets to points not owned by this process
5886: Level: developer
5888: .seealso: PetscSectionCreate()
5889: @*/
5890: PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection)
5891: {
5892: PetscInt *neg = NULL, *tmpOff = NULL;
5893: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
5897: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
5898: PetscSectionGetChart(s, &pStart, &pEnd);
5899: PetscSectionSetChart(*gsection, pStart, pEnd);
5900: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
5901: if (nroots >= 0) {
5902: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
5903: PetscCalloc1(nroots, &neg);
5904: if (nroots > pEnd-pStart) {
5905: PetscCalloc1(nroots, &tmpOff);
5906: } else {
5907: tmpOff = &(*gsection)->atlasDof[-pStart];
5908: }
5909: }
5910: /* Mark ghost points with negative dof */
5911: for (p = pStart; p < pEnd; ++p) {
5912: PetscInt value;
5914: DMLabelGetValue(label, p, &value);
5915: if (value != labelValue) continue;
5916: PetscSectionGetDof(s, p, &dof);
5917: PetscSectionSetDof(*gsection, p, dof);
5918: PetscSectionGetConstraintDof(s, p, &cdof);
5919: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
5920: if (neg) neg[p] = -(dof+1);
5921: }
5922: PetscSectionSetUpBC(*gsection);
5923: if (nroots >= 0) {
5924: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5925: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5926: if (nroots > pEnd-pStart) {
5927: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
5928: }
5929: }
5930: /* Calculate new sizes, get proccess offset, and calculate point offsets */
5931: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5932: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
5933: (*gsection)->atlasOff[p] = off;
5934: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
5935: }
5936: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
5937: globalOff -= off;
5938: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
5939: (*gsection)->atlasOff[p] += globalOff;
5940: if (neg) neg[p] = -((*gsection)->atlasOff[p]+1);
5941: }
5942: /* Put in negative offsets for ghost points */
5943: if (nroots >= 0) {
5944: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
5945: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
5946: if (nroots > pEnd-pStart) {
5947: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
5948: }
5949: }
5950: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
5951: PetscFree(neg);
5952: return(0);
5953: }
5957: /*@
5958: DMPlexCheckSymmetry - Check that the adjacency information in the mesh is symmetric.
5960: Input Parameters:
5961: + dm - The DMPlex object
5963: Note: This is a useful diagnostic when creating meshes programmatically.
5965: Level: developer
5967: .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces()
5968: @*/
5969: PetscErrorCode DMPlexCheckSymmetry(DM dm)
5970: {
5971: PetscSection coneSection, supportSection;
5972: const PetscInt *cone, *support;
5973: PetscInt coneSize, c, supportSize, s;
5974: PetscInt pStart, pEnd, p, csize, ssize;
5975: PetscErrorCode ierr;
5979: DMPlexGetConeSection(dm, &coneSection);
5980: DMPlexGetSupportSection(dm, &supportSection);
5981: /* Check that point p is found in the support of its cone points, and vice versa */
5982: DMPlexGetChart(dm, &pStart, &pEnd);
5983: for (p = pStart; p < pEnd; ++p) {
5984: DMPlexGetConeSize(dm, p, &coneSize);
5985: DMPlexGetCone(dm, p, &cone);
5986: for (c = 0; c < coneSize; ++c) {
5987: PetscBool dup = PETSC_FALSE;
5988: PetscInt d;
5989: for (d = c-1; d >= 0; --d) {
5990: if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;}
5991: }
5992: DMPlexGetSupportSize(dm, cone[c], &supportSize);
5993: DMPlexGetSupport(dm, cone[c], &support);
5994: for (s = 0; s < supportSize; ++s) {
5995: if (support[s] == p) break;
5996: }
5997: if ((s >= supportSize) || (dup && (support[s+1] != p))) {
5998: PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p);
5999: for (s = 0; s < coneSize; ++s) {
6000: PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]);
6001: }
6002: PetscPrintf(PETSC_COMM_SELF, "\n");
6003: PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]);
6004: for (s = 0; s < supportSize; ++s) {
6005: PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]);
6006: }
6007: PetscPrintf(PETSC_COMM_SELF, "\n");
6008: if (dup) {
6009: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]);
6010: } else {
6011: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]);
6012: }
6013: }
6014: }
6015: DMPlexGetSupportSize(dm, p, &supportSize);
6016: DMPlexGetSupport(dm, p, &support);
6017: for (s = 0; s < supportSize; ++s) {
6018: DMPlexGetConeSize(dm, support[s], &coneSize);
6019: DMPlexGetCone(dm, support[s], &cone);
6020: for (c = 0; c < coneSize; ++c) {
6021: if (cone[c] == p) break;
6022: }
6023: if (c >= coneSize) {
6024: PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p);
6025: for (c = 0; c < supportSize; ++c) {
6026: PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]);
6027: }
6028: PetscPrintf(PETSC_COMM_SELF, "\n");
6029: PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]);
6030: for (c = 0; c < coneSize; ++c) {
6031: PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]);
6032: }
6033: PetscPrintf(PETSC_COMM_SELF, "\n");
6034: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]);
6035: }
6036: }
6037: }
6038: PetscSectionGetStorageSize(coneSection, &csize);
6039: PetscSectionGetStorageSize(supportSection, &ssize);
6040: if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize);
6041: return(0);
6042: }
6046: /*@
6047: DMPlexCheckSkeleton - Check that each cell has the correct number of vertices
6049: Input Parameters:
6050: + dm - The DMPlex object
6051: . isSimplex - Are the cells simplices or tensor products
6052: - cellHeight - Normally 0
6054: Note: This is a useful diagnostic when creating meshes programmatically.
6056: Level: developer
6058: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces()
6059: @*/
6060: PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight)
6061: {
6062: PetscInt dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c;
6067: DMPlexGetDimension(dm, &dim);
6068: switch (dim) {
6069: case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break;
6070: case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break;
6071: case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break;
6072: default:
6073: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim);
6074: }
6075: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
6076: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
6077: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
6078: cMax = cMax >= 0 ? cMax : cEnd;
6079: for (c = cStart; c < cMax; ++c) {
6080: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
6082: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6083: for (cl = 0; cl < closureSize*2; cl += 2) {
6084: const PetscInt p = closure[cl];
6085: if ((p >= vStart) && (p < vEnd)) ++coneSize;
6086: }
6087: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6088: if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d vertices != %d", c, coneSize, numCorners);
6089: }
6090: for (c = cMax; c < cEnd; ++c) {
6091: PetscInt *closure = NULL, closureSize, cl, coneSize = 0;
6093: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6094: for (cl = 0; cl < closureSize*2; cl += 2) {
6095: const PetscInt p = closure[cl];
6096: if ((p >= vStart) && (p < vEnd)) ++coneSize;
6097: }
6098: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6099: if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has %d vertices > %d", c, coneSize, numHybridCorners);
6100: }
6101: return(0);
6102: }
6106: /*@
6107: DMPlexCheckFaces - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type
6109: Input Parameters:
6110: + dm - The DMPlex object
6111: . isSimplex - Are the cells simplices or tensor products
6112: - cellHeight - Normally 0
6114: Note: This is a useful diagnostic when creating meshes programmatically.
6116: Level: developer
6118: .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton()
6119: @*/
6120: PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight)
6121: {
6122: PetscInt pMax[4];
6123: PetscInt dim, vStart, vEnd, cStart, cEnd, c, h;
6128: DMPlexGetDimension(dm, &dim);
6129: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
6130: DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);
6131: for (h = cellHeight; h < dim; ++h) {
6132: DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);
6133: for (c = cStart; c < cEnd; ++c) {
6134: const PetscInt *cone, *ornt, *faces;
6135: PetscInt numFaces, faceSize, coneSize,f;
6136: PetscInt *closure = NULL, closureSize, cl, numCorners = 0;
6138: if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue;
6139: DMPlexGetConeSize(dm, c, &coneSize);
6140: DMPlexGetCone(dm, c, &cone);
6141: DMPlexGetConeOrientation(dm, c, &ornt);
6142: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6143: for (cl = 0; cl < closureSize*2; cl += 2) {
6144: const PetscInt p = closure[cl];
6145: if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p;
6146: }
6147: DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);
6148: if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces);
6149: for (f = 0; f < numFaces; ++f) {
6150: PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;
6152: DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);
6153: for (cl = 0; cl < fclosureSize*2; cl += 2) {
6154: const PetscInt p = fclosure[cl];
6155: if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;
6156: }
6157: if (fnumCorners != faceSize) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d (%d) of cell %d has %d vertices but should have %d", cone[f], f, c, fnumCorners, faceSize);
6158: for (v = 0; v < fnumCorners; ++v) {
6159: if (fclosure[v] != faces[f*faceSize+v]) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d (%d) of cell %d vertex %d, %d != %d", cone[f], f, c, v, fclosure[v], faces[f*faceSize+v]);
6160: }
6161: DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);
6162: }
6163: DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);
6164: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
6165: }
6166: }
6167: return(0);
6168: }
6172: /* Pointwise interpolation
6173: Just code FEM for now
6174: u^f = I u^c
6175: sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j
6176: u^f_i = sum_j psi^f_i I phi^c_j u^c_j
6177: I_{ij} = psi^f_i phi^c_j
6178: */
6179: PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
6180: {
6181: PetscSection gsc, gsf;
6182: PetscInt m, n;
6183: void *ctx;
6187: /*
6188: Loop over coarse cells
6189: Loop over coarse basis functions
6190: Loop over fine cells in coarse cell
6191: Loop over fine dual basis functions
6192: Evaluate coarse basis on fine dual basis quad points
6193: Sum
6194: Update local element matrix
6195: Accumulate to interpolation matrix
6197: Can extend PetscFEIntegrateJacobian_Basic() to do a specialized cell loop
6198: */
6199: DMGetDefaultGlobalSection(dmFine, &gsf);
6200: PetscSectionGetConstrainedStorageSize(gsf, &m);
6201: DMGetDefaultGlobalSection(dmCoarse, &gsc);
6202: PetscSectionGetConstrainedStorageSize(gsc, &n);
6203: /* We need to preallocate properly */
6204: MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);
6205: MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
6206: MatSetType(*interpolation, dmCoarse->mattype);
6207: DMGetApplicationContext(dmFine, &ctx);
6208: DMPlexComputeInterpolatorFEM(dmCoarse, dmFine, *interpolation, ctx);
6209: /* Use naive scaling */
6210: DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);
6211: return(0);
6212: }
6216: PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, VecScatter *ctx)
6217: {
6221: DMPlexComputeInjectorFEM(dmCoarse, dmFine, ctx, NULL);
6222: return(0);
6223: }
6227: /* Pointwise interpolation
6228: Just code FEM for now
6229: u^f = I u^c
6230: sum_k u^f_k phi^f_k = I sum_l u^c_l phi^c_l
6231: u^f_i = sum_l int psi^f_i I phi^c_l u^c_l
6232: I_{ij} = int psi^f_i phi^c_j
6233: */
6234: PetscErrorCode DMCreateDefaultSection_Plex(DM dm)
6235: {
6236: PetscSection section;
6237: IS *bcPoints;
6238: PetscInt *bcFields, *numComp, *numDof;
6239: PetscInt depth, dim, numBd, numBC = 0, numFields, bd, bc, f;
6243: /* Handle boundary conditions */
6244: DMPlexGetDepth(dm, &depth);
6245: DMPlexGetDimension(dm, &dim);
6246: DMPlexGetNumBoundary(dm, &numBd);
6247: for (bd = 0; bd < numBd; ++bd) {
6248: PetscBool isEssential;
6249: DMPlexGetBoundary(dm, bd, &isEssential, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
6250: if (isEssential) ++numBC;
6251: }
6252: PetscMalloc2(numBC,&bcFields,numBC,&bcPoints);
6253: for (bd = 0, bc = 0; bd < numBd; ++bd) {
6254: const char *bdLabel;
6255: DMLabel label;
6256: const PetscInt *values;
6257: PetscInt bd2, field, numValues;
6258: PetscBool isEssential, duplicate = PETSC_FALSE;
6260: DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);
6261: if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
6262: DMPlexGetLabel(dm, bdLabel, &label);
6263: /* Only want to do this for FEM, and only once */
6264: for (bd2 = 0; bd2 < bd; ++bd2) {
6265: const char *bdname;
6266: DMPlexGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL);
6267: PetscStrcmp(bdname, bdLabel, &duplicate);
6268: if (duplicate) break;
6269: }
6270: if (!duplicate) {
6271: DMPlexLabelComplete(dm, label);
6272: DMPlexLabelAddCells(dm, label);
6273: }
6274: /* Filter out cells, if you actually want to constraint cells you need to do things by hand right now */
6275: if (isEssential) {
6276: IS tmp;
6277: PetscInt *newidx;
6278: const PetscInt *idx;
6279: PetscInt cStart, cEnd, n, p, newn = 0;
6281: bcFields[bc] = field;
6282: DMPlexGetStratumIS(dm, bdLabel, values[0], &tmp);
6283: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
6284: ISGetLocalSize(tmp, &n);
6285: ISGetIndices(tmp, &idx);
6286: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn;
6287: PetscMalloc1(newn,&newidx);
6288: newn = 0;
6289: for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p];
6290: ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);
6291: ISRestoreIndices(tmp, &idx);
6292: ISDestroy(&tmp);
6293: }
6294: }
6295: /* Handle discretization */
6296: DMGetNumFields(dm, &numFields);
6297: PetscMalloc2(numFields,&numComp,numFields*(dim+1),&numDof);
6298: for (f = 0; f < numFields; ++f) {
6299: PetscFE fe;
6300: const PetscInt *numFieldDof;
6301: PetscInt d;
6303: DMGetField(dm, f, (PetscObject *) &fe);
6304: PetscFEGetNumComponents(fe, &numComp[f]);
6305: PetscFEGetNumDof(fe, &numFieldDof);
6306: for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d];
6307: }
6308: for (f = 0; f < numFields; ++f) {
6309: PetscInt d;
6310: for (d = 1; d < dim; ++d) {
6311: if ((numDof[f*(dim+1)+d] > 0) && (depth < dim)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
6312: }
6313: }
6314: DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, NULL, §ion);
6315: for (f = 0; f < numFields; ++f) {
6316: PetscFE fe;
6317: const char *name;
6319: DMGetField(dm, f, (PetscObject *) &fe);
6320: PetscObjectGetName((PetscObject) fe, &name);
6321: PetscSectionSetFieldName(section, f, name);
6322: }
6323: DMSetDefaultSection(dm, section);
6324: PetscSectionDestroy(§ion);
6325: for (bc = 0; bc < numBC; ++bc) {ISDestroy(&bcPoints[bc]);}
6326: PetscFree2(bcFields,bcPoints);
6327: PetscFree2(numComp,numDof);
6328: return(0);
6329: }
6333: /*@
6334: DMPlexGetCoarseDM - Get the coarse mesh from which this was obtained by refinement
6336: Input Parameter:
6337: . dm - The DMPlex object
6339: Output Parameter:
6340: . cdm - The coarse DM
6342: Level: intermediate
6344: .seealso: DMPlexSetCoarseDM()
6345: @*/
6346: PetscErrorCode DMPlexGetCoarseDM(DM dm, DM *cdm)
6347: {
6351: *cdm = ((DM_Plex *) dm->data)->coarseMesh;
6352: return(0);
6353: }
6357: /*@
6358: DMPlexSetCoarseDM - Set the coarse mesh from which this was obtained by refinement
6360: Input Parameters:
6361: + dm - The DMPlex object
6362: - cdm - The coarse DM
6364: Level: intermediate
6366: .seealso: DMPlexGetCoarseDM()
6367: @*/
6368: PetscErrorCode DMPlexSetCoarseDM(DM dm, DM cdm)
6369: {
6370: DM_Plex *mesh;
6376: mesh = (DM_Plex *) dm->data;
6377: DMDestroy(&mesh->coarseMesh);
6378: mesh->coarseMesh = cdm;
6379: PetscObjectReference((PetscObject) mesh->coarseMesh);
6380: return(0);
6381: }