Actual source code: ctetgenerate.c

petsc-3.10.5 2019-03-28
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_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: }