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

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

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

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

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

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

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

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

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

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

215:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
217:           if (glabel)  DMLabelSetValue(glabel,  edges[0], out.edgemarkerlist[e]);
218:           if (glabel2) DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);
219:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
220:         }
221:       }
222:     }
223:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
224:   }
225: #if 0 /* Do not currently support holes */
226:   DMPlexCopyHoles(*dm, boundary);
227: #endif
228:   FiniOutput_Triangle(&out);
229:   return 0;
230: }

232: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
233: {
234:   MPI_Comm             comm;
235:   PetscInt             dim       = 2;
236:   const char          *labelName = "marker";
237:   struct triangulateio in;
238:   struct triangulateio out;
239:   DMLabel              label;
240:   PetscInt             vStart, vEnd, v, gcStart, cStart, cEnd, c, depth, depthGlobal;
241:   PetscMPIInt          rank;
242:   double               *maxVolumes;

244:   PetscObjectGetComm((PetscObject)dm,&comm);
245:   MPI_Comm_rank(comm, &rank);
246:   InitInput_Triangle(&in);
247:   InitOutput_Triangle(&out);
248:   DMPlexGetDepth(dm, &depth);
249:   MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
250:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
251:   DMGetLabel(dm, labelName, &label);

253:   in.numberofpoints = vEnd - vStart;
254:   if (in.numberofpoints > 0) {
255:     PetscSection coordSection;
256:     Vec          coordinates;
257:     PetscScalar *array;

259:     PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
260:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
261:     DMGetCoordinatesLocal(dm, &coordinates);
262:     DMGetCoordinateSection(dm, &coordSection);
263:     VecGetArray(coordinates, &array);
264:     for (v = vStart; v < vEnd; ++v) {
265:       const PetscInt idx = v - vStart;
266:       PetscInt       off, d, val;

268:       PetscSectionGetOffset(coordSection, v, &off);
269:       for (d = 0; d < dim; ++d) {
270:         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
271:       }
272:       if (label) {
273:         DMLabelGetValue(label, v, &val);
274:         in.pointmarkerlist[idx] = val;
275:       }
276:     }
277:     VecRestoreArray(coordinates, &array);
278:   }
279:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
280:   DMPlexGetGhostCellStratum(dm, &gcStart, NULL);
281:   if (gcStart >= 0) cEnd = gcStart;

283:   in.numberofcorners   = 3;
284:   in.numberoftriangles = cEnd - cStart;

286: #if !defined(PETSC_USE_REAL_DOUBLE)
287:   PetscMalloc1(cEnd - cStart,&maxVolumes);
288:   for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
289: #else
290:   maxVolumes = inmaxVolumes;
291: #endif

293:   in.trianglearealist  = (double*) maxVolumes;
294:   if (in.numberoftriangles > 0) {
295:     PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
296:     for (c = cStart; c < cEnd; ++c) {
297:       const PetscInt idx      = c - cStart;
298:       PetscInt      *closure = NULL;
299:       PetscInt       closureSize;

301:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
303:       for (v = 0; v < 3; ++v) {
304:         in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
305:       }
306:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
307:     }
308:   }
309:   /* TODO: Segment markers are missing on input */
310: #if 0 /* Do not currently support holes */
311:   PetscReal *holeCoords;
312:   PetscInt   h, d;

314:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
315:   if (in.numberofholes > 0) {
316:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
317:     for (h = 0; h < in.numberofholes; ++h) {
318:       for (d = 0; d < dim; ++d) {
319:         in.holelist[h*dim+d] = holeCoords[h*dim+d];
320:       }
321:     }
322:   }
323: #endif
324:   if (rank == 0) {
325:     char args[32];

327:     /* Take away 'Q' for verbose output */
328:     PetscStrcpy(args, "pqezQra");
329:     triangulate(args, &in, &out, NULL);
330:   }
331:   PetscFree(in.pointlist);
332:   PetscFree(in.pointmarkerlist);
333:   PetscFree(in.segmentlist);
334:   PetscFree(in.segmentmarkerlist);
335:   PetscFree(in.trianglelist);

337:   {
338:     DMLabel          rlabel      = NULL;
339:     const PetscInt   numCorners  = 3;
340:     const PetscInt   numCells    = out.numberoftriangles;
341:     const PetscInt   numVertices = out.numberofpoints;
342:     PetscInt         *cells;
343:     PetscReal        *meshCoords;
344:     PetscBool        interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;

346:     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
347:       meshCoords = (PetscReal *) out.pointlist;
348:     } else {
349:       PetscInt i;

351:       PetscMalloc1(dim * numVertices,&meshCoords);
352:       for (i = 0; i < dim * numVertices; i++) {
353:         meshCoords[i] = (PetscReal) out.pointlist[i];
354:       }
355:     }
356:     if (sizeof (PetscInt) == sizeof (out.trianglelist[0])) {
357:       cells = (PetscInt *) out.trianglelist;
358:     } else {
359:       PetscInt i;

361:       PetscMalloc1(numCells * numCorners, &cells);
362:       for (i = 0; i < numCells * numCorners; i++) {
363:         cells[i] = (PetscInt) out.trianglelist[i];
364:       }
365:     }

367:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
368:     if (label) {
369:       DMCreateLabel(*dmRefined, labelName);
370:       DMGetLabel(*dmRefined, labelName, &rlabel);
371:     }
372:     if (sizeof (PetscReal) != sizeof (out.pointlist[0])) {
373:       PetscFree(meshCoords);
374:     }
375:     if (sizeof (PetscInt) != sizeof (out.trianglelist[0])) {
376:       PetscFree(cells);
377:     }
378:     /* Set labels */
379:     for (v = 0; v < numVertices; ++v) {
380:       if (out.pointmarkerlist[v]) {
381:         if (rlabel) DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);
382:       }
383:     }
384:     if (interpolate) {
385:       PetscInt e;

387:       for (e = 0; e < out.numberofedges; e++) {
388:         if (out.edgemarkerlist[e]) {
389:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
390:           const PetscInt *edges;
391:           PetscInt        numEdges;

393:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
395:           if (rlabel) DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);
396:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
397:         }
398:       }
399:     }
400:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
401:   }
402: #if 0 /* Do not currently support holes */
403:   DMPlexCopyHoles(*dm, boundary);
404: #endif
405:   FiniOutput_Triangle(&out);
406: #if !defined(PETSC_USE_REAL_DOUBLE)
407:   PetscFree(maxVolumes);
408: #endif
409:   return 0;
410: }