Actual source code: ctetgenerate.c

petsc-master 2020-07-09
Report Typos and Errors
  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: }