Actual source code: plexgenerate.c
petsc-3.6.4 2016-04-12
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
5: PetscErrorCode DMPlexInvertCell_Internal(PetscInt dim, PetscInt numCorners, PetscInt cone[])
6: {
7: int tmpc;
10: if (dim != 3) return(0);
11: switch (numCorners) {
12: case 4:
13: tmpc = cone[0];
14: cone[0] = cone[1];
15: cone[1] = tmpc;
16: break;
17: case 8:
18: tmpc = cone[1];
19: cone[1] = cone[3];
20: cone[3] = tmpc;
21: break;
22: default: break;
23: }
24: return(0);
25: }
29: /*@C
30: DMPlexInvertCell - This flips tetrahedron and hexahedron orientation since Plex stores them internally with outward normals. Other cells are left untouched.
32: Input Parameters:
33: + numCorners - The number of vertices in a cell
34: - cone - The incoming cone
36: Output Parameter:
37: . cone - The inverted cone (in-place)
39: Level: developer
41: .seealso: DMPlexGenerate()
42: @*/
43: PetscErrorCode DMPlexInvertCell(PetscInt dim, PetscInt numCorners, int cone[])
44: {
45: int tmpc;
48: if (dim != 3) return(0);
49: switch (numCorners) {
50: case 4:
51: tmpc = cone[0];
52: cone[0] = cone[1];
53: cone[1] = tmpc;
54: break;
55: case 8:
56: tmpc = cone[1];
57: cone[1] = cone[3];
58: cone[3] = tmpc;
59: break;
60: default: break;
61: }
62: return(0);
63: }
67: /* This is to fix the tetrahedron orientation from TetGen */
68: PETSC_UNUSED static PetscErrorCode DMPlexInvertCells_Internal(PetscInt dim, PetscInt numCells, PetscInt numCorners, int cells[])
69: {
70: PetscInt bound = numCells*numCorners, coff;
74: for (coff = 0; coff < bound; coff += numCorners) {
75: DMPlexInvertCell(dim, numCorners, &cells[coff]);
76: }
77: return(0);
78: }
82: /*@
83: DMPlexTriangleSetOptions - Set the options used for the Triangle mesh generator
85: Not Collective
87: Inputs Parameters:
88: + dm - The DMPlex object
89: - opts - The command line options
91: Level: developer
93: .keywords: mesh, points
94: .seealso: DMPlexTetgenSetOptions(), DMPlexGenerate()
95: @*/
96: PetscErrorCode DMPlexTriangleSetOptions(DM dm, const char *opts)
97: {
98: DM_Plex *mesh = (DM_Plex*) dm->data;
104: PetscFree(mesh->triangleOpts);
105: PetscStrallocpy(opts, &mesh->triangleOpts);
106: return(0);
107: }
111: /*@
112: DMPlexTetgenSetOptions - Set the options used for the Tetgen mesh generator
114: Not Collective
116: Inputs Parameters:
117: + dm - The DMPlex object
118: - opts - The command line options
120: Level: developer
122: .keywords: mesh, points
123: .seealso: DMPlexTriangleSetOptions(), DMPlexGenerate()
124: @*/
125: PetscErrorCode DMPlexTetgenSetOptions(DM dm, const char *opts)
126: {
127: DM_Plex *mesh = (DM_Plex*) dm->data;
133: PetscFree(mesh->tetgenOpts);
134: PetscStrallocpy(opts, &mesh->tetgenOpts);
135: return(0);
136: }
138: #if defined(PETSC_HAVE_TRIANGLE)
139: #include <triangle.h>
143: PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx)
144: {
146: inputCtx->numberofpoints = 0;
147: inputCtx->numberofpointattributes = 0;
148: inputCtx->pointlist = NULL;
149: inputCtx->pointattributelist = NULL;
150: inputCtx->pointmarkerlist = NULL;
151: inputCtx->numberofsegments = 0;
152: inputCtx->segmentlist = NULL;
153: inputCtx->segmentmarkerlist = NULL;
154: inputCtx->numberoftriangleattributes = 0;
155: inputCtx->trianglelist = NULL;
156: inputCtx->numberofholes = 0;
157: inputCtx->holelist = NULL;
158: inputCtx->numberofregions = 0;
159: inputCtx->regionlist = NULL;
160: return(0);
161: }
165: PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
166: {
168: outputCtx->numberofpoints = 0;
169: outputCtx->pointlist = NULL;
170: outputCtx->pointattributelist = NULL;
171: outputCtx->pointmarkerlist = NULL;
172: outputCtx->numberoftriangles = 0;
173: outputCtx->trianglelist = NULL;
174: outputCtx->triangleattributelist = NULL;
175: outputCtx->neighborlist = NULL;
176: outputCtx->segmentlist = NULL;
177: outputCtx->segmentmarkerlist = NULL;
178: outputCtx->numberofedges = 0;
179: outputCtx->edgelist = NULL;
180: outputCtx->edgemarkerlist = NULL;
181: return(0);
182: }
186: PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
187: {
189: free(outputCtx->pointlist);
190: free(outputCtx->pointmarkerlist);
191: free(outputCtx->segmentlist);
192: free(outputCtx->segmentmarkerlist);
193: free(outputCtx->edgelist);
194: free(outputCtx->edgemarkerlist);
195: free(outputCtx->trianglelist);
196: free(outputCtx->neighborlist);
197: return(0);
198: }
202: PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
203: {
204: MPI_Comm comm;
205: DM_Plex *mesh = (DM_Plex *) boundary->data;
206: PetscInt dim = 2;
207: const PetscBool createConvexHull = PETSC_FALSE;
208: const PetscBool constrained = PETSC_FALSE;
209: const char *labelName = "marker";
210: struct triangulateio in;
211: struct triangulateio out;
212: DMLabel label;
213: PetscInt vStart, vEnd, v, eStart, eEnd, e;
214: PetscMPIInt rank;
215: PetscErrorCode ierr;
218: PetscObjectGetComm((PetscObject)boundary,&comm);
219: MPI_Comm_rank(comm, &rank);
220: InitInput_Triangle(&in);
221: InitOutput_Triangle(&out);
222: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
223: DMPlexGetLabel(boundary, labelName, &label);
225: in.numberofpoints = vEnd - vStart;
226: if (in.numberofpoints > 0) {
227: PetscSection coordSection;
228: Vec coordinates;
229: PetscScalar *array;
231: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
232: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
233: DMGetCoordinatesLocal(boundary, &coordinates);
234: DMGetCoordinateSection(boundary, &coordSection);
235: VecGetArray(coordinates, &array);
236: for (v = vStart; v < vEnd; ++v) {
237: const PetscInt idx = v - vStart;
238: PetscInt off, d;
240: PetscSectionGetOffset(coordSection, v, &off);
241: for (d = 0; d < dim; ++d) {
242: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
243: }
244: if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
245: }
246: VecRestoreArray(coordinates, &array);
247: }
248: DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
249: in.numberofsegments = eEnd - eStart;
250: if (in.numberofsegments > 0) {
251: PetscMalloc1(in.numberofsegments*2, &in.segmentlist);
252: PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);
253: for (e = eStart; e < eEnd; ++e) {
254: const PetscInt idx = e - eStart;
255: const PetscInt *cone;
257: DMPlexGetCone(boundary, e, &cone);
259: in.segmentlist[idx*2+0] = cone[0] - vStart;
260: in.segmentlist[idx*2+1] = cone[1] - vStart;
262: if (label) {DMLabelGetValue(label, e, &in.segmentmarkerlist[idx]);}
263: }
264: }
265: #if 0 /* Do not currently support holes */
266: PetscReal *holeCoords;
267: PetscInt h, d;
269: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
270: if (in.numberofholes > 0) {
271: PetscMalloc1(in.numberofholes*dim, &in.holelist);
272: for (h = 0; h < in.numberofholes; ++h) {
273: for (d = 0; d < dim; ++d) {
274: in.holelist[h*dim+d] = holeCoords[h*dim+d];
275: }
276: }
277: }
278: #endif
279: if (!rank) {
280: char args[32];
282: /* Take away 'Q' for verbose output */
283: PetscStrcpy(args, "pqezQ");
284: if (createConvexHull) {PetscStrcat(args, "c");}
285: if (constrained) {PetscStrcpy(args, "zepDQ");}
286: if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);}
287: else {triangulate(args, &in, &out, NULL);}
288: }
289: PetscFree(in.pointlist);
290: PetscFree(in.pointmarkerlist);
291: PetscFree(in.segmentlist);
292: PetscFree(in.segmentmarkerlist);
293: PetscFree(in.holelist);
295: {
296: DMLabel glabel = NULL;
297: const PetscInt numCorners = 3;
298: const PetscInt numCells = out.numberoftriangles;
299: const PetscInt numVertices = out.numberofpoints;
300: const int *cells = out.trianglelist;
301: const double *meshCoords = out.pointlist;
303: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
304: if (label) {DMPlexCreateLabel(*dm, labelName); DMPlexGetLabel(*dm, labelName, &glabel);}
305: /* Set labels */
306: for (v = 0; v < numVertices; ++v) {
307: if (out.pointmarkerlist[v]) {
308: if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
309: }
310: }
311: if (interpolate) {
312: for (e = 0; e < out.numberofedges; e++) {
313: if (out.edgemarkerlist[e]) {
314: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
315: const PetscInt *edges;
316: PetscInt numEdges;
318: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
319: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
320: if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
321: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
322: }
323: }
324: }
325: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
326: }
327: #if 0 /* Do not currently support holes */
328: DMPlexCopyHoles(*dm, boundary);
329: #endif
330: FiniOutput_Triangle(&out);
331: return(0);
332: }
336: PetscErrorCode DMPlexRefine_Triangle(DM dm, double *maxVolumes, DM *dmRefined)
337: {
338: MPI_Comm comm;
339: PetscInt dim = 2;
340: const char *labelName = "marker";
341: struct triangulateio in;
342: struct triangulateio out;
343: DMLabel label;
344: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
345: PetscMPIInt rank;
346: PetscErrorCode ierr;
349: PetscObjectGetComm((PetscObject)dm,&comm);
350: MPI_Comm_rank(comm, &rank);
351: InitInput_Triangle(&in);
352: InitOutput_Triangle(&out);
353: DMPlexGetDepth(dm, &depth);
354: MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
355: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
356: DMPlexGetLabel(dm, labelName, &label);
358: in.numberofpoints = vEnd - vStart;
359: if (in.numberofpoints > 0) {
360: PetscSection coordSection;
361: Vec coordinates;
362: PetscScalar *array;
364: PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
365: PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
366: DMGetCoordinatesLocal(dm, &coordinates);
367: DMGetCoordinateSection(dm, &coordSection);
368: VecGetArray(coordinates, &array);
369: for (v = vStart; v < vEnd; ++v) {
370: const PetscInt idx = v - vStart;
371: PetscInt off, d;
373: PetscSectionGetOffset(coordSection, v, &off);
374: for (d = 0; d < dim; ++d) {
375: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
376: }
377: if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
378: }
379: VecRestoreArray(coordinates, &array);
380: }
381: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
383: in.numberofcorners = 3;
384: in.numberoftriangles = cEnd - cStart;
386: in.trianglearealist = (double*) maxVolumes;
387: if (in.numberoftriangles > 0) {
388: PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
389: for (c = cStart; c < cEnd; ++c) {
390: const PetscInt idx = c - cStart;
391: PetscInt *closure = NULL;
392: PetscInt closureSize;
394: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
395: if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
396: for (v = 0; v < 3; ++v) {
397: in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
398: }
399: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
400: }
401: }
402: /* TODO: Segment markers are missing on input */
403: #if 0 /* Do not currently support holes */
404: PetscReal *holeCoords;
405: PetscInt h, d;
407: DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
408: if (in.numberofholes > 0) {
409: PetscMalloc1(in.numberofholes*dim, &in.holelist);
410: for (h = 0; h < in.numberofholes; ++h) {
411: for (d = 0; d < dim; ++d) {
412: in.holelist[h*dim+d] = holeCoords[h*dim+d];
413: }
414: }
415: }
416: #endif
417: if (!rank) {
418: char args[32];
420: /* Take away 'Q' for verbose output */
421: PetscStrcpy(args, "pqezQra");
422: triangulate(args, &in, &out, NULL);
423: }
424: PetscFree(in.pointlist);
425: PetscFree(in.pointmarkerlist);
426: PetscFree(in.segmentlist);
427: PetscFree(in.segmentmarkerlist);
428: PetscFree(in.trianglelist);
430: {
431: DMLabel rlabel = NULL;
432: const PetscInt numCorners = 3;
433: const PetscInt numCells = out.numberoftriangles;
434: const PetscInt numVertices = out.numberofpoints;
435: const int *cells = out.trianglelist;
436: const double *meshCoords = out.pointlist;
437: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
439: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
440: if (label) {DMPlexCreateLabel(*dmRefined, labelName); DMPlexGetLabel(*dmRefined, labelName, &rlabel);}
441: /* Set labels */
442: for (v = 0; v < numVertices; ++v) {
443: if (out.pointmarkerlist[v]) {
444: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
445: }
446: }
447: if (interpolate) {
448: PetscInt e;
450: for (e = 0; e < out.numberofedges; e++) {
451: if (out.edgemarkerlist[e]) {
452: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
453: const PetscInt *edges;
454: PetscInt numEdges;
456: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
457: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
458: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
459: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
460: }
461: }
462: }
463: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
464: }
465: #if 0 /* Do not currently support holes */
466: DMPlexCopyHoles(*dm, boundary);
467: #endif
468: FiniOutput_Triangle(&out);
469: return(0);
470: }
471: #endif /* PETSC_HAVE_TRIANGLE */
473: #if defined(PETSC_HAVE_TETGEN)
474: #include <tetgen.h>
477: PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
478: {
479: MPI_Comm comm;
480: DM_Plex *mesh = (DM_Plex *) boundary->data;
481: const PetscInt dim = 3;
482: const char *labelName = "marker";
483: ::tetgenio in;
484: ::tetgenio out;
485: DMLabel label;
486: PetscInt vStart, vEnd, v, fStart, fEnd, f;
487: PetscMPIInt rank;
491: PetscObjectGetComm((PetscObject)boundary,&comm);
492: MPI_Comm_rank(comm, &rank);
493: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
494: DMPlexGetLabel(boundary, labelName, &label);
496: in.numberofpoints = vEnd - vStart;
497: if (in.numberofpoints > 0) {
498: PetscSection coordSection;
499: Vec coordinates;
500: PetscScalar *array;
502: in.pointlist = new double[in.numberofpoints*dim];
503: in.pointmarkerlist = new int[in.numberofpoints];
505: DMGetCoordinatesLocal(boundary, &coordinates);
506: DMGetCoordinateSection(boundary, &coordSection);
507: VecGetArray(coordinates, &array);
508: for (v = vStart; v < vEnd; ++v) {
509: const PetscInt idx = v - vStart;
510: PetscInt off, d;
512: PetscSectionGetOffset(coordSection, v, &off);
513: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
514: if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
515: }
516: VecRestoreArray(coordinates, &array);
517: }
518: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
520: in.numberoffacets = fEnd - fStart;
521: if (in.numberoffacets > 0) {
522: in.facetlist = new tetgenio::facet[in.numberoffacets];
523: in.facetmarkerlist = new int[in.numberoffacets];
524: for (f = fStart; f < fEnd; ++f) {
525: const PetscInt idx = f - fStart;
526: PetscInt *points = NULL, numPoints, p, numVertices = 0, v;
528: in.facetlist[idx].numberofpolygons = 1;
529: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
530: in.facetlist[idx].numberofholes = 0;
531: in.facetlist[idx].holelist = NULL;
533: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
534: for (p = 0; p < numPoints*2; p += 2) {
535: const PetscInt point = points[p];
536: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
537: }
539: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
540: poly->numberofvertices = numVertices;
541: poly->vertexlist = new int[poly->numberofvertices];
542: for (v = 0; v < numVertices; ++v) {
543: const PetscInt vIdx = points[v] - vStart;
544: poly->vertexlist[v] = vIdx;
545: }
546: if (label) {DMLabelGetValue(label, f, &in.facetmarkerlist[idx]);}
547: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
548: }
549: }
550: if (!rank) {
551: char args[32];
553: /* Take away 'Q' for verbose output */
554: PetscStrcpy(args, "pqezQ");
555: if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
556: else {::tetrahedralize(args, &in, &out);}
557: }
558: {
559: DMLabel glabel = NULL;
560: const PetscInt numCorners = 4;
561: const PetscInt numCells = out.numberoftetrahedra;
562: const PetscInt numVertices = out.numberofpoints;
563: const double *meshCoords = out.pointlist;
564: int *cells = out.tetrahedronlist;
566: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
567: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
568: if (label) {DMPlexCreateLabel(*dm, labelName); DMPlexGetLabel(*dm, labelName, &glabel);}
569: /* Set labels */
570: for (v = 0; v < numVertices; ++v) {
571: if (out.pointmarkerlist[v]) {
572: if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
573: }
574: }
575: if (interpolate) {
576: PetscInt e;
578: for (e = 0; e < out.numberofedges; e++) {
579: if (out.edgemarkerlist[e]) {
580: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
581: const PetscInt *edges;
582: PetscInt numEdges;
584: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
585: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
586: if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
587: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
588: }
589: }
590: for (f = 0; f < out.numberoftrifaces; f++) {
591: if (out.trifacemarkerlist[f]) {
592: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
593: const PetscInt *faces;
594: PetscInt numFaces;
596: DMPlexGetJoin(*dm, 3, vertices, &numFaces, &faces);
597: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
598: if (glabel) {DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);}
599: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
600: }
601: }
602: }
603: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
604: }
605: return(0);
606: }
610: PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
611: {
612: MPI_Comm comm;
613: const PetscInt dim = 3;
614: const char *labelName = "marker";
615: ::tetgenio in;
616: ::tetgenio out;
617: DMLabel label;
618: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
619: PetscMPIInt rank;
623: PetscObjectGetComm((PetscObject)dm,&comm);
624: MPI_Comm_rank(comm, &rank);
625: DMPlexGetDepth(dm, &depth);
626: MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
627: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
628: DMPlexGetLabel(dm, labelName, &label);
630: in.numberofpoints = vEnd - vStart;
631: if (in.numberofpoints > 0) {
632: PetscSection coordSection;
633: Vec coordinates;
634: PetscScalar *array;
636: in.pointlist = new double[in.numberofpoints*dim];
637: in.pointmarkerlist = new int[in.numberofpoints];
639: DMGetCoordinatesLocal(dm, &coordinates);
640: DMGetCoordinateSection(dm, &coordSection);
641: VecGetArray(coordinates, &array);
642: for (v = vStart; v < vEnd; ++v) {
643: const PetscInt idx = v - vStart;
644: PetscInt off, d;
646: PetscSectionGetOffset(coordSection, v, &off);
647: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
648: if (label) {DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);}
649: }
650: VecRestoreArray(coordinates, &array);
651: }
652: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
654: in.numberofcorners = 4;
655: in.numberoftetrahedra = cEnd - cStart;
656: in.tetrahedronvolumelist = (double*) maxVolumes;
657: if (in.numberoftetrahedra > 0) {
658: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
659: for (c = cStart; c < cEnd; ++c) {
660: const PetscInt idx = c - cStart;
661: PetscInt *closure = NULL;
662: PetscInt closureSize;
664: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
665: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
666: for (v = 0; v < 4; ++v) {
667: in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
668: }
669: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
670: }
671: }
672: /* TODO: Put in boundary faces with markers */
673: if (!rank) {
674: char args[32];
676: /* Take away 'Q' for verbose output */
677: /*PetscStrcpy(args, "qezQra"); */
678: PetscStrcpy(args, "qezraVVVV");
679: ::tetrahedralize(args, &in, &out);
680: }
681: in.tetrahedronvolumelist = NULL;
683: {
684: DMLabel rlabel = NULL;
685: const PetscInt numCorners = 4;
686: const PetscInt numCells = out.numberoftetrahedra;
687: const PetscInt numVertices = out.numberofpoints;
688: const double *meshCoords = out.pointlist;
689: int *cells = out.tetrahedronlist;
691: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
693: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
694: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
695: if (label) {DMPlexCreateLabel(*dmRefined, labelName); DMPlexGetLabel(*dmRefined, labelName, &rlabel);}
696: /* Set labels */
697: for (v = 0; v < numVertices; ++v) {
698: if (out.pointmarkerlist[v]) {
699: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
700: }
701: }
702: if (interpolate) {
703: PetscInt e, f;
705: for (e = 0; e < out.numberofedges; e++) {
706: if (out.edgemarkerlist[e]) {
707: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
708: const PetscInt *edges;
709: PetscInt numEdges;
711: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
712: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
713: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
714: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
715: }
716: }
717: for (f = 0; f < out.numberoftrifaces; f++) {
718: if (out.trifacemarkerlist[f]) {
719: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
720: const PetscInt *faces;
721: PetscInt numFaces;
723: DMPlexGetJoin(*dmRefined, 3, vertices, &numFaces, &faces);
724: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
725: if (rlabel) {DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);}
726: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
727: }
728: }
729: }
730: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
731: }
732: return(0);
733: }
734: #endif /* PETSC_HAVE_TETGEN */
736: #if defined(PETSC_HAVE_CTETGEN)
737: #include <ctetgen.h>
741: PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
742: {
743: MPI_Comm comm;
744: const PetscInt dim = 3;
745: const char *labelName = "marker";
746: PLC *in, *out;
747: DMLabel label;
748: PetscInt verbose = 0, vStart, vEnd, v, fStart, fEnd, f;
749: PetscMPIInt rank;
753: PetscObjectGetComm((PetscObject)boundary,&comm);
754: PetscOptionsGetInt(((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
755: MPI_Comm_rank(comm, &rank);
756: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
757: DMPlexGetLabel(boundary, labelName, &label);
758: PLCCreate(&in);
759: PLCCreate(&out);
761: in->numberofpoints = vEnd - vStart;
762: if (in->numberofpoints > 0) {
763: PetscSection coordSection;
764: Vec coordinates;
765: PetscScalar *array;
767: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
768: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
769: DMGetCoordinatesLocal(boundary, &coordinates);
770: DMGetCoordinateSection(boundary, &coordSection);
771: VecGetArray(coordinates, &array);
772: for (v = vStart; v < vEnd; ++v) {
773: const PetscInt idx = v - vStart;
774: PetscInt off, d, m = -1;
776: PetscSectionGetOffset(coordSection, v, &off);
777: for (d = 0; d < dim; ++d) {
778: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
779: }
780: if (label) {DMLabelGetValue(label, v, &m);}
782: in->pointmarkerlist[idx] = (int) m;
783: }
784: VecRestoreArray(coordinates, &array);
785: }
786: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
788: in->numberoffacets = fEnd - fStart;
789: if (in->numberoffacets > 0) {
790: PetscMalloc1(in->numberoffacets, &in->facetlist);
791: PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);
792: for (f = fStart; f < fEnd; ++f) {
793: const PetscInt idx = f - fStart;
794: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
795: polygon *poly;
797: in->facetlist[idx].numberofpolygons = 1;
799: PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);
801: in->facetlist[idx].numberofholes = 0;
802: in->facetlist[idx].holelist = NULL;
804: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
805: for (p = 0; p < numPoints*2; p += 2) {
806: const PetscInt point = points[p];
807: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
808: }
810: poly = in->facetlist[idx].polygonlist;
811: poly->numberofvertices = numVertices;
812: PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
813: for (v = 0; v < numVertices; ++v) {
814: const PetscInt vIdx = points[v] - vStart;
815: poly->vertexlist[v] = vIdx;
816: }
817: if (label) {DMLabelGetValue(label, f, &m);}
818: in->facetmarkerlist[idx] = (int) m;
819: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
820: }
821: }
822: if (!rank) {
823: TetGenOpts t;
825: TetGenOptsInitialize(&t);
826: t.in = boundary; /* Should go away */
827: t.plc = 1;
828: t.quality = 1;
829: t.edgesout = 1;
830: t.zeroindex = 1;
831: t.quiet = 1;
832: t.verbose = verbose;
833: TetGenCheckOpts(&t);
834: TetGenTetrahedralize(&t, in, out);
835: }
836: {
837: DMLabel glabel = NULL;
838: const PetscInt numCorners = 4;
839: const PetscInt numCells = out->numberoftetrahedra;
840: const PetscInt numVertices = out->numberofpoints;
841: const double *meshCoords = out->pointlist;
842: int *cells = out->tetrahedronlist;
844: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
845: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
846: if (label) {DMPlexCreateLabel(*dm, labelName); DMPlexGetLabel(*dm, labelName, &glabel);}
847: /* Set labels */
848: for (v = 0; v < numVertices; ++v) {
849: if (out->pointmarkerlist[v]) {
850: if (glabel) {DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);}
851: }
852: }
853: if (interpolate) {
854: PetscInt e;
856: for (e = 0; e < out->numberofedges; e++) {
857: if (out->edgemarkerlist[e]) {
858: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
859: const PetscInt *edges;
860: PetscInt numEdges;
862: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
863: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
864: if (glabel) {DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);}
865: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
866: }
867: }
868: for (f = 0; f < out->numberoftrifaces; f++) {
869: if (out->trifacemarkerlist[f]) {
870: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
871: const PetscInt *faces;
872: PetscInt numFaces;
874: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
875: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
876: if (glabel) {DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);}
877: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
878: }
879: }
880: }
881: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
882: }
884: PLCDestroy(&in);
885: PLCDestroy(&out);
886: return(0);
887: }
891: PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
892: {
893: MPI_Comm comm;
894: const PetscInt dim = 3;
895: const char *labelName = "marker";
896: PLC *in, *out;
897: DMLabel label;
898: PetscInt verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
899: PetscMPIInt rank;
903: PetscObjectGetComm((PetscObject)dm,&comm);
904: PetscOptionsGetInt(((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
905: MPI_Comm_rank(comm, &rank);
906: DMPlexGetDepth(dm, &depth);
907: MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
908: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
909: DMPlexGetLabel(dm, labelName, &label);
910: PLCCreate(&in);
911: PLCCreate(&out);
913: in->numberofpoints = vEnd - vStart;
914: if (in->numberofpoints > 0) {
915: PetscSection coordSection;
916: Vec coordinates;
917: PetscScalar *array;
919: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
920: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
921: DMGetCoordinatesLocal(dm, &coordinates);
922: DMGetCoordinateSection(dm, &coordSection);
923: VecGetArray(coordinates, &array);
924: for (v = vStart; v < vEnd; ++v) {
925: const PetscInt idx = v - vStart;
926: PetscInt off, d, m = -1;
928: PetscSectionGetOffset(coordSection, v, &off);
929: for (d = 0; d < dim; ++d) {
930: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
931: }
932: if (label) {DMLabelGetValue(label, v, &m);}
934: in->pointmarkerlist[idx] = (int) m;
935: }
936: VecRestoreArray(coordinates, &array);
937: }
938: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
940: in->numberofcorners = 4;
941: in->numberoftetrahedra = cEnd - cStart;
942: in->tetrahedronvolumelist = maxVolumes;
943: if (in->numberoftetrahedra > 0) {
944: PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
945: for (c = cStart; c < cEnd; ++c) {
946: const PetscInt idx = c - cStart;
947: PetscInt *closure = NULL;
948: PetscInt closureSize;
950: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
951: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
952: for (v = 0; v < 4; ++v) {
953: in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
954: }
955: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
956: }
957: }
958: if (!rank) {
959: TetGenOpts t;
961: TetGenOptsInitialize(&t);
963: t.in = dm; /* Should go away */
964: t.refine = 1;
965: t.varvolume = 1;
966: t.quality = 1;
967: t.edgesout = 1;
968: t.zeroindex = 1;
969: t.quiet = 1;
970: t.verbose = verbose; /* Change this */
972: TetGenCheckOpts(&t);
973: TetGenTetrahedralize(&t, in, out);
974: }
975: {
976: DMLabel rlabel = NULL;
977: const PetscInt numCorners = 4;
978: const PetscInt numCells = out->numberoftetrahedra;
979: const PetscInt numVertices = out->numberofpoints;
980: const double *meshCoords = out->pointlist;
981: int *cells = out->tetrahedronlist;
982: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
984: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
985: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
986: if (label) {DMPlexCreateLabel(*dmRefined, labelName); DMPlexGetLabel(*dmRefined, labelName, &rlabel);}
987: /* Set labels */
988: for (v = 0; v < numVertices; ++v) {
989: if (out->pointmarkerlist[v]) {
990: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);}
991: }
992: }
993: if (interpolate) {
994: PetscInt e, f;
996: for (e = 0; e < out->numberofedges; e++) {
997: if (out->edgemarkerlist[e]) {
998: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
999: const PetscInt *edges;
1000: PetscInt numEdges;
1002: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
1003: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
1004: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);}
1005: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
1006: }
1007: }
1008: for (f = 0; f < out->numberoftrifaces; f++) {
1009: if (out->trifacemarkerlist[f]) {
1010: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
1011: const PetscInt *faces;
1012: PetscInt numFaces;
1014: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1015: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
1016: if (rlabel) {DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);}
1017: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1018: }
1019: }
1020: }
1021: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
1022: }
1023: PLCDestroy(&in);
1024: PLCDestroy(&out);
1025: return(0);
1026: }
1027: #endif /* PETSC_HAVE_CTETGEN */
1031: /*@C
1032: DMPlexGenerate - Generates a mesh.
1034: Not Collective
1036: Input Parameters:
1037: + boundary - The DMPlex boundary object
1038: . name - The mesh generation package name
1039: - interpolate - Flag to create intermediate mesh elements
1041: Output Parameter:
1042: . mesh - The DMPlex object
1044: Level: intermediate
1046: .keywords: mesh, elements
1047: .seealso: DMPlexCreate(), DMRefine()
1048: @*/
1049: PetscErrorCode DMPlexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
1050: {
1051: PetscInt dim;
1052: char genname[1024];
1053: PetscBool isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1059: DMGetDimension(boundary, &dim);
1060: PetscOptionsGetString(((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);
1061: if (flg) name = genname;
1062: if (name) {
1063: PetscStrcmp(name, "triangle", &isTriangle);
1064: PetscStrcmp(name, "tetgen", &isTetgen);
1065: PetscStrcmp(name, "ctetgen", &isCTetgen);
1066: }
1067: switch (dim) {
1068: case 1:
1069: if (!name || isTriangle) {
1070: #if defined(PETSC_HAVE_TRIANGLE)
1071: DMPlexGenerate_Triangle(boundary, interpolate, mesh);
1072: #else
1073: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
1074: #endif
1075: } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1076: break;
1077: case 2:
1078: if (!name || isCTetgen) {
1079: #if defined(PETSC_HAVE_CTETGEN)
1080: DMPlexGenerate_CTetgen(boundary, interpolate, mesh);
1081: #else
1082: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1083: #endif
1084: } else if (isTetgen) {
1085: #if defined(PETSC_HAVE_TETGEN)
1086: DMPlexGenerate_Tetgen(boundary, interpolate, mesh);
1087: #else
1088: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
1089: #endif
1090: } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1091: break;
1092: default:
1093: SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
1094: }
1095: return(0);
1096: }
1100: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
1101: {
1102: PetscReal refinementLimit;
1103: PetscInt dim, cStart, cEnd;
1104: char genname[1024], *name = NULL;
1105: PetscBool isUniform, isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1109: DMPlexGetRefinementUniform(dm, &isUniform);
1110: if (isUniform) {
1111: CellRefiner cellRefiner;
1113: DMPlexGetCellRefiner_Internal(dm, &cellRefiner);
1114: DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);
1115: DMPlexCopyBoundary(dm, *dmRefined);
1116: return(0);
1117: }
1118: DMPlexGetRefinementLimit(dm, &refinementLimit);
1119: if (refinementLimit == 0.0) return(0);
1120: DMGetDimension(dm, &dim);
1121: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1122: PetscOptionsGetString(((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);
1123: if (flg) name = genname;
1124: if (name) {
1125: PetscStrcmp(name, "triangle", &isTriangle);
1126: PetscStrcmp(name, "tetgen", &isTetgen);
1127: PetscStrcmp(name, "ctetgen", &isCTetgen);
1128: }
1129: switch (dim) {
1130: case 2:
1131: if (!name || isTriangle) {
1132: #if defined(PETSC_HAVE_TRIANGLE)
1133: double *maxVolumes;
1134: PetscInt c;
1136: PetscMalloc1(cEnd - cStart, &maxVolumes);
1137: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1138: DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);
1139: PetscFree(maxVolumes);
1140: #else
1141: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
1142: #endif
1143: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1144: break;
1145: case 3:
1146: if (!name || isCTetgen) {
1147: #if defined(PETSC_HAVE_CTETGEN)
1148: PetscReal *maxVolumes;
1149: PetscInt c;
1151: PetscMalloc1(cEnd - cStart, &maxVolumes);
1152: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1153: DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);
1154: #else
1155: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1156: #endif
1157: } else if (isTetgen) {
1158: #if defined(PETSC_HAVE_TETGEN)
1159: double *maxVolumes;
1160: PetscInt c;
1162: PetscMalloc1(cEnd - cStart, &maxVolumes);
1163: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1164: DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);
1165: #else
1166: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
1167: #endif
1168: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1169: break;
1170: default:
1171: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
1172: }
1173: DMPlexCopyBoundary(dm, *dmRefined);
1174: return(0);
1175: }
1179: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
1180: {
1181: DM cdm = dm;
1182: PetscInt r;
1183: PetscBool isUniform;
1187: DMPlexGetRefinementUniform(dm, &isUniform);
1188: if (!isUniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Non-uniform refinement is incompatible with the hierarchy");
1189: for (r = 0; r < nlevels; ++r) {
1190: CellRefiner cellRefiner;
1192: DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);
1193: DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);
1194: DMPlexCopyBoundary(cdm, dmRefined[r]);
1195: DMPlexSetCoarseDM(dmRefined[r], cdm);
1196: cdm = dmRefined[r];
1197: }
1198: return(0);
1199: }
1203: PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened)
1204: {
1205: DM_Plex *mesh = (DM_Plex*) dm->data;
1209: PetscObjectReference((PetscObject) mesh->coarseMesh);
1210: *dmCoarsened = mesh->coarseMesh;
1211: return(0);
1212: }