Actual source code: tetgenerate.cxx

petsc-3.11.4 2019-09-28
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>

  3: #include <tetgen.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_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
 19: {
 20:   MPI_Comm       comm;
 21:   DM_Plex       *mesh      = (DM_Plex *) boundary->data;
 22:   const PetscInt dim       = 3;
 23:   const char    *labelName = "marker";
 24:   ::tetgenio     in;
 25:   ::tetgenio     out;
 26:   DMLabel        label;
 27:   PetscInt       vStart, vEnd, v, fStart, fEnd, f;
 28:   PetscMPIInt    rank;

 32:   PetscObjectGetComm((PetscObject)boundary,&comm);
 33:   MPI_Comm_rank(comm, &rank);
 34:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
 35:   DMGetLabel(boundary, labelName, &label);

 37:   in.numberofpoints = vEnd - vStart;
 38:   if (in.numberofpoints > 0) {
 39:     PetscSection coordSection;
 40:     Vec          coordinates;
 41:     PetscScalar *array;

 43:     in.pointlist       = new double[in.numberofpoints*dim];
 44:     in.pointmarkerlist = new int[in.numberofpoints];

 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;

 53:       PetscSectionGetOffset(coordSection, v, &off);
 54:       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
 55:       if (label) {
 56:         PetscInt val;

 58:         DMLabelGetValue(label, v, &val);
 59:         in.pointmarkerlist[idx] = (int) val;
 60:       }
 61:     }
 62:     VecRestoreArray(coordinates, &array);
 63:   }
 64:   DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);

 66:   in.numberoffacets = fEnd - fStart;
 67:   if (in.numberoffacets > 0) {
 68:     in.facetlist       = new tetgenio::facet[in.numberoffacets];
 69:     in.facetmarkerlist = new int[in.numberoffacets];
 70:     for (f = fStart; f < fEnd; ++f) {
 71:       const PetscInt idx     = f - fStart;
 72:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v;

 74:       in.facetlist[idx].numberofpolygons = 1;
 75:       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
 76:       in.facetlist[idx].numberofholes    = 0;
 77:       in.facetlist[idx].holelist         = NULL;

 79:       DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
 80:       for (p = 0; p < numPoints*2; p += 2) {
 81:         const PetscInt point = points[p];
 82:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
 83:       }

 85:       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
 86:       poly->numberofvertices = numVertices;
 87:       poly->vertexlist       = new int[poly->numberofvertices];
 88:       for (v = 0; v < numVertices; ++v) {
 89:         const PetscInt vIdx = points[v] - vStart;
 90:         poly->vertexlist[v] = vIdx;
 91:       }
 92:       if (label) {
 93:         PetscInt val;

 95:         DMLabelGetValue(label, f, &val);
 96:         in.facetmarkerlist[idx] = (int) val;
 97:       }
 98:       DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
 99:     }
100:   }
101:   if (!rank) {
102:     char args[32];

104:     /* Take away 'Q' for verbose output */
105:     PetscStrcpy(args, "pqezQ");
106:     if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
107:     else                  {::tetrahedralize(args, &in, &out);}
108:   }
109:   {
110:     DMLabel        glabel      = NULL;
111:     const PetscInt numCorners  = 4;
112:     const PetscInt numCells    = out.numberoftetrahedra;
113:     const PetscInt numVertices = out.numberofpoints;
114:     const double   *meshCoords = out.pointlist;
115:     int            *cells      = out.tetrahedronlist;

117:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
118:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
119:     if (label) {DMCreateLabel(*dm, labelName); DMGetLabel(*dm, labelName, &glabel);}
120:     /* Set labels */
121:     for (v = 0; v < numVertices; ++v) {
122:       if (out.pointmarkerlist[v]) {
123:         if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
124:       }
125:     }
126:     if (interpolate) {
127: #if 0
128:       PetscInt e;

130:       /* This check is never actually executed for ctetgen (which never returns edgemarkers) and seems to be broken for
131:        * tetgen */
132:       for (e = 0; e < out.numberofedges; e++) {
133:         if (out.edgemarkerlist[e]) {
134:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
135:           const PetscInt *edges;
136:           PetscInt        numEdges;

138:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
139:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
140:           if (glabel) {DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);}
141:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
142:         }
143:       }
144: #endif
145:       for (f = 0; f < out.numberoftrifaces; f++) {
146:         if (out.trifacemarkerlist[f]) {
147:           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
148:           const PetscInt *faces;
149:           PetscInt        numFaces;

151:           DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
152:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
153:           if (glabel) {DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);}
154:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
155:         }
156:       }
157:     }
158:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
159:   }
160:   return(0);
161: }

