Actual source code: tetgenerate.cxx
petsc-3.11.4 2019-09-28
1: #include <petsc/private/dmpleximpl.h>
3: #include <tetgen.h>
5: /* This is to fix the tetrahedron orientation from TetGen */
6: static PetscErrorCode DMPlexInvertCells_Internal(PetscInt dim, PetscInt numCells, PetscInt numCorners, int cells[])
7: {
8: PetscInt bound = numCells*numCorners, coff;
12: for (coff = 0; coff < bound; coff += numCorners) {
13: DMPlexInvertCell(dim, numCorners, &cells[coff]);
14: }
15: return(0);
16: }
18: PETSC_EXTERN PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
19: {
20: MPI_Comm comm;
21: DM_Plex *mesh = (DM_Plex *) boundary->data;
22: const PetscInt dim = 3;
23: const char *labelName = "marker";
24: ::tetgenio in;
25: ::tetgenio out;
26: DMLabel label;
27: PetscInt vStart, vEnd, v, fStart, fEnd, f;
28: PetscMPIInt rank;
32: PetscObjectGetComm((PetscObject)boundary,&comm);
33: MPI_Comm_rank(comm, &rank);
34: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
35: DMGetLabel(boundary, labelName, &label);
37: in.numberofpoints = vEnd - vStart;
38: if (in.numberofpoints > 0) {
39: PetscSection coordSection;
40: Vec coordinates;
41: PetscScalar *array;
43: in.pointlist = new double[in.numberofpoints*dim];
44: in.pointmarkerlist = new int[in.numberofpoints];
46: DMGetCoordinatesLocal(boundary, &coordinates);
47: DMGetCoordinateSection(boundary, &coordSection);
48: VecGetArray(coordinates, &array);
49: for (v = vStart; v < vEnd; ++v) {
50: const PetscInt idx = v - vStart;
51: PetscInt off, d;
53: PetscSectionGetOffset(coordSection, v, &off);
54: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
55: if (label) {
56: PetscInt val;
58: DMLabelGetValue(label, v, &val);
59: in.pointmarkerlist[idx] = (int) val;
60: }
61: }
62: VecRestoreArray(coordinates, &array);
63: }
64: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
66: in.numberoffacets = fEnd - fStart;
67: if (in.numberoffacets > 0) {
68: in.facetlist = new tetgenio::facet[in.numberoffacets];
69: in.facetmarkerlist = new int[in.numberoffacets];
70: for (f = fStart; f < fEnd; ++f) {
71: const PetscInt idx = f - fStart;
72: PetscInt *points = NULL, numPoints, p, numVertices = 0, v;
74: in.facetlist[idx].numberofpolygons = 1;
75: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
76: in.facetlist[idx].numberofholes = 0;
77: in.facetlist[idx].holelist = NULL;
79: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
80: for (p = 0; p < numPoints*2; p += 2) {
81: const PetscInt point = points[p];
82: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
83: }
85: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
86: poly->numberofvertices = numVertices;
87: poly->vertexlist = new int[poly->numberofvertices];
88: for (v = 0; v < numVertices; ++v) {
89: const PetscInt vIdx = points[v] - vStart;
90: poly->vertexlist[v] = vIdx;
91: }
92: if (label) {
93: PetscInt val;
95: DMLabelGetValue(label, f, &val);
96: in.facetmarkerlist[idx] = (int) val;
97: }
98: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
99: }
100: }
101: if (!rank) {
102: char args[32];
104: /* Take away 'Q' for verbose output */
105: PetscStrcpy(args, "pqezQ");
106: if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
107: else {::tetrahedralize(args, &in, &out);}
108: }
109: {
110: DMLabel glabel = NULL;
111: const PetscInt numCorners = 4;
112: const PetscInt numCells = out.numberoftetrahedra;
113: const PetscInt numVertices = out.numberofpoints;
114: const double *meshCoords = out.pointlist;
115: int *cells = out.tetrahedronlist;
117: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
118: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
119: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
120: /* Set labels */
121: for (v = 0; v < numVertices; ++v) {
122: if (out.pointmarkerlist[v]) {
123: if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
124: }
125: }
126: if (interpolate) {
127: #if 0
128: PetscInt e;
130: /* This check is never actually executed for ctetgen (which never returns edgemarkers) and seems to be broken for
131: * tetgen */
132: for (e = 0; e < out.numberofedges; e++) {
133: if (out.edgemarkerlist[e]) {
134: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
135: const PetscInt *edges;
136: PetscInt numEdges;
138: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
139: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
140: if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
141: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
142: }
143: }
144: #endif
145: for (f = 0; f < out.numberoftrifaces; f++) {
146: if (out.trifacemarkerlist[f]) {
147: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
148: const PetscInt *faces;
149: PetscInt numFaces;
151: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
152: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
153: if (glabel) {DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);}
154: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
155: }
156: }
157: }
158: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
159: }
160: return(0);
161: }
163: PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
164: {
165: MPI_Comm comm;
166: const PetscInt dim = 3;
167: const char *labelName = "marker";
168: ::tetgenio in;
169: ::tetgenio out;
170: DMLabel label;
171: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
172: PetscMPIInt rank;
176: PetscObjectGetComm((PetscObject)dm,&comm);
177: MPI_Comm_rank(comm, &rank);
178: DMPlexGetDepth(dm, &depth);
179: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
180: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
181: DMGetLabel(dm, labelName, &label);
183: in.numberofpoints = vEnd - vStart;
184: if (in.numberofpoints > 0) {
185: PetscSection coordSection;
186: Vec coordinates;
187: PetscScalar *array;
189: in.pointlist = new double[in.numberofpoints*dim];
190: in.pointmarkerlist = new int[in.numberofpoints];
192: DMGetCoordinatesLocal(dm, &coordinates);
193: DMGetCoordinateSection(dm, &coordSection);
194: VecGetArray(coordinates, &array);
195: for (v = vStart; v < vEnd; ++v) {
196: const PetscInt idx = v - vStart;
197: PetscInt off, d;
199: PetscSectionGetOffset(coordSection, v, &off);
200: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
201: if (label) {
202: PetscInt val;
203:
204: DMLabelGetValue(label, v, &val);
205: in.pointmarkerlist[idx] = (int) val;
206: }
207: }
208: VecRestoreArray(coordinates, &array);
209: }
210: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
212: in.numberofcorners = 4;
213: in.numberoftetrahedra = cEnd - cStart;
214: in.tetrahedronvolumelist = (double*) maxVolumes;
215: if (in.numberoftetrahedra > 0) {
216: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
217: for (c = cStart; c < cEnd; ++c) {
218: const PetscInt idx = c - cStart;
219: PetscInt *closure = NULL;
220: PetscInt closureSize;
222: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
223: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
224: for (v = 0; v < 4; ++v) {
225: in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
226: }
227: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
228: }
229: }
230: /* TODO: Put in boundary faces with markers */
231: if (!rank) {
232: char args[32];
234: #if 1
235: /* Take away 'Q' for verbose output */
236: PetscStrcpy(args, "qezQra");
237: #else
238: PetscStrcpy(args, "qezraVVVV");
239: #endif
240: ::tetrahedralize(args, &in, &out);
241: }
242: in.tetrahedronvolumelist = NULL;
244: {
245: DMLabel rlabel = NULL;
246: const PetscInt numCorners = 4;
247: const PetscInt numCells = out.numberoftetrahedra;
248: const PetscInt numVertices = out.numberofpoints;
249: const double *meshCoords = out.pointlist;
250: int *cells = out.tetrahedronlist;
252: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
254: DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
255: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
256: if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
257: /* Set labels */
258: for (v = 0; v < numVertices; ++v) {
259: if (out.pointmarkerlist[v]) {
260: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
261: }
262: }
263: if (interpolate) {
264: PetscInt f;
265: #if 0
266: PetscInt e;
268: for (e = 0; e < out.numberofedges; e++) {
269: if (out.edgemarkerlist[e]) {
270: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
271: const PetscInt *edges;
272: PetscInt numEdges;
274: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
275: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
276: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
277: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
278: }
279: }
280: #endif
281: for (f = 0; f < out.numberoftrifaces; f++) {
282: if (out.trifacemarkerlist[f]) {
283: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
284: const PetscInt *faces;
285: PetscInt numFaces;
287: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
288: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
289: if (rlabel) {DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);}
290: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
291: }
292: }
293: }
294: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
295: }
296: return(0);
297: }