Actual source code: tetgenerate.cxx

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

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

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

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

 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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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