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