Actual source code: plexadapt.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #if defined(PETSC_HAVE_PRAGMATIC)
      3: #include <pragmatic/cpragmatic.h>
      4: #endif

      6: static PetscErrorCode DMPlexLabelToVolumeConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscReal refRatio, PetscReal maxVolumes[])
      7: {
      8:   PetscInt       dim, c;

     12:   DMGetDimension(dm, &dim);
     13:   refRatio = refRatio == PETSC_DEFAULT ? (PetscReal) ((PetscInt) 1 << dim) : refRatio;
     14:   for (c = cStart; c < cEnd; c++) {
     15:     PetscReal vol;
     16:     PetscInt  closureSize = 0, cl;
     17:     PetscInt *closure     = NULL;
     18:     PetscBool anyRefine   = PETSC_FALSE;
     19:     PetscBool anyCoarsen  = PETSC_FALSE;
     20:     PetscBool anyKeep     = PETSC_FALSE;

     22:     DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);
     23:     maxVolumes[c - cStart] = vol;
     24:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
     25:     for (cl = 0; cl < closureSize*2; cl += 2) {
     26:       const PetscInt point = closure[cl];
     27:       PetscInt       refFlag;

     29:       DMLabelGetValue(adaptLabel, point, &refFlag);
     30:       switch (refFlag) {
     31:       case DM_ADAPT_REFINE:
     32:         anyRefine  = PETSC_TRUE;break;
     33:       case DM_ADAPT_COARSEN:
     34:         anyCoarsen = PETSC_TRUE;break;
     35:       case DM_ADAPT_KEEP:
     36:         anyKeep    = PETSC_TRUE;break;
     37:       case DM_ADAPT_DETERMINE:
     38:         break;
     39:       default:
     40:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMPlex does not support refinement flag %D\n", refFlag);
     41:       }
     42:       if (anyRefine) break;
     43:     }
     44:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
     45:     if (anyRefine) {
     46:       maxVolumes[c - cStart] = vol / refRatio;
     47:     } else if (anyKeep) {
     48:       maxVolumes[c - cStart] = vol;
     49:     } else if (anyCoarsen) {
     50:       maxVolumes[c - cStart] = vol * refRatio;
     51:     }
     52:   }
     53:   return(0);
     54: }

     56: static PetscErrorCode DMPlexLabelToMetricConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscInt vStart, PetscInt vEnd, PetscReal refRatio, Vec *metricVec)
     57: {
     58:   DM              udm, coordDM;
     59:   PetscSection    coordSection;
     60:   Vec             coordinates, mb, mx;
     61:   Mat             A;
     62:   PetscScalar    *metric, *eqns;
     63:   const PetscReal coarseRatio = refRatio == PETSC_DEFAULT ? PetscSqr(0.5) : 1/refRatio;
     64:   PetscInt        dim, Nv, Neq, c, v;
     65:   PetscErrorCode  ierr;

     68:   DMPlexUninterpolate(dm, &udm);
     69:   DMGetDimension(dm, &dim);
     70:   DMGetCoordinateDM(dm, &coordDM);
     71:   DMGetLocalSection(coordDM, &coordSection);
     72:   DMGetCoordinatesLocal(dm, &coordinates);
     73:   Nv   = vEnd - vStart;
     74:   VecCreateSeq(PETSC_COMM_SELF, Nv*PetscSqr(dim), metricVec);
     75:   VecGetArray(*metricVec, &metric);
     76:   Neq  = (dim*(dim+1))/2;
     77:   PetscMalloc1(PetscSqr(Neq), &eqns);
     78:   MatCreateSeqDense(PETSC_COMM_SELF, Neq, Neq, eqns, &A);
     79:   MatCreateVecs(A, &mx, &mb);
     80:   VecSet(mb, 1.0);
     81:   for (c = cStart; c < cEnd; ++c) {
     82:     const PetscScalar *sol;
     83:     PetscScalar       *cellCoords = NULL;
     84:     PetscReal          e[3], vol;
     85:     const PetscInt    *cone;
     86:     PetscInt           coneSize, cl, i, j, d, r;

     88:     DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);
     89:     /* Only works for simplices */
     90:     for (i = 0, r = 0; i < dim+1; ++i) {
     91:       for (j = 0; j < i; ++j, ++r) {
     92:         for (d = 0; d < dim; ++d) e[d] = PetscRealPart(cellCoords[i*dim+d] - cellCoords[j*dim+d]);
     93:         /* FORTRAN ORDERING */
     94:         switch (dim) {
     95:         case 2:
     96:           eqns[0*Neq+r] = PetscSqr(e[0]);
     97:           eqns[1*Neq+r] = 2.0*e[0]*e[1];
     98:           eqns[2*Neq+r] = PetscSqr(e[1]);
     99:           break;
    100:         case 3:
    101:           eqns[0*Neq+r] = PetscSqr(e[0]);
    102:           eqns[1*Neq+r] = 2.0*e[0]*e[1];
    103:           eqns[2*Neq+r] = 2.0*e[0]*e[2];
    104:           eqns[3*Neq+r] = PetscSqr(e[1]);
    105:           eqns[4*Neq+r] = 2.0*e[1]*e[2];
    106:           eqns[5*Neq+r] = PetscSqr(e[2]);
    107:           break;
    108:         }
    109:       }
    110:     }
    111:     MatSetUnfactored(A);
    112:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);
    113:     MatLUFactor(A, NULL, NULL, NULL);
    114:     MatSolve(A, mb, mx);
    115:     VecGetArrayRead(mx, &sol);
    116:     DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);
    117:     DMPlexGetCone(udm, c, &cone);
    118:     DMPlexGetConeSize(udm, c, &coneSize);
    119:     for (cl = 0; cl < coneSize; ++cl) {
    120:       const PetscInt v = cone[cl] - vStart;

    122:       if (dim == 2) {
    123:         metric[v*4+0] += vol*coarseRatio*sol[0];
    124:         metric[v*4+1] += vol*coarseRatio*sol[1];
    125:         metric[v*4+2] += vol*coarseRatio*sol[1];
    126:         metric[v*4+3] += vol*coarseRatio*sol[2];
    127:       } else {
    128:         metric[v*9+0] += vol*coarseRatio*sol[0];
    129:         metric[v*9+1] += vol*coarseRatio*sol[1];
    130:         metric[v*9+3] += vol*coarseRatio*sol[1];
    131:         metric[v*9+2] += vol*coarseRatio*sol[2];
    132:         metric[v*9+6] += vol*coarseRatio*sol[2];
    133:         metric[v*9+4] += vol*coarseRatio*sol[3];
    134:         metric[v*9+5] += vol*coarseRatio*sol[4];
    135:         metric[v*9+7] += vol*coarseRatio*sol[4];
    136:         metric[v*9+8] += vol*coarseRatio*sol[5];
    137:       }
    138:     }
    139:     VecRestoreArrayRead(mx, &sol);
    140:   }
    141:   for (v = 0; v < Nv; ++v) {
    142:     const PetscInt *support;
    143:     PetscInt        supportSize, s;
    144:     PetscReal       vol, totVol = 0.0;

    146:     DMPlexGetSupport(udm, v+vStart, &support);
    147:     DMPlexGetSupportSize(udm, v+vStart, &supportSize);
    148:     for (s = 0; s < supportSize; ++s) {DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL); totVol += vol;}
    149:     for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol;
    150:   }
    151:   PetscFree(eqns);
    152:   VecRestoreArray(*metricVec, &metric);
    153:   VecDestroy(&mx);
    154:   VecDestroy(&mb);
    155:   MatDestroy(&A);
    156:   DMDestroy(&udm);
    157:   return(0);
    158: }

    160: /*
    161:    Contains the list of registered DMPlexGenerators routines
    162: */
    163: extern PlexGeneratorFunctionList DMPlexGenerateList;

    165: PetscErrorCode DMPlexRefine_Internal(DM dm, DMLabel adaptLabel, DM *dmRefined)
    166: {
    167:   PlexGeneratorFunctionList fl;
    168:   PetscErrorCode          (*refine)(DM,PetscReal*,DM*);
    169:   PetscErrorCode          (*adapt)(DM,DMLabel,DM*);
    170:   PetscErrorCode          (*refinementFunc)(const PetscReal [], PetscReal *);
    171:   char                      genname[PETSC_MAX_PATH_LEN], *name = NULL;
    172:   PetscReal                 refinementLimit;
    173:   PetscReal                *maxVolumes;
    174:   PetscInt                  dim, cStart, cEnd, c;
    175:   PetscBool                 flg, flg2, localized;
    176:   PetscErrorCode            ierr;

    179:   DMGetCoordinatesLocalized(dm, &localized);
    180:   DMPlexGetRefinementLimit(dm, &refinementLimit);
    181:   DMPlexGetRefinementFunction(dm, &refinementFunc);
    182:   if (refinementLimit == 0.0 && !refinementFunc && !adaptLabel) return(0);
    183:   DMGetDimension(dm, &dim);
    184:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
    185:   PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_adaptor", genname, sizeof(genname), &flg);
    186:   if (flg) name = genname;
    187:   else {
    188:     PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, sizeof(genname), &flg2);
    189:     if (flg2) name = genname;
    190:   }

    192:   fl = DMPlexGenerateList;
    193:   if (name) {
    194:     while (fl) {
    195:       PetscStrcmp(fl->name,name,&flg);
    196:       if (flg) {
    197:         refine = fl->refine;
    198:         adapt  = fl->adaptlabel;
    199:         goto gotit;
    200:       }
    201:       fl = fl->next;
    202:     }
    203:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Grid refiner %s not registered",name);
    204:   } else {
    205:     while (fl) {
    206:       if (dim-1 == fl->dim) {
    207:         refine = fl->refine;
    208:         adapt  = fl->adaptlabel;
    209:         goto gotit;
    210:       }
    211:       fl = fl->next;
    212:     }
    213:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No grid refiner of dimension %D registered",dim);
    214:   }

    216:   gotit:
    217:   switch (dim) {
    218:     case 2:
    219:     case 3:
    220:       if (adapt) {
    221:         (*adapt)(dm, adaptLabel, dmRefined);
    222:       } else {
    223:         PetscMalloc1(cEnd - cStart, &maxVolumes);
    224:         if (adaptLabel) {
    225:           DMPlexLabelToVolumeConstraint(dm, adaptLabel, cStart, cEnd, PETSC_DEFAULT, maxVolumes);
    226:         } else if (refinementFunc) {
    227:           for (c = cStart; c < cEnd; ++c) {
    228:             PetscReal vol, centroid[3];

    230:             DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
    231:             (*refinementFunc)(centroid, &maxVolumes[c-cStart]);
    232:           }
    233:         } else {
    234:           for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
    235:         }
    236:         (*refine)(dm, maxVolumes, dmRefined);
    237:         PetscFree(maxVolumes);
    238:       }
    239:       break;
    240:     default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %D is not supported.", dim);
    241:   }
    242:   DMCopyBoundary(dm, *dmRefined);
    243:   if (localized) {DMLocalizeCoordinates(*dmRefined);}
    244:   return(0);
    245: }

    247: PetscErrorCode DMPlexCoarsen_Internal(DM dm, DMLabel adaptLabel, DM *dmCoarsened)
    248: {
    249:   Vec            metricVec;
    250:   PetscInt       cStart, cEnd, vStart, vEnd;
    251:   DMLabel        bdLabel = NULL;
    252:   char           bdLabelName[PETSC_MAX_PATH_LEN];
    253:   PetscBool      localized, flg;

    257:   DMGetCoordinatesLocalized(dm, &localized);
    258:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
    259:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    260:   DMPlexLabelToMetricConstraint(dm, adaptLabel, cStart, cEnd, vStart, vEnd, PETSC_DEFAULT, &metricVec);
    261:   PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_bd_label", bdLabelName, sizeof(bdLabelName), &flg);
    262:   if (flg) {DMGetLabel(dm, bdLabelName, &bdLabel);}
    263:   DMAdaptMetric_Plex(dm, metricVec, bdLabel, dmCoarsened);
    264:   VecDestroy(&metricVec);
    265:   if (localized) {DMLocalizeCoordinates(*dmCoarsened);}
    266:   return(0);
    267: }

    269: PetscErrorCode DMAdaptLabel_Plex(DM dm, DMLabel adaptLabel, DM *dmAdapted)
    270: {
    271:   IS              flagIS;
    272:   const PetscInt *flags;
    273:   PetscInt        defFlag, minFlag, maxFlag, numFlags, f;
    274:   PetscErrorCode  ierr;

    277:   DMLabelGetDefaultValue(adaptLabel, &defFlag);
    278:   minFlag = defFlag;
    279:   maxFlag = defFlag;
    280:   DMLabelGetValueIS(adaptLabel, &flagIS);
    281:   ISGetLocalSize(flagIS, &numFlags);
    282:   ISGetIndices(flagIS, &flags);
    283:   for (f = 0; f < numFlags; ++f) {
    284:     const PetscInt flag = flags[f];

    286:     minFlag = PetscMin(minFlag, flag);
    287:     maxFlag = PetscMax(maxFlag, flag);
    288:   }
    289:   ISRestoreIndices(flagIS, &flags);
    290:   ISDestroy(&flagIS);
    291:   {
    292:     PetscInt minMaxFlag[2], minMaxFlagGlobal[2];

    294:     minMaxFlag[0] =  minFlag;
    295:     minMaxFlag[1] = -maxFlag;
    296:     MPI_Allreduce(minMaxFlag, minMaxFlagGlobal, 2, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
    297:     minFlag =  minMaxFlagGlobal[0];
    298:     maxFlag = -minMaxFlagGlobal[1];
    299:   }
    300:   if (minFlag == maxFlag) {
    301:     switch (minFlag) {
    302:     case DM_ADAPT_DETERMINE:
    303:       *dmAdapted = NULL;break;
    304:     case DM_ADAPT_REFINE:
    305:       DMPlexSetRefinementUniform(dm, PETSC_TRUE);
    306:       DMRefine(dm, MPI_COMM_NULL, dmAdapted);break;
    307:     case DM_ADAPT_COARSEN:
    308:       DMCoarsen(dm, MPI_COMM_NULL, dmAdapted);break;
    309:     default:
    310:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP,"DMPlex does not support refinement flag %D\n", minFlag);
    311:     }
    312:   } else {
    313:     DMPlexSetRefinementUniform(dm, PETSC_FALSE);
    314:     DMPlexRefine_Internal(dm, adaptLabel, dmAdapted);
    315:   }
    316:   return(0);
    317: }

    319: /*
    320:   DMAdaptMetric_Plex - Generates a new mesh conforming to a metric field.

    322:   Input Parameters:
    323: + dm - The DM object
    324: . vertexMetric - The metric to which the mesh is adapted, defined vertex-wise in a LOCAL vector
    325: - bdLabel - Label for boundary tags which are preserved in dmNew, or NULL. Should not be named "_boundary_".

    327:   Output Parameter:
    328: . dmNew  - the new DM

    330:   Level: advanced

    332: .seealso: DMCoarsen(), DMRefine()
    333: */
    334: PetscErrorCode DMAdaptMetric_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew)
    335: {
    336: #if defined(PETSC_HAVE_PRAGMATIC)
    337:   MPI_Comm           comm;
    338:   const char        *bdName = "_boundary_";
    339: #if 0
    340:   DM                 odm = dm;
    341: #endif
    342:   DM                 udm, cdm;
    343:   DMLabel            bdLabelFull;
    344:   const char        *bdLabelName;
    345:   IS                 bdIS, globalVertexNum;
    346:   PetscSection       coordSection;
    347:   Vec                coordinates;
    348:   const PetscScalar *coords, *met;
    349:   const PetscInt    *bdFacesFull, *gV;
    350:   PetscInt          *bdFaces, *bdFaceIds, *l2gv;
    351:   PetscReal         *x, *y, *z, *metric;
    352:   PetscInt          *cells;
    353:   PetscInt           dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v;
    354:   PetscInt           off, maxConeSize, numBdFaces, f, bdSize;
    355:   PetscBool          flg;
    356:   DMLabel            bdLabelNew;
    357:   PetscReal         *coordsNew;
    358:   PetscInt          *bdTags;
    359:   PetscReal         *xNew[3] = {NULL, NULL, NULL};
    360:   PetscInt          *cellsNew;
    361:   PetscInt           d, numCellsNew, numVerticesNew;
    362:   PetscInt           numCornersNew, fStart, fEnd;
    363:   PetscMPIInt        numProcs;
    364:   PetscErrorCode     ierr;

    367:   /* Check for FEM adjacency flags */
    368:   PetscObjectGetComm((PetscObject) dm, &comm);
    369:   MPI_Comm_size(comm, &numProcs);
    370:   if (bdLabel) {
    371:     PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);
    372:     PetscStrcmp(bdLabelName, bdName, &flg);
    373:     if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
    374:   }
    375:   /* Add overlap for Pragmatic */
    376: #if 0
    377:   /* Check for overlap by looking for cell in the SF */
    378:   if (!overlapped) {
    379:     DMPlexDistributeOverlap(odm, 1, NULL, &dm);
    380:     if (!dm) {dm = odm; PetscObjectReference((PetscObject) dm);}
    381:   }
    382: #endif
    383:   /* Get mesh information */
    384:   DMGetDimension(dm, &dim);
    385:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
    386:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    387:   DMPlexUninterpolate(dm, &udm);
    388:   DMPlexGetMaxSizes(udm, &maxConeSize, NULL);
    389:   numCells    = cEnd - cStart;
    390:   numVertices = vEnd - vStart;
    391:   PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);
    392:   for (c = 0, coff = 0; c < numCells; ++c) {
    393:     const PetscInt *cone;
    394:     PetscInt        coneSize, cl;

    396:     DMPlexGetConeSize(udm, c, &coneSize);
    397:     DMPlexGetCone(udm, c, &cone);
    398:     for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
    399:   }
    400:   PetscCalloc1(numVertices, &l2gv);
    401:   DMPlexGetVertexNumbering(udm, &globalVertexNum);
    402:   ISGetIndices(globalVertexNum, &gV);
    403:   for (v = 0, numLocVertices = 0; v < numVertices; ++v) {
    404:     if (gV[v] >= 0) ++numLocVertices;
    405:     l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v];
    406:   }
    407:   ISRestoreIndices(globalVertexNum, &gV);
    408:   DMDestroy(&udm);
    409:   DMGetCoordinateDM(dm, &cdm);
    410:   DMGetLocalSection(cdm, &coordSection);
    411:   DMGetCoordinatesLocal(dm, &coordinates);
    412:   VecGetArrayRead(coordinates, &coords);
    413:   for (v = vStart; v < vEnd; ++v) {
    414:     PetscSectionGetOffset(coordSection, v, &off);
    415:     x[v-vStart] = PetscRealPart(coords[off+0]);
    416:     if (dim > 1) y[v-vStart] = PetscRealPart(coords[off+1]);
    417:     if (dim > 2) z[v-vStart] = PetscRealPart(coords[off+2]);
    418:   }
    419:   VecRestoreArrayRead(coordinates, &coords);
    420:   /* Get boundary mesh */
    421:   DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);
    422:   DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);
    423:   DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);
    424:   DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);
    425:   ISGetIndices(bdIS, &bdFacesFull);
    426:   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
    427:     PetscInt *closure = NULL;
    428:     PetscInt  closureSize, cl;

    430:     DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
    431:     for (cl = 0; cl < closureSize*2; cl += 2) {
    432:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
    433:     }
    434:     DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
    435:   }
    436:   PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);
    437:   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
    438:     PetscInt *closure = NULL;
    439:     PetscInt  closureSize, cl;

    441:     DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
    442:     for (cl = 0; cl < closureSize*2; cl += 2) {
    443:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart;
    444:     }
    445:     DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
    446:     if (bdLabel) {DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);}
    447:     else         {bdFaceIds[f] = 1;}
    448:   }
    449:   ISDestroy(&bdIS);
    450:   DMLabelDestroy(&bdLabelFull);
    451:   /* Get metric */
    452:   VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");
    453:   VecGetArrayRead(vertexMetric, &met);
    454:   for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]);
    455:   VecRestoreArrayRead(vertexMetric, &met);
    456: #if 0
    457:   /* Destroy overlap mesh */
    458:   DMDestroy(&dm);
    459: #endif
    460:   /* Create new mesh */
    461:   switch (dim) {
    462:   case 2:
    463:     pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);break;
    464:   case 3:
    465:     pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);break;
    466:   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim);
    467:   }
    468:   pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
    469:   pragmatic_set_metric(metric);
    470:   pragmatic_adapt(((DM_Plex *) dm->data)->remeshBd ? 1 : 0);
    471:   PetscFree(l2gv);
    472:   /* Read out mesh */
    473:   pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew);
    474:   PetscMalloc1(numVerticesNew*dim, &coordsNew);
    475:   switch (dim) {
    476:   case 2:
    477:     numCornersNew = 3;
    478:     PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);
    479:     pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]);
    480:     break;
    481:   case 3:
    482:     numCornersNew = 4;
    483:     PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);
    484:     pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]);
    485:     break;
    486:   default:
    487:     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim);
    488:   }
    489:   for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];}
    490:   PetscMalloc1(numCellsNew*(dim+1), &cellsNew);
    491:   pragmatic_get_elements(cellsNew);
    492:   DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, dmNew);
    493:   /* Read out boundary label */
    494:   pragmatic_get_boundaryTags(&bdTags);
    495:   DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);
    496:   DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);
    497:   DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);
    498:   DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);
    499:   DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);
    500:   for (c = cStart; c < cEnd; ++c) {
    501:     /* Only for simplicial meshes */
    502:     coff = (c-cStart)*(dim+1);
    503:     /* d is the local cell number of the vertex opposite to the face we are marking */
    504:     for (d = 0; d < dim+1; ++d) {
    505:       if (bdTags[coff+d]) {
    506:         const PetscInt  perm[4][4] = {{-1, -1, -1, -1}, {-1, -1, -1, -1}, {1, 2, 0, -1}, {3, 2, 1, 0}}; /* perm[d] = face opposite */
    507:         const PetscInt *cone;

    509:         /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
    510:         DMPlexGetCone(*dmNew, c, &cone);
    511:         DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);
    512:       }
    513:     }
    514:   }
    515:   /* Cleanup */
    516:   switch (dim) {
    517:   case 2: PetscFree2(xNew[0], xNew[1]);break;
    518:   case 3: PetscFree3(xNew[0], xNew[1], xNew[2]);break;
    519:   }
    520:   PetscFree(cellsNew);
    521:   PetscFree5(x, y, z, metric, cells);
    522:   PetscFree2(bdFaces, bdFaceIds);
    523:   PetscFree(coordsNew);
    524:   pragmatic_finalize();
    525:   return(0);
    526: #else
    528:   SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic.");
    529: #endif
    530: }