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: {
     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: }