Actual source code: trigenerate.c

petsc-master 2021-02-24
Report Typos and Errors
  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: {
 11:   inputCtx->numberofpoints             = 0;
 12:   inputCtx->numberofpointattributes    = 0;
 13:   inputCtx->pointlist                  = NULL;
 14:   inputCtx->pointattributelist         = NULL;
 15:   inputCtx->pointmarkerlist            = NULL;
 16:   inputCtx->numberofsegments           = 0;
 17:   inputCtx->segmentlist                = NULL;
 18:   inputCtx->segmentmarkerlist          = NULL;
 19:   inputCtx->numberoftriangleattributes = 0;
 20:   inputCtx->trianglelist               = NULL;
 21:   inputCtx->numberofholes              = 0;
 22:   inputCtx->holelist                   = NULL;
 23:   inputCtx->numberofregions            = 0;
 24:   inputCtx->regionlist                 = NULL;
 25:   return(0);
 26: }

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

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

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

 78:   PetscObjectGetComm((PetscObject)boundary,&comm);
 79:   MPI_Comm_rank(comm, &rank);
 80:   InitInput_Triangle(&in);
 81:   InitOutput_Triangle(&out);
 82:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
 83:   DMGetLabel(boundary, labelName,  &label);
 84:   DMGetLabel(boundary, labelName2, &label2);

 86:   in.numberofpoints = vEnd - vStart;
 87:   if (in.numberofpoints > 0) {
 88:     PetscSection coordSection;
 89:     Vec          coordinates;
 90:     PetscScalar *array;

 92:     PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
 93:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
 94:     DMGetCoordinatesLocal(boundary, &coordinates);
 95:     DMGetCoordinateSection(boundary, &coordSection);
 96:     VecGetArray(coordinates, &array);
 97:     for (v = vStart; v < vEnd; ++v) {
 98:       const PetscInt idx = v - vStart;
 99:       PetscInt       val, off, d;

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

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

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

127:       if (label) {
128:         DMLabelGetValue(label, e, &val);
129:         in.segmentmarkerlist[idx] = val;
130:       }
131:     }
132:   }
133: #if 0 /* Do not currently support holes */
134:   PetscReal *holeCoords;
135:   PetscInt   h, d;

137:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
138:   if (in.numberofholes > 0) {
139:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
140:     for (h = 0; h < in.numberofholes; ++h) {
141:       for (d = 0; d < dim; ++d) {
142:         in.holelist[h*dim+d] = holeCoords[h*dim+d];
143:       }
144:     }
145:   }
146: #endif
147:   if (!rank) {
148:     char args[32];

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

163:   {
164:     DMLabel          glabel      = NULL;
165:     DMLabel          glabel2     = NULL;
166:     const PetscInt   numCorners  = 3;
167:     const PetscInt   numCells    = out.numberoftriangles;
168:     const PetscInt   numVertices = out.numberofpoints;
169:     PetscInt         *cells;
170:     PetscReal        *meshCoords;

172:     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
173:       meshCoords = (PetscReal *) out.pointlist;
174:     } else {
175:       PetscInt i;

177:       PetscMalloc1(dim * numVertices,&meshCoords);
178:       for (i = 0; i < dim * numVertices; i++) {
179:         meshCoords[i] = (PetscReal) out.pointlist[i];
180:       }
181:     }
182:     if (sizeof (PetscInt) == sizeof (out.trianglelist[0])) {
183:       cells = (PetscInt *) out.trianglelist;
184:     } else {
185:       PetscInt i;

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

220:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
221:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
222:           if (glabel)  {DMLabelSetValue(glabel,  edges[0], out.edgemarkerlist[e]);}
223:           if (glabel2) {DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);}
224:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
225:         }
226:       }
227:     }
228:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
229:   }
230: #if 0 /* Do not currently support holes */
231:   DMPlexCopyHoles(*dm, boundary);
232: #endif
233:   FiniOutput_Triangle(&out);
234:   return(0);
235: }

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

251:   PetscObjectGetComm((PetscObject)dm,&comm);
252:   MPI_Comm_rank(comm, &rank);
253:   InitInput_Triangle(&in);
254:   InitOutput_Triangle(&out);
255:   DMPlexGetDepth(dm, &depth);
256:   MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
257:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
258:   DMGetLabel(dm, labelName, &label);

