Actual source code: tetgenerate.cxx
petsc-3.14.6 2021-03-30
1: #include <petsc/private/dmpleximpl.h>
3: #if defined(PETSC_HAVE_TETGEN_TETLIBRARY_NEEDED)
4: #define TETLIBRARY
5: #endif
6: #include <tetgen.h>
8: /* This is to fix the tetrahedron orientation from TetGen */
9: static PetscErrorCode DMPlexInvertCells_Tetgen(PetscInt numCells, PetscInt numCorners, PetscInt cells[])
10: {
11: PetscInt bound = numCells*numCorners, coff;
14: #define SWAP(a,b) do { PetscInt tmp = (a); (a) = (b); (b) = tmp; } while (0)
15: for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff],cells[coff+1]);
16: #undef SWAP
17: return(0);
18: }
20: PETSC_EXTERN PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
21: {
22: MPI_Comm comm;
23: DM_Plex *mesh = (DM_Plex *) boundary->data;
24: const PetscInt dim = 3;
25: const char *labelName = "marker";
26: ::tetgenio in;
27: ::tetgenio out;
28: DMLabel label;
29: PetscInt vStart, vEnd, v, fStart, fEnd, f;
30: PetscMPIInt rank;
34: PetscObjectGetComm((PetscObject)boundary,&comm);
35: MPI_Comm_rank(comm, &rank);
36: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
37: DMGetLabel(boundary, labelName, &label);
39: in.numberofpoints = vEnd - vStart;
40: if (in.numberofpoints > 0) {
41: PetscSection coordSection;
42: Vec coordinates;
43: PetscScalar *array;
45: in.pointlist = new double[in.numberofpoints*dim];
46: in.pointmarkerlist = new int[in.numberofpoints];
48: DMGetCoordinatesLocal(boundary, &coordinates);
49: DMGetCoordinateSection(boundary, &coordSection);
50: VecGetArray(coordinates, &array);
51: for (v = vStart; v < vEnd; ++v) {
52: const PetscInt idx = v - vStart;
53: PetscInt off, d;
55: PetscSectionGetOffset(coordSection, v, &off);
56: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
57: if (label) {
58: PetscInt val;
60: DMLabelGetValue(label, v, &val);
61: in.pointmarkerlist[idx] = (int) val;
62: }
63: }
64: VecRestoreArray(coordinates, &array);
65: }
66: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
68: in.numberoffacets = fEnd - fStart;
69: if (in.numberoffacets > 0) {
70: in.facetlist = new tetgenio::facet[in.numberoffacets];
71: in.facetmarkerlist = new int[in.numberoffacets];
72: for (f = fStart; f < fEnd; ++f) {
73: const PetscInt idx = f - fStart;
74: PetscInt *points = NULL, numPoints, p, numVertices = 0, v;
76: in.facetlist[idx].numberofpolygons = 1;
77: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
78: in.facetlist[idx].numberofholes = 0;
79: in.facetlist[idx].holelist = NULL;
81: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
82: for (p = 0; p < numPoints*2; p += 2) {
83: const PetscInt point = points[p];
84: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
85: }
87: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
88: poly->numberofvertices = numVertices;
89: poly->vertexlist = new int[poly->numberofvertices];
90: for (v = 0; v < numVertices; ++v) {
91: const PetscInt vIdx = points[v] - vStart;
92: poly->vertexlist[v] = vIdx;
93: }
94: if (label) {
95: PetscInt val;
97: DMLabelGetValue(label, f, &val);
98: in.facetmarkerlist[idx] = (int) val;
99: }
100: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
101: }
102: }
103: if (!rank) {
104: char args[32];
106: /* Take away 'Q' for verbose output */
107: PetscStrcpy(args, "pqezQ");
108: if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
109: else {::tetrahedralize(args, &in, &out);}
110: }
111: {
112: DMLabel glabel = NULL;
113: const PetscInt numCorners = 4;
114: const PetscInt numCells = out.numberoftetrahedra;
115: const PetscInt numVertices = out.numberofpoints;
116: PetscReal *meshCoords = NULL;
117: PetscInt *cells = NULL;
119: if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
120: meshCoords = (PetscReal *) out.pointlist;
121: } else {
122: PetscInt i;
124: meshCoords = new PetscReal[dim * numVertices];
125: for (i = 0; i < dim * numVertices; i++) {
126: meshCoords[i] = (PetscReal) out.pointlist[i];
127: }
128: }
129: if (sizeof (PetscInt) == sizeof (out.tetrahedronlist[0])) {
130: cells = (PetscInt *) out.tetrahedronlist;
131: } else {
132: PetscInt i;
134: cells = new PetscInt[numCells * numCorners];
135: for (i = 0; i < numCells * numCorners; i++) {
136: cells[i] = (PetscInt) out.tetrahedronlist[i];
137: }
138: }
140: DMPlexInvertCells_Tetgen(numCells, numCorners, cells);
141: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
142: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
143: /* Set labels */
144: for (v = 0; v < numVertices; ++v) {
145: if (out.pointmarkerlist[v]) {
146: if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
147: }
148: }
149: if (interpolate) {
150: #if 0
151: PetscInt e;
153: /* This check is never actually executed for ctetgen (which never returns edgemarkers) and seems to be broken for
154: * tetgen */
155: for (e = 0; e < out.numberofedges; e++) {
156: if (out.edgemarkerlist[e]) {
157: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
158: const PetscInt *edges;
159: PetscInt numEdges;
161: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
162: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
163: if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
164: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
165: }
166: }
167: #endif
168: for (f = 0; f < out.numberoftrifaces; f++) {
169: if (out.trifacemarkerlist[f]) {
170: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
171: const PetscInt *faces;
172: PetscInt numFaces;
174: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
175: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
176: if (glabel) {DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);}
177: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
178: }
179: }
180: }
181: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
182: }
183: return(0);
184: }
186: PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
187: {
188: MPI_Comm comm;
189: const PetscInt dim = 3;
190: const char *labelName = "marker";
191: ::tetgenio in;
192: ::tetgenio out;
193: DMLabel label;
194: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
195: PetscMPIInt rank;
199: PetscObjectGetComm((PetscObject)dm,&comm);
200: MPI_Comm_rank(comm, &rank);
201: DMPlexGetDepth(dm, &depth);
202: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
203: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
204: DMGetLabel(dm, labelName, &label);
206: in.numberofpoints = vEnd - vStart;
207: if (in.numberofpoints > 0) {
208: PetscSection coordSection;
209: Vec coordinates;
210: PetscScalar *array;
212: in.pointlist = new double[in.numberofpoints*dim];
213: in.pointmarkerlist = new int[in.numberofpoints];
215: DMGetCoordinatesLocal(dm, &coordinates);
216: DMGetCoordinateSection(dm, &coordSection);
217: VecGetArray(coordinates, &array);
218: for (v = vStart; v < vEnd; ++v) {
219: const PetscInt idx = v - vStart;
220: PetscInt off, d;
222: PetscSectionGetOffset(coordSection, v, &off);
223: for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
224: if (label) {
225: PetscInt val;
227: DMLabelGetValue(label, v, &val);
228: in.pointmarkerlist[idx] = (int) val;
229: }
230: }
231: VecRestoreArray(coordinates, &array);
232: }
233: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
235: in.numberofcorners = 4;
236: in.numberoftetrahedra = cEnd - cStart;
237: in.tetrahedronvolumelist = (double*) maxVolumes;
238: if (in.numberoftetrahedra > 0) {
239: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
240: for (c = cStart; c < cEnd; ++c) {
241: const PetscInt idx = c - cStart;
242: PetscInt *closure = NULL;
243: PetscInt closureSize;
245: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
246: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
247: for (v = 0; v < 4; ++v) {
248: in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
249: }
250: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
251: }
252: }
253: /* TODO: Put in boundary faces with markers */
254: if (!rank) {
255: char args[32];
257: #if 1
258: /* Take away 'Q' for verbose output */
259: PetscStrcpy(args, "qezQra");
260: #else
261: PetscStrcpy(args, "qezraVVVV");
262: #endif
263: ::tetrahedralize(args, &in, &out);
264: }
265: in.tetrahedronvolumelist = NULL;
267: {
268: DMLabel rlabel = NULL;
269: const PetscInt numCorners = 4;
270: const PetscInt numCells = out.numberoftetrahedra;
271: const PetscInt numVertices = out.numberofpoints;
272: PetscReal *meshCoords = NULL;
273: PetscInt *cells = NULL;
274: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
276: if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
277: meshCoords = (PetscReal *) out.pointlist;
278: } else {
279: PetscInt i;
281: meshCoords = new PetscReal[dim * numVertices];
282: for (i = 0; i < dim * numVertices; i++) {
283: meshCoords[i] = (PetscReal) out.pointlist[i];
284: }
285: }
286: if (sizeof (PetscInt) == sizeof (out.tetrahedronlist[0])) {
287: cells = (PetscInt *) out.tetrahedronlist;
288: } else {
289: PetscInt i;
291: cells = new PetscInt[numCells * numCorners];
292: for (i = 0; i < numCells * numCorners; i++) {
293: cells[i] = (PetscInt) out.tetrahedronlist[i];
294: }
295: }
297: DMPlexInvertCells_Tetgen(numCells, numCorners, cells);
298: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
299: if (label) {
300: DMCreateLabel(*dmRefined, labelName);
301: DMGetLabel(*dmRefined, labelName, &rlabel);
302: }
303: /* Set labels */
304: for (v = 0; v < numVertices; ++v) {
305: if (out.pointmarkerlist[v]) {
306: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
307: }
308: }
309: if (interpolate) {
310: PetscInt f;
311: #if 0
312: PetscInt e;
314: for (e = 0; e < out.numberofedges; e++) {
315: if (out.edgemarkerlist[e]) {
316: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
317: const PetscInt *edges;
318: PetscInt numEdges;
320: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
321: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
322: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
323: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
324: }
325: }
326: #endif
327: for (f = 0; f < out.numberoftrifaces; f++) {
328: if (out.trifacemarkerlist[f]) {
329: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
330: const PetscInt *faces;
331: PetscInt numFaces;
333: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
334: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
335: if (rlabel) {DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);}
336: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
337: }
338: }
339: }
340: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
341: }
342: return(0);
343: }