163: PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
164: {
165:   MPI_Comm       comm;
166:   const PetscInt dim       = 3;
167:   const char    *labelName = "marker";
168:   ::tetgenio     in;
169:   ::tetgenio     out;
170:   DMLabel        label;
171:   PetscInt       vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
172:   PetscMPIInt    rank;

176:   PetscObjectGetComm((PetscObject)dm,&comm);
177:   MPI_Comm_rank(comm, &rank);
178:   DMPlexGetDepth(dm, &depth);
179:   MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
180:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
181:   DMGetLabel(dm, labelName, &label);

183:   in.numberofpoints = vEnd - vStart;
184:   if (in.numberofpoints > 0) {
185:     PetscSection coordSection;
186:     Vec          coordinates;
187:     PetscScalar *array;

189:     in.pointlist       = new double[in.numberofpoints*dim];
190:     in.pointmarkerlist = new int[in.numberofpoints];

192:     DMGetCoordinatesLocal(dm, &coordinates);
193:     DMGetCoordinateSection(dm, &coordSection);
194:     VecGetArray(coordinates, &array);
195:     for (v = vStart; v < vEnd; ++v) {
196:       const PetscInt idx = v - vStart;
197:       PetscInt       off, d;

199:       PetscSectionGetOffset(coordSection, v, &off);
200:       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
201:       if (label) {
202:         PetscInt val;
203: 
204:         DMLabelGetValue(label, v, &val);
205:         in.pointmarkerlist[idx] = (int) val;
206:       }
207:     }
208:     VecRestoreArray(coordinates, &array);
209:   }
210:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

212:   in.numberofcorners       = 4;
213:   in.numberoftetrahedra    = cEnd - cStart;
214:   in.tetrahedronvolumelist = (double*) maxVolumes;
215:   if (in.numberoftetrahedra > 0) {
216:     in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
217:     for (c = cStart; c < cEnd; ++c) {
218:       const PetscInt idx      = c - cStart;
219:       PetscInt      *closure = NULL;
220:       PetscInt       closureSize;

222:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
223:       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
224:       for (v = 0; v < 4; ++v) {
225:         in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
226:       }
227:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
228:     }
229:   }
230:   /* TODO: Put in boundary faces with markers */
231:   if (!rank) {
232:     char args[32];

234: #if 1
235:     /* Take away 'Q' for verbose output */
236:     PetscStrcpy(args, "qezQra");
237: #else
238:     PetscStrcpy(args, "qezraVVVV");
239: #endif
240:     ::tetrahedralize(args, &in, &out);
241:   }
242:   in.tetrahedronvolumelist = NULL;

244:   {
245:     DMLabel        rlabel      = NULL;
246:     const PetscInt numCorners  = 4;
247:     const PetscInt numCells    = out.numberoftetrahedra;
248:     const PetscInt numVertices = out.numberofpoints;
249:     const double   *meshCoords = out.pointlist;
250:     int            *cells      = out.tetrahedronlist;

252:     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;

254:     DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);
255:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
256:     if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
257:     /* Set labels */
258:     for (v = 0; v < numVertices; ++v) {
259:       if (out.pointmarkerlist[v]) {
260:         if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
261:       }
262:     }
263:     if (interpolate) {
264:       PetscInt f;
265: #if 0
266:       PetscInt e;

268:       for (e = 0; e < out.numberofedges; e++) {
269:         if (out.edgemarkerlist[e]) {
270:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
271:           const PetscInt *edges;
272:           PetscInt        numEdges;

274:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
275:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
276:           if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
277:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
278:         }
279:       }
280: #endif
281:       for (f = 0; f < out.numberoftrifaces; f++) {
282:         if (out.trifacemarkerlist[f]) {
283:           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
284:           const PetscInt *faces;
285:           PetscInt        numFaces;

287:           DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
288:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
289:           if (rlabel) {DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);}
290:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
291:         }
292:       }
293:     }
294:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
295:   }
296:   return(0);
297: }