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