260:   in.numberofpoints = vEnd - vStart;
261:   if (in.numberofpoints > 0) {
262:     PetscSection coordSection;
263:     Vec          coordinates;
264:     PetscScalar *array;

266:     PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
267:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
268:     DMGetCoordinatesLocal(dm, &coordinates);
269:     DMGetCoordinateSection(dm, &coordSection);
270:     VecGetArray(coordinates, &array);
271:     for (v = vStart; v < vEnd; ++v) {
272:       const PetscInt idx = v - vStart;
273:       PetscInt       off, d, val;

275:       PetscSectionGetOffset(coordSection, v, &off);
276:       for (d = 0; d < dim; ++d) {
277:         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
278:       }
279:       if (label) {
280:         DMLabelGetValue(label, v, &val);
281:         in.pointmarkerlist[idx] = val;
282:       }
283:     }
284:     VecRestoreArray(coordinates, &array);
285:   }
286:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
287:   DMPlexGetGhostCellStratum(dm, &gcStart, NULL);
288:   if (gcStart >= 0) cEnd = gcStart;

290:   in.numberofcorners   = 3;
291:   in.numberoftriangles = cEnd - cStart;

293: #if !defined(PETSC_USE_REAL_DOUBLE)
294:   PetscMalloc1(cEnd - cStart,&maxVolumes);
295:   for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
296: #else
297:   maxVolumes = inmaxVolumes;
298: #endif

300:   in.trianglearealist  = (double*) maxVolumes;
301:   if (in.numberoftriangles > 0) {
302:     PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
303:     for (c = cStart; c < cEnd; ++c) {
304:       const PetscInt idx      = c - cStart;
305:       PetscInt      *closure = NULL;
306:       PetscInt       closureSize;

308:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
309:       if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
310:       for (v = 0; v < 3; ++v) {
311:         in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
312:       }
313:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
314:     }
315:   }
316:   /* TODO: Segment markers are missing on input */
317: #if 0 /* Do not currently support holes */
318:   PetscReal *holeCoords;
319:   PetscInt   h, d;

321:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
322:   if (in.numberofholes > 0) {
323:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
324:     for (h = 0; h < in.numberofholes; ++h) {
325:       for (d = 0; d < dim; ++d) {
326:         in.holelist[h*dim+d] = holeCoords[h*dim+d];
327:       }
328:     }
329:   }
330: #endif
331:   if (!rank) {
332:     char args[32];

334:     /* Take away 'Q' for verbose output */
335:     PetscStrcpy(args, "pqezQra");
336:     triangulate(args, &in, &out, NULL);
337:   }
338:   PetscFree(in.pointlist);
339:   PetscFree(in.pointmarkerlist);
340:   PetscFree(in.segmentlist);
341:   PetscFree(in.segmentmarkerlist);
342:   PetscFree(in.trianglelist);

344:   {
345:     DMLabel          rlabel      = NULL;
346:     const PetscInt   numCorners  = 3;
347:     const PetscInt   numCells    = out.numberoftriangles;
348:     const PetscInt   numVertices = out.numberofpoints;
349:     PetscInt         *cells;
350:     PetscReal        *meshCoords;
351:     PetscBool        interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;

353:     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
354:       meshCoords = (PetscReal *) out.pointlist;
355:     } else {
356:       PetscInt i;

358:       PetscMalloc1(dim * numVertices,&meshCoords);
359:       for (i = 0; i < dim * numVertices; i++) {
360:         meshCoords[i] = (PetscReal) out.pointlist[i];
361:       }
362:     }
363:     if (sizeof (PetscInt) == sizeof (out.trianglelist[0])) {
364:       cells = (PetscInt *) out.trianglelist;
365:     } else {
366:       PetscInt i;

368:       PetscMalloc1(numCells * numCorners, &cells);
369:       for (i = 0; i < numCells * numCorners; i++) {
370:         cells[i] = (PetscInt) out.trianglelist[i];
371:       }
372:     }

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

394:       for (e = 0; e < out.numberofedges; e++) {
395:         if (out.edgemarkerlist[e]) {
396:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
397:           const PetscInt *edges;
398:           PetscInt        numEdges;

400:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
401:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
402:           if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
403:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
404:         }
405:       }
406:     }
407:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
408:   }
409: #if 0 /* Do not currently support holes */
410:   DMPlexCopyHoles(*dm, boundary);
411: #endif
412:   FiniOutput_Triangle(&out);
413: #if !defined(PETSC_USE_REAL_DOUBLE)
414:   PetscFree(maxVolumes);
415: #endif
416:   return(0);
417: }