Actual source code: trigenerate.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>

  3: #include <triangle.h>

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

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

 44: static PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx)
 45: {
 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;
 72:   PetscErrorCode       ierr;

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

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

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

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

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

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

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

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

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

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

169:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
170:     if (label)  {DMCreateLabel(*dm, labelName);  DMGetLabel(*dm, labelName,  &glabel);}
171:     if (label2) {DMCreateLabel(*dm, labelName2); DMGetLabel(*dm, labelName2, &glabel2);}
172:     /* Set labels */
173:     for (v = 0; v < numVertices; ++v) {
174:       if (out.pointmarkerlist[v]) {
175:         if (glabel) {DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);}
176:       }
177:     }
178:     if (interpolate) {
179:       for (e = 0; e < out.numberofedges; e++) {
180:         if (out.edgemarkerlist[e]) {
181:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
182:           const PetscInt *edges;
183:           PetscInt        numEdges;

185:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
186:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
187:           if (glabel)  {DMLabelSetValue(glabel,  edges[0], out.edgemarkerlist[e]);}
188:           if (glabel2) {DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);}
189:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
190:         }
191:       }
192:     }
193:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
194:   }
195: #if 0 /* Do not currently support holes */
196:   DMPlexCopyHoles(*dm, boundary);
197: #endif
198:   FiniOutput_Triangle(&out);
199:   return(0);
200: }

202: PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined)
203: {
204:   MPI_Comm             comm;
205:   PetscInt             dim       = 2;
206:   const char          *labelName = "marker";
207:   struct triangulateio in;
208:   struct triangulateio out;
209:   DMLabel              label;
210:   PetscInt             vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
211:   PetscMPIInt          rank;
212:   PetscErrorCode       ierr;
213:   double               *maxVolumes;

216:   PetscObjectGetComm((PetscObject)dm,&comm);
217:   MPI_Comm_rank(comm, &rank);
218:   InitInput_Triangle(&in);
219:   InitOutput_Triangle(&out);
220:   DMPlexGetDepth(dm, &depth);
221:   MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
222:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
223:   DMGetLabel(dm, labelName, &label);

225:   in.numberofpoints = vEnd - vStart;
226:   if (in.numberofpoints > 0) {
227:     PetscSection coordSection;
228:     Vec          coordinates;
229:     PetscScalar *array;

231:     PetscMalloc1(in.numberofpoints*dim, &in.pointlist);
232:     PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);
233:     DMGetCoordinatesLocal(dm, &coordinates);
234:     DMGetCoordinateSection(dm, &coordSection);
235:     VecGetArray(coordinates, &array);
236:     for (v = vStart; v < vEnd; ++v) {
237:       const PetscInt idx = v - vStart;
238:       PetscInt       off, d, val;

240:       PetscSectionGetOffset(coordSection, v, &off);
241:       for (d = 0; d < dim; ++d) {
242:         in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
243:       }
244:       if (label) {
245:         DMLabelGetValue(label, v, &val);
246:         in.pointmarkerlist[idx] = val;
247:       }
248:     }
249:     VecRestoreArray(coordinates, &array);
250:   }
251:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

253:   in.numberofcorners   = 3;
254:   in.numberoftriangles = cEnd - cStart;

256: #if !defined(PETSC_USE_REAL_DOUBLE)
257:   PetscMalloc1(cEnd - cStart,&maxVolumes);
258:   for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c];
259: #else
260:   maxVolumes = inmaxVolumes;
261: #endif

263:   in.trianglearealist  = (double*) maxVolumes;
264:   if (in.numberoftriangles > 0) {
265:     PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);
266:     for (c = cStart; c < cEnd; ++c) {
267:       const PetscInt idx      = c - cStart;
268:       PetscInt      *closure = NULL;
269:       PetscInt       closureSize;

271:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
272:       if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
273:       for (v = 0; v < 3; ++v) {
274:         in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
275:       }
276:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
277:     }
278:   }
279:   /* TODO: Segment markers are missing on input */
280: #if 0 /* Do not currently support holes */
281:   PetscReal *holeCoords;
282:   PetscInt   h, d;

284:   DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);
285:   if (in.numberofholes > 0) {
286:     PetscMalloc1(in.numberofholes*dim, &in.holelist);
287:     for (h = 0; h < in.numberofholes; ++h) {
288:       for (d = 0; d < dim; ++d) {
289:         in.holelist[h*dim+d] = holeCoords[h*dim+d];
290:       }
291:     }
292:   }
293: #endif
294:   if (!rank) {
295:     char args[32];

297:     /* Take away 'Q' for verbose output */
298:     PetscStrcpy(args, "pqezQra");
299:     triangulate(args, &in, &out, NULL);
300:   }
301:   PetscFree(in.pointlist);
302:   PetscFree(in.pointmarkerlist);
303:   PetscFree(in.segmentlist);
304:   PetscFree(in.segmentmarkerlist);
305:   PetscFree(in.trianglelist);

307:   {
308:     DMLabel        rlabel      = NULL;
309:     const PetscInt numCorners  = 3;
310:     const PetscInt numCells    = out.numberoftriangles;
311:     const PetscInt numVertices = out.numberofpoints;
312:     const int     *cells      = out.trianglelist;
313:     const double  *meshCoords = out.pointlist;
314:     PetscBool      interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;

316:     DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
317:     if (label) {DMCreateLabel(*dmRefined, labelName); DMGetLabel(*dmRefined, labelName, &rlabel);}
318:     /* Set labels */
319:     for (v = 0; v < numVertices; ++v) {
320:       if (out.pointmarkerlist[v]) {
321:         if (rlabel) {DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);}
322:       }
323:     }
324:     if (interpolate) {
325:       PetscInt e;

327:       for (e = 0; e < out.numberofedges; e++) {
328:         if (out.edgemarkerlist[e]) {
329:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
330:           const PetscInt *edges;
331:           PetscInt        numEdges;

333:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
334:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
335:           if (rlabel) {DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);}
336:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
337:         }
338:       }
339:     }
340:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
341:   }
342: #if 0 /* Do not currently support holes */
343:   DMPlexCopyHoles(*dm, boundary);
344: #endif
345:   FiniOutput_Triangle(&out);
346: #if !defined(PETSC_USE_REAL_DOUBLE)
347:   PetscFree(maxVolumes);
348: #endif
349:   return(0);
350: }