Actual source code: plexgenerate.c
petsc-3.7.7 2017-09-25
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: DMGetLabel(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) {DMCreateLabel(*dm, labelName); DMGetLabel(*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: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
355: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
356: DMGetLabel(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) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*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: DMGetLabel(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) {DMCreateLabel(*dm, labelName); DMGetLabel(*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: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
627: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
628: DMGetLabel(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) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*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(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
755: MPI_Comm_rank(comm, &rank);
756: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
757: DMGetLabel(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: double *meshCoords;
842: int *cells = out->tetrahedronlist;
844: if (sizeof (PetscReal) == sizeof (double)) {
845: meshCoords = (double *) out->pointlist;
846: }
847: else {
848: PetscInt i;
850: PetscMalloc1(3 * numVertices,&meshCoords);
851: for (i = 0; i < 3 * numVertices; i++) {
852: meshCoords[i] = (PetscReal) out->pointlist[i];
853: }
854: }
856: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
857: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
858: if (sizeof (PetscReal) != sizeof (double)) {
859: PetscFree(meshCoords);
860: }
861: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
862: /* Set labels */
863: for (v = 0; v < numVertices; ++v) {
864: if (out->pointmarkerlist[v]) {
865: if (glabel) {DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);}
866: }
867: }
868: if (interpolate) {
869: PetscInt e;
871: for (e = 0; e < out->numberofedges; e++) {
872: if (out->edgemarkerlist[e]) {
873: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
874: const PetscInt *edges;
875: PetscInt numEdges;
877: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
878: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
879: if (glabel) {DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);}
880: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
881: }
882: }
883: for (f = 0; f < out->numberoftrifaces; f++) {
884: if (out->trifacemarkerlist[f]) {
885: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
886: const PetscInt *faces;
887: PetscInt numFaces;
889: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
890: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
891: if (glabel) {DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);}
892: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
893: }
894: }
895: }
896: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
897: }
899: PLCDestroy(&in);
900: PLCDestroy(&out);
901: return(0);
902: }
906: PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
907: {
908: MPI_Comm comm;
909: const PetscInt dim = 3;
910: const char *labelName = "marker";
911: PLC *in, *out;
912: DMLabel label;
913: PetscInt verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
914: PetscMPIInt rank;
918: PetscObjectGetComm((PetscObject)dm,&comm);
919: PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
920: MPI_Comm_rank(comm, &rank);
921: DMPlexGetDepth(dm, &depth);
922: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
923: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
924: DMGetLabel(dm, labelName, &label);
925: PLCCreate(&in);
926: PLCCreate(&out);
928: in->numberofpoints = vEnd - vStart;
929: if (in->numberofpoints > 0) {
930: PetscSection coordSection;
931: Vec coordinates;
932: PetscScalar *array;
934: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
935: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
936: DMGetCoordinatesLocal(dm, &coordinates);
937: DMGetCoordinateSection(dm, &coordSection);
938: VecGetArray(coordinates, &array);
939: for (v = vStart; v < vEnd; ++v) {
940: const PetscInt idx = v - vStart;
941: PetscInt off, d, m = -1;
943: PetscSectionGetOffset(coordSection, v, &off);
944: for (d = 0; d < dim; ++d) {
945: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
946: }
947: if (label) {DMLabelGetValue(label, v, &m);}
949: in->pointmarkerlist[idx] = (int) m;
950: }
951: VecRestoreArray(coordinates, &array);
952: }
953: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
955: in->numberofcorners = 4;
956: in->numberoftetrahedra = cEnd - cStart;
957: in->tetrahedronvolumelist = maxVolumes;
958: if (in->numberoftetrahedra > 0) {
959: PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
960: for (c = cStart; c < cEnd; ++c) {
961: const PetscInt idx = c - cStart;
962: PetscInt *closure = NULL;
963: PetscInt closureSize;
965: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
966: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
967: for (v = 0; v < 4; ++v) {
968: in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
969: }
970: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
971: }
972: }
973: if (!rank) {
974: TetGenOpts t;
976: TetGenOptsInitialize(&t);
978: t.in = dm; /* Should go away */
979: t.refine = 1;
980: t.varvolume = 1;
981: t.quality = 1;
982: t.edgesout = 1;
983: t.zeroindex = 1;
984: t.quiet = 1;
985: t.verbose = verbose; /* Change this */
987: TetGenCheckOpts(&t);
988: TetGenTetrahedralize(&t, in, out);
989: }
990: {
991: DMLabel rlabel = NULL;
992: const PetscInt numCorners = 4;
993: const PetscInt numCells = out->numberoftetrahedra;
994: const PetscInt numVertices = out->numberofpoints;
995: double *meshCoords;
996: int *cells = out->tetrahedronlist;
997: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
999: if (sizeof (PetscReal) == sizeof (double)) {
1000: meshCoords = (double *) out->pointlist;
1001: }
1002: else {
1003: PetscInt i;
1005: PetscMalloc1(3 * numVertices,&meshCoords);
1006: for (i = 0; i < 3 * numVertices; i++) {
1007: meshCoords[i] = (PetscReal) out->pointlist[i];
1008: }
1009: }
1011: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
1012: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
1013: if (sizeof (PetscReal) != sizeof (double)) {
1014: PetscFree(meshCoords);
1015: }
1016: if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
1017: /* Set labels */
1018: for (v = 0; v < numVertices; ++v) {
1019: if (out->pointmarkerlist[v]) {
1020: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);}
1021: }
1022: }
1023: if (interpolate) {
1024: PetscInt e, f;
1026: for (e = 0; e < out->numberofedges; e++) {
1027: if (out->edgemarkerlist[e]) {
1028: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
1029: const PetscInt *edges;
1030: PetscInt numEdges;
1032: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
1033: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
1034: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);}
1035: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
1036: }
1037: }
1038: for (f = 0; f < out->numberoftrifaces; f++) {
1039: if (out->trifacemarkerlist[f]) {
1040: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
1041: const PetscInt *faces;
1042: PetscInt numFaces;
1044: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1045: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
1046: if (rlabel) {DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);}
1047: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
1048: }
1049: }
1050: }
1051: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
1052: }
1053: PLCDestroy(&in);
1054: PLCDestroy(&out);
1055: return(0);
1056: }
1057: #endif /* PETSC_HAVE_CTETGEN */
1061: /*@C
1062: DMPlexGenerate - Generates a mesh.
1064: Not Collective
1066: Input Parameters:
1067: + boundary - The DMPlex boundary object
1068: . name - The mesh generation package name
1069: - interpolate - Flag to create intermediate mesh elements
1071: Output Parameter:
1072: . mesh - The DMPlex object
1074: Level: intermediate
1076: .keywords: mesh, elements
1077: .seealso: DMPlexCreate(), DMRefine()
1078: @*/
1079: PetscErrorCode DMPlexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
1080: {
1081: PetscInt dim;
1082: char genname[1024];
1083: PetscBool isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
1089: DMGetDimension(boundary, &dim);
1090: PetscOptionsGetString(((PetscObject) boundary)->options,((PetscObject) boundary)->prefix, "-dm_plex_generator", genname, 1024, &flg);
1091: if (flg) name = genname;
1092: if (name) {
1093: PetscStrcmp(name, "triangle", &isTriangle);
1094: PetscStrcmp(name, "tetgen", &isTetgen);
1095: PetscStrcmp(name, "ctetgen", &isCTetgen);
1096: }
1097: switch (dim) {
1098: case 1:
1099: if (!name || isTriangle) {
1100: #if defined(PETSC_HAVE_TRIANGLE)
1101: DMPlexGenerate_Triangle(boundary, interpolate, mesh);
1102: #else
1103: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
1104: #endif
1105: } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1106: break;
1107: case 2:
1108: if (!name || isCTetgen) {
1109: #if defined(PETSC_HAVE_CTETGEN)
1110: DMPlexGenerate_CTetgen(boundary, interpolate, mesh);
1111: #else
1112: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1113: #endif
1114: } else if (isTetgen) {
1115: #if defined(PETSC_HAVE_TETGEN)
1116: DMPlexGenerate_Tetgen(boundary, interpolate, mesh);
1117: #else
1118: SETERRQ(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --download-tetgen.");
1119: #endif
1120: } else SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1121: break;
1122: default:
1123: SETERRQ1(PetscObjectComm((PetscObject)boundary), PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
1124: }
1125: return(0);
1126: }
1130: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined)
1131: {
1132: PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *);
1133: PetscReal refinementLimit;
1134: PetscInt dim, cStart, cEnd;
1135: char genname[1024], *name = NULL;
1136: PetscBool isUniform, isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg, localized;
1137: PetscErrorCode ierr;
1140: DMGetCoordinatesLocalized(dm, &localized);
1141: DMPlexGetRefinementUniform(dm, &isUniform);
1142: if (isUniform) {
1143: CellRefiner cellRefiner;
1145: DMPlexGetCellRefiner_Internal(dm, &cellRefiner);
1146: DMPlexRefineUniform_Internal(dm, cellRefiner, dmRefined);
1147: DMCopyBoundary(dm, *dmRefined);
1148: if (localized) {
1149: DMLocalizeCoordinates(*dmRefined);
1150: }
1151: return(0);
1152: }
1153: DMPlexGetRefinementLimit(dm, &refinementLimit);
1154: DMPlexGetRefinementFunction(dm, &refinementFunc);
1155: if (refinementLimit == 0.0 && !refinementFunc) return(0);
1156: DMGetDimension(dm, &dim);
1157: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1158: PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);
1159: if (flg) name = genname;
1160: if (name) {
1161: PetscStrcmp(name, "triangle", &isTriangle);
1162: PetscStrcmp(name, "tetgen", &isTetgen);
1163: PetscStrcmp(name, "ctetgen", &isCTetgen);
1164: }
1165: switch (dim) {
1166: case 2:
1167: if (!name || isTriangle) {
1168: #if defined(PETSC_HAVE_TRIANGLE)
1169: double *maxVolumes;
1170: PetscInt c;
1172: PetscMalloc1(cEnd - cStart, &maxVolumes);
1173: if (refinementFunc) {
1174: for (c = cStart; c < cEnd; ++c) {
1175: PetscReal vol, centroid[3];
1176: PetscReal maxVol;
1178: DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
1179: (*refinementFunc)(centroid, &maxVol);
1180: maxVolumes[c - cStart] = (double) maxVol;
1181: }
1182: } else {
1183: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1184: }
1185: DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);
1186: PetscFree(maxVolumes);
1187: #else
1188: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
1189: #endif
1190: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
1191: break;
1192: case 3:
1193: if (!name || isCTetgen) {
1194: #if defined(PETSC_HAVE_CTETGEN)
1195: PetscReal *maxVolumes;
1196: PetscInt c;
1198: PetscMalloc1(cEnd - cStart, &maxVolumes);
1199: if (refinementFunc) {
1200: for (c = cStart; c < cEnd; ++c) {
1201: PetscReal vol, centroid[3];
1203: DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
1204: (*refinementFunc)(centroid, &maxVolumes[c-cStart]);
1205: }
1206: } else {
1207: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1208: }
1209: DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);
1210: #else
1211: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
1212: #endif
1213: } else if (isTetgen) {
1214: #if defined(PETSC_HAVE_TETGEN)
1215: double *maxVolumes;
1216: PetscInt c;
1218: PetscMalloc1(cEnd - cStart, &maxVolumes);
1219: if (refinementFunc) {
1220: for (c = cStart; c < cEnd; ++c) {
1221: PetscReal vol, centroid[3];
1223: DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
1224: (*refinementFunc)(centroid, &maxVolumes[c-cStart]);
1225: }
1226: } else {
1227: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
1228: }
1229: DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);
1230: #else
1231: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
1232: #endif
1233: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
1234: break;
1235: default:
1236: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
1237: }
1238: DMCopyBoundary(dm, *dmRefined);
1239: if (localized) {
1240: DMLocalizeCoordinates(*dmRefined);
1241: }
1242: return(0);
1243: }
1247: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[])
1248: {
1249: DM cdm = dm;
1250: PetscInt r;
1251: PetscBool isUniform, localized;
1255: DMPlexGetRefinementUniform(dm, &isUniform);
1256: DMGetCoordinatesLocalized(dm, &localized);
1257: if (!isUniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Non-uniform refinement is incompatible with the hierarchy");
1258: for (r = 0; r < nlevels; ++r) {
1259: CellRefiner cellRefiner;
1261: DMPlexGetCellRefiner_Internal(cdm, &cellRefiner);
1262: DMPlexRefineUniform_Internal(cdm, cellRefiner, &dmRefined[r]);
1263: DMCopyBoundary(cdm, dmRefined[r]);
1264: if (localized) {
1265: DMLocalizeCoordinates(dmRefined[r]);
1266: }
1267: DMSetCoarseDM(dmRefined[r], cdm);
1268: DMPlexSetRegularRefinement(dmRefined[r], PETSC_TRUE);
1269: cdm = dmRefined[r];
1270: }
1271: return(0);
1272: }