Actual source code: trigenerate.c

  1: #include <petsc/private/dmpleximpl.h>

  3: #if !defined(ANSI_DECLARATORS)
  4:   #define ANSI_DECLARATORS
  5: #endif
  6: #include <triangle.h>

  8: static PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx)
  9: {
 10:   inputCtx->numberofpoints             = 0;
 11:   inputCtx->numberofpointattributes    = 0;
 12:   inputCtx->pointlist                  = NULL;
 13:   inputCtx->pointattributelist         = NULL;
 14:   inputCtx->pointmarkerlist            = NULL;
 15:   inputCtx->numberofsegments           = 0;
 16:   inputCtx->segmentlist                = NULL;
 17:   inputCtx->segmentmarkerlist          = NULL;
 18:   inputCtx->numberoftriangleattributes = 0;
 19:   inputCtx->trianglelist               = NULL;
 20:   inputCtx->numberofholes              = 0;
 21:   inputCtx->holelist                   = NULL;
 22:   inputCtx->numberofregions            = 0;
 23:   inputCtx->regionlist                 = NULL;
 24:   return 0;
 25: }

 27: static PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx)
 28: {
 29:   outputCtx->numberofpoints        = 0;
 30:   outputCtx->pointlist             = NULL;
 31:   outputCtx->pointattributelist    = NULL;
 32:   outputCtx->pointmarkerlist       = NULL;
 33:   outputCtx->numberoftriangles     = 0;
 34:   outputCtx->trianglelist          = NULL;
 35:   outputCtx->triangleattributelist = NULL;
 36:   outputCtx->neighborlist          = NULL;
 37:   outputCtx->segmentlist           = NULL;
 38:   outputCtx->segmentmarkerlist     = NULL;
 39:   outputCtx->numberofedges         = 0;
 40:   outputCtx->edgelist              = NULL;
 41:   outputCtx->edgemarkerlist        = NULL;
 42:   return 0;
 43: }

 45: static PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
 46: {
 47:   free(outputCtx->pointlist);
 48:   free(outputCtx->pointmarkerlist);
 49:   free(outputCtx->segmentlist);
 50:   free(outputCtx->segmentmarkerlist);
 51:   free(outputCtx->edgelist);
 52:   free(outputCtx->edgemarkerlist);
 53:   free(outputCtx->trianglelist);
 54:   free(outputCtx->neighborlist);
 55:   return 0;
 56: }

 58: PETSC_EXTERN PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
 59: {
 60:   MPI_Comm             comm;
 61:   DM_Plex             *mesh             = (DM_Plex *)boundary->data;
 62:   PetscInt             dim              = 2;
 63:   const PetscBool      createConvexHull = PETSC_FALSE;
 64:   const PetscBool      constrained      = PETSC_FALSE;
 65:   const char          *labelName        = "marker";
 66:   const char          *labelName2       = "Face Sets";
 67:   struct triangulateio in;
 68:   struct triangulateio out;
 69:   DMLabel              label, label2;
 70:   PetscInt             vStart, vEnd, v, eStart, eEnd, e;
 71:   PetscMPIInt          rank;

 73:   PetscObjectGetComm((PetscObject)boundary, &comm);
 74:   MPI_Comm_rank(comm, &rank);
 75:   InitInput_Triangle(&in);
 76:   InitOutput_Triangle(&out);
 77:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
 78:   DMGetLabel(boundary, labelName, &label);
 79:   DMGetLabel(boundary, labelName2, &label2);

 81:   in.numberofpoints = vEnd - vStart;
 82:   if (in.numberofpoints > 0) {
 83:     PetscSection coordSection;
 84:     Vec          coordinates;
 85:     PetscScalar *array;

 87:     PetscMalloc1(in.numberofpoints * dim, &in.pointlist);
 88:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
 89:     DMGetCoordinatesLocal(boundary, &coordinates);
 90:     DMGetCoordinateSection(boundary, &coordSection);
 91:     VecGetArray(coordinates, &array);
 92:     for (v = vStart; v < vEnd; ++v) {
 93:       const PetscInt idx = v - vStart;
 94:       PetscInt       val, off, d;

 96:       PetscSectionGetOffset(coordSection, v, &off);
 97:       for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
 98:       if (label) {
 99:         DMLabelGetValue(label, v, &val);
100:         in.pointmarkerlist[idx] = val;
101:       }
102:     }
103:     VecRestoreArray(coordinates, &array);
104:   }
105:   DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);
106:   in.numberofsegments = eEnd - eStart;
107:   if (in.numberofsegments > 0) {
108:     PetscMalloc1(in.numberofsegments * 2, &in.segmentlist);
109:     PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);
110:     for (e = eStart; e < eEnd; ++e) {
111:       const PetscInt  idx = e - eStart;
112:       const PetscInt *cone;
113:       PetscInt        val;

115:       DMPlexGetCone(boundary, e, &cone);

117:       in.segmentlist[idx * 2 + 0] = cone[0] - vStart;
118:       in.segmentlist[idx * 2 + 1] = cone[1] - vStart;

120:       if (label) {
121:         DMLabelGetValue(label, e, &val);
122:         in.segmentmarkerlist[idx] = val;
123:       }
124:     }
125:   }
126: #if 0 /* Do not currently support holes */
127:   PetscReal *holeCoords;
128:   PetscInt   h, d;

130:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
131:   if (in.numberofholes > 0) {
132:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
133:     for (h = 0; h < in.numberofholes; ++h) {
134:       for (d = 0; d < dim; ++d) {
135:         in.holelist[h*dim+d] = holeCoords[h*dim+d];
136:       }
137:     }
138:   }
139: #endif
140:   if (rank == 0) {
141:     char args[32];

143:     /* Take away 'Q' for verbose output */
144:     PetscStrcpy(args, "pqezQ");
145:     if (createConvexHull) PetscStrcat(args, "c");
146:     if (constrained) PetscStrcpy(args, "zepDQ");
147:     if (mesh->triangleOpts) {
148:       triangulate(mesh->triangleOpts, &in, &out, NULL);
149:     } else {
150:       triangulate(args, &in, &out, NULL);
151:     }
152:   }
153:   PetscFree(in.pointlist);
154:   PetscFree(in.pointmarkerlist);
155:   PetscFree(in.segmentlist);
156:   PetscFree(in.segmentmarkerlist);
157:   PetscFree(in.holelist);

159:   {
160:     DMLabel        glabel      = NULL;
161:     DMLabel        glabel2     = NULL;
162:     const PetscInt numCorners  = 3;
163:     const PetscInt numCells    = out.numberoftriangles;
164:     const PetscInt numVertices = out.numberofpoints;
165:     PetscInt      *cells;
166:     PetscReal     *meshCoords;

168:     if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
169:       meshCoords = (PetscReal *)out.pointlist;
170:     } else {
171:       PetscInt i;

173:       PetscMalloc1(dim * numVertices, &meshCoords);
174:       for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
175:     }
176:     if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
177:       cells = (PetscInt *)out.trianglelist;
178:     } else {
179:       PetscInt i;

181:       PetscMalloc1(numCells * numCorners, &cells);
182:       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
183:     }
184:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
185:     if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscFree(meshCoords);
186:     if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscFree(cells);
187:     if (label) {
188:       DMCreateLabel(*dm, labelName);
189:       DMGetLabel(*dm, labelName, &glabel);
190:     }
191:     if (label2) {
192:       DMCreateLabel(*dm, labelName2);
193:       DMGetLabel(*dm, labelName2, &glabel2);
194:     }
195:     /* Set labels */
196:     for (v = 0; v < numVertices; ++v) {
197:       if (out.pointmarkerlist[v]) {
198:         if (glabel) DMLabelSetValue(glabel, v + numCells, out.pointmarkerlist[v]);
199:       }
200:     }
201:     if (interpolate) {
202:       for (e = 0; e < out.numberofedges; e++) {
203:         if (out.edgemarkerlist[e]) {
204:           const PetscInt  vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
205:           const PetscInt *edges;
206:           PetscInt        numEdges;

208:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
210:           if (glabel) DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);
211:           if (glabel2) DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);
212:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
213:         }
214:       }
215:     }
216:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
217:   }
218: #if 0 /* Do not currently support holes */
219:   DMPlexCopyHoles(*dm, boundary);
220: #endif
221:   FiniOutput_Triangle(&out);
222:   return 0;
223: }

