Actual source code: ctetgenerate.c
petsc-3.14.6 2021-03-30
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, PetscInt cells[])
7: {
8: PetscInt bound = numCells*numCorners, coff;
11: #define SWAP(a,b) do { PetscInt 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: PetscReal *meshCoords;
118: PetscInt *cells;
120: if (sizeof (PetscReal) == sizeof (out->pointlist[0])) {
121: meshCoords = (PetscReal *) out->pointlist;
122: } else {
123: PetscInt i;
125: PetscMalloc1(dim * numVertices,&meshCoords);
126: for (i = 0; i < dim * numVertices; i++) {
127: meshCoords[i] = (PetscReal) out->pointlist[i];
128: }
129: }
130: if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) {
131: cells = (PetscInt *) out->tetrahedronlist;
132: } else {
133: PetscInt i;
135: PetscMalloc1(numCells * numCorners, &cells);
136: for (i = 0; i < numCells * numCorners; i++) {
137: cells[i] = (PetscInt) out->tetrahedronlist[i];
138: }
139: }
141: DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
142: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
143: if (sizeof (PetscReal) != sizeof (out->pointlist[0])) {
144: PetscFree(meshCoords);
145: }
146: if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) {
147: PetscFree(cells);
148: }
149: if (label) {
150: DMCreateLabel(*dm, labelName);
151: DMGetLabel(*dm, labelName, &glabel);
152: }
153: /* Set labels */
154: for (v = 0; v < numVertices; ++v) {
155: if (out->pointmarkerlist[v]) {
156: if (glabel) {DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);}
157: }
158: }
159: if (interpolate) {
160: PetscInt e;
162: for (e = 0; e < out->numberofedges; e++) {
163: if (out->edgemarkerlist[e]) {
164: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
165: const PetscInt *edges;
166: PetscInt numEdges;
168: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
169: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
170: if (glabel) {DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);}
171: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
172: }
173: }
174: for (f = 0; f < out->numberoftrifaces; f++) {
175: if (out->trifacemarkerlist[f]) {
176: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
177: const PetscInt *faces;
178: PetscInt numFaces;
180: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
181: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
182: if (glabel) {DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);}
183: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
184: }
185: }
186: }
187: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
188: }
190: PLCDestroy(&in);
191: PLCDestroy(&out);
192: return(0);
193: }
195: PETSC_EXTERN PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
196: {
197: MPI_Comm comm;
198: const PetscInt dim = 3;
199: const char *labelName = "marker";
200: PLC *in, *out;
201: DMLabel label;
202: PetscInt verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
203: PetscMPIInt rank;
207: PetscObjectGetComm((PetscObject)dm,&comm);
208: PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
209: MPI_Comm_rank(comm, &rank);
210: DMPlexGetDepth(dm, &depth);
211: MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
212: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
213: DMGetLabel(dm, labelName, &label);
214: PLCCreate(&in);
215: PLCCreate(&out);
217: in->numberofpoints = vEnd - vStart;
218: if (in->numberofpoints > 0) {
219: PetscSection coordSection;
220: Vec coordinates;
221: PetscScalar *array;
223: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
224: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
225: DMGetCoordinatesLocal(dm, &coordinates);
226: DMGetCoordinateSection(dm, &coordSection);
227: VecGetArray(coordinates, &array);
228: for (v = vStart; v < vEnd; ++v) {
229: const PetscInt idx = v - vStart;
230: PetscInt off, d, m = -1;
232: PetscSectionGetOffset(coordSection, v, &off);
233: for (d = 0; d < dim; ++d) {
234: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
235: }
236: if (label) {DMLabelGetValue(label, v, &m);}
238: in->pointmarkerlist[idx] = (int) m;
239: }
240: VecRestoreArray(coordinates, &array);
241: }
242: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
244: in->numberofcorners = 4;
245: in->numberoftetrahedra = cEnd - cStart;
246: in->tetrahedronvolumelist = maxVolumes;
247: if (in->numberoftetrahedra > 0) {
248: PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
249: for (c = cStart; c < cEnd; ++c) {
250: const PetscInt idx = c - cStart;
251: PetscInt *closure = NULL;
252: PetscInt closureSize;
254: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
255: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
256: for (v = 0; v < 4; ++v) {
257: in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
258: }
259: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
260: }
261: }
262: if (!rank) {
263: TetGenOpts t;
265: TetGenOptsInitialize(&t);
267: t.in = dm; /* Should go away */
268: t.refine = 1;
269: t.varvolume = 1;
270: t.quality = 1;
271: t.edgesout = 1;
272: t.zeroindex = 1;
273: t.quiet = 1;
274: t.verbose = verbose; /* Change this */
276: TetGenCheckOpts(&t);
277: TetGenTetrahedralize(&t, in, out);
278: }
279: in->tetrahedronvolumelist = NULL;
280: {
281: DMLabel rlabel = NULL;
282: const PetscInt numCorners = 4;
283: const PetscInt numCells = out->numberoftetrahedra;
284: const PetscInt numVertices = out->numberofpoints;
285: PetscReal *meshCoords;
286: PetscInt *cells;
287: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
289: if (sizeof (PetscReal) == sizeof (out->pointlist[0])) {
290: meshCoords = (PetscReal *) out->pointlist;
291: } else {
292: PetscInt i;
294: PetscMalloc1(dim * numVertices,&meshCoords);
295: for (i = 0; i < dim * numVertices; i++) {
296: meshCoords[i] = (PetscReal) out->pointlist[i];
297: }
298: }
299: if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) {
300: cells = (PetscInt *) out->tetrahedronlist;
301: } else {
302: PetscInt i;
304: PetscMalloc1(numCells * numCorners, &cells);
305: for (i = 0; i < numCells * numCorners; i++) {
306: cells[i] = (PetscInt) out->tetrahedronlist[i];
307: }
308: }
310: DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
311: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
312: if (sizeof (PetscReal) != sizeof (out->pointlist[0])) {
313: PetscFree(meshCoords);
314: }
315: if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) {
316: PetscFree(cells);
317: }
318: if (label) {
319: DMCreateLabel(*dmRefined, labelName);
320: DMGetLabel(*dmRefined, labelName, &rlabel);
321: }
322: /* Set labels */
323: for (v = 0; v < numVertices; ++v) {
324: if (out->pointmarkerlist[v]) {
325: if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);}
326: }
327: }
328: if (interpolate) {
329: PetscInt e, f;
331: for (e = 0; e < out->numberofedges; e++) {
332: if (out->edgemarkerlist[e]) {
333: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
334: const PetscInt *edges;
335: PetscInt numEdges;
337: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
338: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
339: if (rlabel) {DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);}
340: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
341: }
342: }
343: for (f = 0; f < out->numberoftrifaces; f++) {
344: if (out->trifacemarkerlist[f]) {
345: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
346: const PetscInt *faces;
347: PetscInt numFaces;
349: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
350: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
351: if (rlabel) {DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);}
352: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
353: }
354: }
355: }
356: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
357: }
358: PLCDestroy(&in);
359: PLCDestroy(&out);
360: return(0);
361: }