Actual source code: ctetgenerate.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/dmpleximpl.h>
3: #include <ctetgen.h>
5: /* This is to fix the tetrahedron orientation from TetGen */
6: static PetscErrorCode DMPlexInvertCells_CTetgen(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_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
18: {
19: MPI_Comm comm;
20: const PetscInt dim = 3;
21: const char *labelName = "marker";
22: PLC *in, *out;
23: DMLabel label;
24: PetscInt verbose = 0, vStart, vEnd, v, fStart, fEnd, f;
25: PetscMPIInt rank;
29: PetscObjectGetComm((PetscObject)boundary,&comm);
30: PetscOptionsGetInt(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
31: MPI_Comm_rank(comm, &rank);
32: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
33: DMGetLabel(boundary, labelName, &label);
34: PLCCreate(&in);
35: PLCCreate(&out);
37: in->numberofpoints = vEnd - vStart;
38: if (in->numberofpoints > 0) {
39: PetscSection coordSection;
40: Vec coordinates;
41: PetscScalar *array;
43: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
44: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
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, m = -1;
52: PetscSectionGetOffset(coordSection, v, &off);
53: for (d = 0; d < dim; ++d) {
54: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
55: }
56: if (label) {DMLabelGetValue(label, v, &m);}
58: in->pointmarkerlist[idx] = (int) m;
59: }
60: VecRestoreArray(coordinates, &array);
61: }
62: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
64: in->numberoffacets = fEnd - fStart;
65: if (in->numberoffacets > 0) {
66: PetscMalloc1(in->numberoffacets, &in->facetlist);
67: PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);
68: for (f = fStart; f < fEnd; ++f) {
69: const PetscInt idx = f - fStart;
70: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
71: polygon *poly;
73: in->facetlist[idx].numberofpolygons = 1;
75: PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);
77: in->facetlist[idx].numberofholes = 0;
78: in->facetlist[idx].holelist = NULL;
80: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
81: for (p = 0; p < numPoints*2; p += 2) {
82: const PetscInt point = points[p];
83: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
84: }
86: poly = in->facetlist[idx].polygonlist;
87: poly->numberofvertices = numVertices;
88: PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
89: for (v = 0; v < numVertices; ++v) {
90: const PetscInt vIdx = points[v] - vStart;
91: poly->vertexlist[v] = vIdx;
92: }
93: if (label) {DMLabelGetValue(label, f, &m);}
94: in->facetmarkerlist[idx] = (int) m;
95: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
96: }
97: }
98: if (!rank) {
99: TetGenOpts t;
101: TetGenOptsInitialize(&t);
102: t.in = boundary; /* Should go away */
103: t.plc = 1;
104: t.quality = 1;
105: t.edgesout = 1;
106: t.zeroindex = 1;
107: t.quiet = 1;
108: t.verbose = verbose;
109: TetGenCheckOpts(&t);
110: TetGenTetrahedralize(&t, in, out);
111: }
112: {
113: DMLabel glabel = NULL;
114: const PetscInt numCorners = 4;
115: const PetscInt numCells = out->numberoftetrahedra;
116: const PetscInt numVertices = out->numberofpoints;
117: double *meshCoords;
118: int *cells = out->tetrahedronlist;
120: if (sizeof (PetscReal) == sizeof (double)) {
121: meshCoords = (double *) out->pointlist;
122: }
123: else {
124: PetscInt i;
126: PetscMalloc1(3 * numVertices,&meshCoords);
127: for (i = 0; i < 3 * numVertices; i++) {
128: meshCoords[i] = (PetscReal) out->pointlist[i];
129: }
130: }
132: DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
133: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
134: if (sizeof (PetscReal) != sizeof (double)) {
135: PetscFree(meshCoords);
136: }
137: if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
138: /* Set labels */
139: for (v = 0; v < numVertices; ++v) {
140: if (out->pointmarkerlist[v]) {
141: if (glabel) {DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);}
142: }
143: }
144: if (interpolate) {
145: PetscInt e;
147: for (e = 0; e < out->numberofedges; e++) {
148: if (out->edgemarkerlist[e]) {
149: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
150: const PetscInt *edges;
151: PetscInt numEdges;
153: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
154: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
155: if (glabel) {DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);}
156: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
157: }
158: }
159: for (f = 0; f < out->numberoftrifaces; f++) {
160: if (out->trifacemarkerlist[f]) {
161: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
162: const PetscInt *faces;
163: PetscInt numFaces;
165: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
166: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
167: if (glabel) {DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);}
168: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
169: }
170: }
171: }
172: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
173: }
175: PLCDestroy(&in);
176: PLCDestroy(&out);
177: return(0);
178: }
180: PETSC_EXTERN PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
181: {
182: MPI_Comm comm;
183: const PetscInt dim = 3;
184: const char *labelName = "marker";
185: PLC *in, *out;
186: DMLabel label;
187: PetscInt verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
188: PetscMPIInt rank;
192: PetscObjectGetComm((PetscObject)dm,&comm);
193: PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
194: MPI_Comm_rank(comm, &rank);
195: DMPlexGetDepth(dm, &depth);
196: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
197: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
198: DMGetLabel(dm, labelName, &label);
199: PLCCreate(&in);
200: PLCCreate(&out);
202: in->numberofpoints = vEnd - vStart;
203: if (in->numberofpoints > 0) {
204: PetscSection coordSection;
205: Vec coordinates;
206: PetscScalar *array;
208: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
209: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
210: DMGetCoordinatesLocal(dm, &coordinates);
211: DMGetCoordinateSection(dm, &coordSection);
212: VecGetArray(coordinates, &array);
213: for (v = vStart; v < vEnd; ++v) {
214: const PetscInt idx = v - vStart;
215: PetscInt off, d, m = -1;
217: PetscSectionGetOffset(coordSection, v, &off);
218: for (d = 0; d < dim; ++d) {
219: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
220: }
221: if (label) {DMLabelGetValue(label, v, &m);}
223: in->pointmarkerlist[idx] = (int) m;
224: }
225: VecRestoreArray(coordinates, &array);
226: }
227: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
229: in->numberofcorners = 4;
230: in->numberoftetrahedra = cEnd - cStart;
231: in->tetrahedronvolumelist = maxVolumes;
232: if (in->numberoftetrahedra > 0) {
233: PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
234: for (c = cStart; c < cEnd; ++c) {
235: const PetscInt idx = c - cStart;
236: PetscInt *closure = NULL;
237: PetscInt closureSize;
239: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
240: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
241: for (v = 0; v < 4; ++v) {
242: in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
243: }
244: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
245: }
246: }
247: if (!rank) {
248: TetGenOpts t;
250: TetGenOptsInitialize(&t);
252: t.in = dm; /* Should go away */
253: t.refine = 1;
254: t.varvolume = 1;
255: t.quality = 1;
256: t.edgesout = 1;
257: t.zeroindex = 1;
258: t.quiet = 1;
259: t.verbose = verbose; /* Change this */
261: TetGenCheckOpts(&t);
262: TetGenTetrahedralize(&t, in, out);
263: }
264: in->tetrahedronvolumelist = NULL;
265: {
266: DMLabel rlabel = NULL;
267: const PetscInt numCorners = 4;
268: const PetscInt numCells = out->numberoftetrahedra;
269: const PetscInt numVertices = out->numberofpoints;
270: double *meshCoords;
271: int *cells = out->tetrahedronlist;
272: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
274: if (sizeof (PetscReal) == sizeof (double)) {
275: meshCoords = (double *) out->pointlist;
276: }
277: else {
278: PetscInt i;
280: PetscMalloc1(3 * numVertices,&meshCoords);
281: for (i = 0; i < 3 * numVertices; i++) {
282: meshCoords[i] = (PetscReal) out->pointlist[i];
283: }
284: }
286: DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
287: DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
288: if (sizeof (PetscReal) != sizeof (double)) {
289: PetscFree(meshCoords);
290: }
291: if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
292: /* Set labels */
293: for (v = 0; v < numVertices; ++v) {
294: if (out->pointmarkerlist[v]) {
295: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);}
296: }
297: }
298: if (interpolate) {
299: PetscInt e, f;
301: for (e = 0; e < out->numberofedges; e++) {
302: if (out->edgemarkerlist[e]) {
303: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
304: const PetscInt *edges;
305: PetscInt numEdges;
307: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
308: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
309: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);}
310: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
311: }
312: }
313: for (f = 0; f < out->numberoftrifaces; f++) {
314: if (out->trifacemarkerlist[f]) {
315: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
316: const PetscInt *faces;
317: PetscInt numFaces;
319: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
320: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
321: if (rlabel) {DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);}
322: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
323: }
324: }
325: }
326: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
327: }
328: PLCDestroy(&in);
329: PLCDestroy(&out);
330: return(0);
331: }