225: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
226: {
227:   MPI_Comm             comm;
228:   PetscInt             dim       = 2;
229:   const char          *labelName = "marker";
230:   struct triangulateio in;
231:   struct triangulateio out;
232:   DMLabel              label;
233:   PetscInt             vStart, vEnd, v, gcStart, cStart, cEnd, c, depth, depthGlobal;
234:   PetscMPIInt          rank;
235:   double              *maxVolumes;

237:   PetscObjectGetComm((PetscObject)dm, &comm);
238:   MPI_Comm_rank(comm, &rank);
239:   InitInput_Triangle(&in);
240:   InitOutput_Triangle(&out);
241:   DMPlexGetDepth(dm, &depth);
242:   MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
243:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
244:   DMGetLabel(dm, labelName, &label);

246:   in.numberofpoints = vEnd - vStart;
247:   if (in.numberofpoints > 0) {
248:     PetscSection coordSection;
249:     Vec          coordinates;
250:     PetscScalar *array;

252:     PetscMalloc1(in.numberofpoints * dim, &in.pointlist);
253:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
254:     DMGetCoordinatesLocal(dm, &coordinates);
255:     DMGetCoordinateSection(dm, &coordSection);
256:     VecGetArray(coordinates, &array);
257:     for (v = vStart; v < vEnd; ++v) {
258:       const PetscInt idx = v - vStart;
259:       PetscInt       off, d, val;

261:       PetscSectionGetOffset(coordSection, v, &off);
262:       for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
263:       if (label) {
264:         DMLabelGetValue(label, v, &val);
265:         in.pointmarkerlist[idx] = val;
266:       }
267:     }
268:     VecRestoreArray(coordinates, &array);
269:   }
270:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
271:   DMPlexGetGhostCellStratum(dm, &gcStart, NULL);
272:   if (gcStart >= 0) cEnd = gcStart;

274:   in.numberofcorners   = 3;
275:   in.numberoftriangles = cEnd - cStart;

277: #if !defined(PETSC_USE_REAL_DOUBLE)
278:   PetscMalloc1(cEnd - cStart, &maxVolumes);
279:   for (c = 0; c < cEnd - cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
280: #else
281:   maxVolumes = inmaxVolumes;
282: #endif

284:   in.trianglearealist = (double *)maxVolumes;
285:   if (in.numberoftriangles > 0) {
286:     PetscMalloc1(in.numberoftriangles * in.numberofcorners, &in.trianglelist);
287:     for (c = cStart; c < cEnd; ++c) {
288:       const PetscInt idx     = c - cStart;
289:       PetscInt      *closure = NULL;
290:       PetscInt       closureSize;

292:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
294:       for (v = 0; v < 3; ++v) in.trianglelist[idx * in.numberofcorners + v] = closure[(v + closureSize - 3) * 2] - vStart;
295:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
296:     }
297:   }
298:   /* TODO: Segment markers are missing on input */
299: #if 0 /* Do not currently support holes */
300:   PetscReal *holeCoords;
301:   PetscInt   h, d;

303:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
304:   if (in.numberofholes > 0) {
305:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
306:     for (h = 0; h < in.numberofholes; ++h) {
307:       for (d = 0; d < dim; ++d) {
308:         in.holelist[h*dim+d] = holeCoords[h*dim+d];
309:       }
310:     }
311:   }
312: #endif
313:   if (rank == 0) {
314:     char args[32];

316:     /* Take away 'Q' for verbose output */
317:     PetscStrcpy(args, "pqezQra");
318:     triangulate(args, &in, &out, NULL);
319:   }
320:   PetscFree(in.pointlist);
321:   PetscFree(in.pointmarkerlist);
322:   PetscFree(in.segmentlist);
323:   PetscFree(in.segmentmarkerlist);
324:   PetscFree(in.trianglelist);

326:   {
327:     DMLabel        rlabel      = NULL;
328:     const PetscInt numCorners  = 3;
329:     const PetscInt numCells    = out.numberoftriangles;
330:     const PetscInt numVertices = out.numberofpoints;
331:     PetscInt      *cells;
332:     PetscReal     *meshCoords;
333:     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;

335:     if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
336:       meshCoords = (PetscReal *)out.pointlist;
337:     } else {
338:       PetscInt i;

340:       PetscMalloc1(dim * numVertices, &meshCoords);
341:       for (i = 0; i < dim * numVertices; i++) meshCoords[i] = (PetscReal)out.pointlist[i];
342:     }
343:     if (sizeof(PetscInt) == sizeof(out.trianglelist[0])) {
344:       cells = (PetscInt *)out.trianglelist;
345:     } else {
346:       PetscInt i;

348:       PetscMalloc1(numCells * numCorners, &cells);
349:       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.trianglelist[i];
350:     }

352:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
353:     if (label) {
354:       DMCreateLabel(*dmRefined, labelName);
355:       DMGetLabel(*dmRefined, labelName, &rlabel);
356:     }
357:     if (sizeof(PetscReal) != sizeof(out.pointlist[0])) PetscFree(meshCoords);
358:     if (sizeof(PetscInt) != sizeof(out.trianglelist[0])) PetscFree(cells);
359:     /* Set labels */
360:     for (v = 0; v < numVertices; ++v) {
361:       if (out.pointmarkerlist[v]) {
362:         if (rlabel) DMLabelSetValue(rlabel, v + numCells, out.pointmarkerlist[v]);
363:       }
364:     }
365:     if (interpolate) {
366:       PetscInt e;

368:       for (e = 0; e < out.numberofedges; e++) {
369:         if (out.edgemarkerlist[e]) {
370:           const PetscInt  vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
371:           const PetscInt *edges;
372:           PetscInt        numEdges;

374:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
376:           if (rlabel) DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);
377:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
378:         }
379:       }
380:     }
381:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
382:   }
383: #if 0 /* Do not currently support holes */
384:   DMPlexCopyHoles(*dm, boundary);
385: #endif
386:   FiniOutput_Triangle(&out);
387: #if !defined(PETSC_USE_REAL_DOUBLE)
388:   PetscFree(maxVolumes);
389: #endif
390:   return 0;
391: }