Actual source code: parmmgadapt.c

      1: #include "../mmgcommon.h" /*I      "petscdmplex.h"   I*/
      2: #include <parmmg/libparmmg.h>

      4: PetscBool  ParMmgCite       = PETSC_FALSE;
      5: const char ParMmgCitation[] = "@techreport{cirrottola:hal-02386837,\n"
      6:                               "  title       = {Parallel unstructured mesh adaptation using iterative remeshing and repartitioning},\n"
      7:                               "  institution = {Inria Bordeaux},\n"
      8:                               "  author      = {L. Cirrottola and A. Froehly},\n"
      9:                               "  number      = {9307},\n"
     10:                               "  note        = {\\url{https://hal.inria.fr/hal-02386837}},\n"
     11:                               "  year        = {2019}\n}\n";

     13: /*
     14:  Coupling code for the ParMmg metric-based mesh adaptation package.
     15: */
     16: PETSC_EXTERN PetscErrorCode DMAdaptMetric_ParMmg_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DMLabel rgLabel, DM *dmNew)
     17: {
     18:   MPI_Comm           comm;
     19:   const char        *bdName = "_boundary_";
     20:   const char        *rgName = "_regions_";
     21:   DM                 udm, cdm;
     22:   DMLabel            bdLabelNew, rgLabelNew;
     23:   const char        *bdLabelName, *rgLabelName;
     24:   IS                 globalVertexNum;
     25:   PetscSection       coordSection;
     26:   Vec                coordinates;
     27:   PetscSF            sf;
     28:   const PetscScalar *coords, *met;
     29:   PetscReal         *vertices, *metric, *verticesNew, *verticesNewLoc, gradationFactor, hausdorffNumber;
     30:   PetscInt          *cells, *cellsNew, *cellTags, *cellTagsNew, *verTags, *verTagsNew;
     31:   PetscInt          *bdFaces, *faceTags, *facesNew, *faceTagsNew;
     32:   PetscInt          *corners, *requiredCells, *requiredVer, *ridges, *requiredFaces;
     33:   PetscInt           cStart, cEnd, c, numCells, fStart, fEnd, f, numFaceTags, vStart, vEnd, v, numVertices;
     34:   PetscInt           numCellsNotShared, *cIsLeaf, numUsedVertices, *vertexNumber, *fIsIncluded;
     35:   PetscInt           dim, off, coff, maxConeSize, bdSize, i, j, k, Neq, verbosity, numIter;
     36:   PetscInt          *interfaces_lv, *interfaces_gv, *interfacesOffset;
     37:   PetscMPIInt        niranks, nrranks, numNgbRanks;
     38:   PetscInt           r, lv, gv;
     39:   PetscInt          *gv_new, *owners, *verticesNewSorted, pStart, pEnd;
     40:   PetscInt           numCellsNew, numVerticesNew, numCornersNew, numFacesNew, numVerticesNewLoc;
     41:   const PetscInt    *gV, *ioffset, *irootloc, *roffset, *rmine, *rremote;
     42:   PetscBool          flg = PETSC_FALSE, noInsert, noSwap, noMove, noSurf, isotropic, uniform;
     43:   const PetscMPIInt *iranks, *rranks;
     44:   PetscMPIInt        numProcs, rank;
     45:   PMMG_pParMesh      parmesh = NULL;

     47:   // DEVELOPER NOTE: ParMmg wants to know the rank of every process which is sharing a given point and
     48:   //                 for this information to be conveyed to every process that is sharing that point.

     50:   PetscFunctionBegin;
     51:   PetscCall(PetscCitationsRegister(ParMmgCitation, &ParMmgCite));
     52:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
     53:   PetscCallMPI(MPI_Comm_size(comm, &numProcs));
     54:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
     55:   if (bdLabel) {
     56:     PetscCall(PetscObjectGetName((PetscObject)bdLabel, &bdLabelName));
     57:     PetscCall(PetscStrcmp(bdLabelName, bdName, &flg));
     58:     PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
     59:   }
     60:   if (rgLabel) {
     61:     PetscCall(PetscObjectGetName((PetscObject)rgLabel, &rgLabelName));
     62:     PetscCall(PetscStrcmp(rgLabelName, rgName, &flg));
     63:     PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for element tags", rgLabelName);
     64:   }

     66:   /* Get mesh information */
     67:   PetscCall(DMGetDimension(dm, &dim));
     68:   PetscCheck(dim == 3, comm, PETSC_ERR_ARG_OUTOFRANGE, "ParMmg only works in 3D.");
     69:   Neq = (dim * (dim + 1)) / 2;
     70:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
     71:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
     72:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
     73:   PetscCall(DMPlexUninterpolate(dm, &udm));
     74:   PetscCall(DMPlexGetMaxSizes(udm, &maxConeSize, NULL));
     75:   numCells    = cEnd - cStart;
     76:   numVertices = vEnd - vStart;

     78:   /* Get parallel data; work out which cells are owned and which are leaves */
     79:   PetscCall(PetscCalloc1(numCells, &cIsLeaf));
     80:   numCellsNotShared = numCells;
     81:   niranks = nrranks = 0;
     82:   if (numProcs > 1) {
     83:     PetscCall(DMGetPointSF(dm, &sf));
     84:     PetscCall(PetscSFSetUp(sf));
     85:     PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, &ioffset, &irootloc));
     86:     PetscCall(PetscSFGetRootRanks(sf, &nrranks, &rranks, &roffset, &rmine, &rremote));
     87:     for (r = 0; r < nrranks; ++r) {
     88:       for (i = roffset[r]; i < roffset[r + 1]; ++i) {
     89:         if (rmine[i] >= cStart && rmine[i] < cEnd) {
     90:           cIsLeaf[rmine[i] - cStart] = 1;
     91:           numCellsNotShared--;
     92:         }
     93:       }
     94:     }
     95:   }

     97:   /*
     98:     Create a vertex numbering for ParMmg starting at 1. Vertices not included in any
     99:     owned cell remain 0 and will be removed. Using this numbering, create cells.
    100:   */
    101:   numUsedVertices = 0;
    102:   PetscCall(PetscCalloc1(numVertices, &vertexNumber));
    103:   PetscCall(PetscMalloc1(numCellsNotShared * maxConeSize, &cells));
    104:   for (c = 0, coff = 0; c < numCells; ++c) {
    105:     const PetscInt *cone;
    106:     PetscInt        coneSize, cl;

    108:     if (!cIsLeaf[c]) {
    109:       PetscCall(DMPlexGetConeSize(udm, cStart + c, &coneSize));
    110:       PetscCall(DMPlexGetCone(udm, cStart + c, &cone));
    111:       for (cl = 0; cl < coneSize; ++cl) {
    112:         if (!vertexNumber[cone[cl] - vStart]) vertexNumber[cone[cl] - vStart] = ++numUsedVertices;
    113:         cells[coff++] = vertexNumber[cone[cl] - vStart];
    114:       }
    115:     }
    116:   }

    118:   /* Get array of vertex coordinates */
    119:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    120:   PetscCall(DMGetLocalSection(cdm, &coordSection));
    121:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
    122:   PetscCall(VecGetArrayRead(coordinates, &coords));
    123:   PetscCall(PetscMalloc2(numUsedVertices * Neq, &metric, dim * numUsedVertices, &vertices));
    124:   for (v = 0; v < vEnd - vStart; ++v) {
    125:     PetscCall(PetscSectionGetOffset(coordSection, v + vStart, &off));
    126:     if (vertexNumber[v]) {
    127:       for (i = 0; i < dim; ++i) vertices[dim * (vertexNumber[v] - 1) + i] = PetscRealPart(coords[off + i]);
    128:     }
    129:   }
    130:   PetscCall(VecRestoreArrayRead(coordinates, &coords));

    132:   /* Get face tags */
    133:   if (!bdLabel) {
    134:     flg = PETSC_TRUE;
    135:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabel));
    136:     PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabel));
    137:   }
    138:   PetscCall(DMLabelGetBounds(bdLabel, &pStart, &pEnd));
    139:   PetscCall(PetscCalloc1(pEnd - pStart, &fIsIncluded));
    140:   for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
    141:     PetscBool hasPoint;
    142:     PetscInt *closure = NULL, closureSize, cl;

    144:     PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint));
    145:     if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;

    147:     /* Only faces adjacent to an owned (non-leaf) cell are included */
    148:     PetscInt        nnbrs;
    149:     const PetscInt *nbrs;
    150:     PetscCall(DMPlexGetSupportSize(dm, f, &nnbrs));
    151:     PetscCall(DMPlexGetSupport(dm, f, &nbrs));
    152:     for (c = 0; c < nnbrs; ++c) fIsIncluded[f - pStart] = fIsIncluded[f - pStart] || !cIsLeaf[nbrs[c]];
    153:     if (!fIsIncluded[f - pStart]) continue;

    155:     numFaceTags++;

    157:     PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
    158:     for (cl = 0; cl < closureSize * 2; cl += 2) {
    159:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
    160:     }
    161:     PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
    162:   }
    163:   PetscCall(PetscMalloc2(bdSize, &bdFaces, numFaceTags, &faceTags));
    164:   for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
    165:     PetscInt *closure = NULL, closureSize, cl;

    167:     if (!fIsIncluded[f - pStart]) continue;

    169:     PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
    170:     for (cl = 0; cl < closureSize * 2; cl += 2) {
    171:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = vertexNumber[closure[cl] - vStart];
    172:     }
    173:     PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
    174:     PetscCall(DMLabelGetValue(bdLabel, f, &faceTags[numFaceTags++]));
    175:   }
    176:   PetscCall(PetscFree(fIsIncluded));

    178:   /* Get cell tags */
    179:   PetscCall(PetscCalloc2(numUsedVertices, &verTags, numCellsNotShared, &cellTags));
    180:   if (rgLabel) {
    181:     for (c = cStart, coff = 0; c < cEnd; ++c) {
    182:       if (!cIsLeaf[c - cStart]) PetscCall(DMLabelGetValue(rgLabel, c, &cellTags[coff++]));
    183:     }
    184:   }
    185:   PetscCall(PetscFree(cIsLeaf));

    187:   /* Get metric, using only the upper triangular part */
    188:   PetscCall(VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view"));
    189:   PetscCall(VecGetArrayRead(vertexMetric, &met));
    190:   PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
    191:   PetscCall(DMPlexMetricIsUniform(dm, &uniform));
    192:   for (v = 0; v < (vEnd - vStart); ++v) {
    193:     PetscInt vv = vertexNumber[v];
    194:     if (!vv--) continue;
    195:     for (i = 0, k = 0; i < dim; ++i) {
    196:       for (j = i; j < dim; ++j, ++k) {
    197:         if (isotropic) {
    198:           if (i == j) {
    199:             if (uniform) metric[Neq * vv + k] = PetscRealPart(met[0]);
    200:             else metric[Neq * vv + k] = PetscRealPart(met[v]);
    201:           } else metric[Neq * vv + k] = 0.0;
    202:         } else metric[Neq * vv + k] = PetscRealPart(met[dim * dim * v + dim * i + j]);
    203:       }
    204:     }
    205:   }
    206:   PetscCall(VecRestoreArrayRead(vertexMetric, &met));

    208:   /* Build ParMmg communicators: the list of vertices between two partitions  */
    209:   numNgbRanks = 0;
    210:   if (numProcs > 1) {
    211:     DM              rankdm;
    212:     PetscSection    rankSection, rankGlobalSection;
    213:     PetscSF         rankSF;
    214:     const PetscInt *degree;
    215:     PetscInt       *rankOfUsedVertices, *rankOfUsedMultiRootLeaves, *usedCopies;
    216:     PetscInt       *rankArray, *rankGlobalArray, *interfacesPerRank;
    217:     PetscInt        offset, mrl, rootDegreeCnt, s, shareCnt, gv;

    219:     PetscCall(PetscSFComputeDegreeBegin(sf, &degree));
    220:     PetscCall(PetscSFComputeDegreeEnd(sf, &degree));
    221:     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    222:     for (i = 0, rootDegreeCnt = 0; i < pEnd - pStart; ++i) rootDegreeCnt += degree[i];

    224:     /* rankOfUsedVertices, point-array: rank+1 if vertex and in use */
    225:     PetscCall(PetscCalloc1(pEnd - pStart, &rankOfUsedVertices));
    226:     for (i = 0; i < pEnd - pStart; ++i) rankOfUsedVertices[i] = -1;
    227:     for (i = vStart; i < vEnd; ++i) {
    228:       if (vertexNumber[i - vStart]) rankOfUsedVertices[i] = rank;
    229:     }

    231:     /* rankOfUsedMultiRootLeaves, multiroot-array: rank if vertex and in use, else -1 */
    232:     PetscCall(PetscMalloc1(rootDegreeCnt, &rankOfUsedMultiRootLeaves));
    233:     PetscCall(PetscSFGatherBegin(sf, MPIU_INT, rankOfUsedVertices, rankOfUsedMultiRootLeaves));
    234:     PetscCall(PetscSFGatherEnd(sf, MPIU_INT, rankOfUsedVertices, rankOfUsedMultiRootLeaves));
    235:     PetscCall(PetscFree(rankOfUsedVertices));

    237:     /* usedCopies, point-array: if vertex, shared by how many processes */
    238:     PetscCall(PetscCalloc1(pEnd - pStart, &usedCopies));
    239:     for (i = 0, mrl = 0; i < vStart - pStart; i++) mrl += degree[i];
    240:     for (i = vStart - pStart; i < vEnd - pStart; ++i) {
    241:       for (j = 0; j < degree[i]; j++, mrl++) {
    242:         if (rankOfUsedMultiRootLeaves[mrl] != -1) usedCopies[i]++;
    243:       }
    244:       if (vertexNumber[i - vStart + pStart]) usedCopies[i]++;
    245:     }
    246:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, usedCopies, usedCopies, MPI_REPLACE));
    247:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, usedCopies, usedCopies, MPI_REPLACE));

    249:     /* Create a section to store ranks of vertices shared by more than one process */
    250:     PetscCall(PetscSectionCreate(comm, &rankSection));
    251:     PetscCall(PetscSectionSetNumFields(rankSection, 1));
    252:     PetscCall(PetscSectionSetChart(rankSection, pStart, pEnd));
    253:     for (i = vStart - pStart; i < vEnd - pStart; ++i) {
    254:       if (usedCopies[i] > 1) PetscCall(PetscSectionSetDof(rankSection, i + pStart, usedCopies[i]));
    255:     }
    256:     PetscCall(PetscSectionSetUp(rankSection));
    257:     PetscCall(PetscSectionCreateGlobalSection(rankSection, sf, PETSC_FALSE, PETSC_FALSE, PETSC_TRUE, &rankGlobalSection));

    259:     PetscCall(PetscSectionGetStorageSize(rankGlobalSection, &s));
    260:     PetscCall(PetscMalloc1(s, &rankGlobalArray));
    261:     for (i = 0, mrl = 0; i < vStart - pStart; i++) mrl += degree[i];
    262:     for (i = vStart - pStart, k = 0; i < vEnd - pStart; ++i) {
    263:       if (usedCopies[i] > 1 && degree[i]) {
    264:         PetscCall(PetscSectionGetOffset(rankSection, k, &offset));
    265:         if (vertexNumber[i - vStart + pStart]) rankGlobalArray[k++] = rank;
    266:         for (j = 0; j < degree[i]; j++, mrl++) {
    267:           if (rankOfUsedMultiRootLeaves[mrl] != -1) rankGlobalArray[k++] = rankOfUsedMultiRootLeaves[mrl];
    268:         }
    269:       } else mrl += degree[i];
    270:     }
    271:     PetscCall(PetscFree(rankOfUsedMultiRootLeaves));
    272:     PetscCall(PetscFree(usedCopies));
    273:     PetscCall(PetscSectionDestroy(&rankGlobalSection));

    275:     /*
    276:       Broadcast the array of ranks.
    277:         (We want all processes to know all the ranks that are looking at each point.
    278:         Above, we tell the roots. Here, the roots tell the leaves.)
    279:     */
    280:     PetscCall(DMClone(dm, &rankdm));
    281:     PetscCall(DMSetLocalSection(rankdm, rankSection));
    282:     PetscCall(DMGetSectionSF(rankdm, &rankSF));
    283:     PetscCall(PetscSectionGetStorageSize(rankSection, &s));
    284:     PetscCall(PetscMalloc1(s, &rankArray));
    285:     PetscCall(PetscSFBcastBegin(rankSF, MPI_INT, rankGlobalArray, rankArray, MPI_REPLACE));
    286:     PetscCall(PetscSFBcastEnd(rankSF, MPI_INT, rankGlobalArray, rankArray, MPI_REPLACE));
    287:     PetscCall(PetscFree(rankGlobalArray));
    288:     PetscCall(DMDestroy(&rankdm));

    290:     /* Count the number of interfaces per rank, not including those on the root */
    291:     PetscCall(PetscCalloc1(numProcs, &interfacesPerRank));
    292:     for (v = vStart; v < vEnd; v++) {
    293:       if (vertexNumber[v - vStart]) {
    294:         PetscCall(PetscSectionGetDof(rankSection, v, &shareCnt));
    295:         if (shareCnt) {
    296:           PetscCall(PetscSectionGetOffset(rankSection, v, &offset));
    297:           for (j = 0; j < shareCnt; j++) interfacesPerRank[rankArray[offset + j]]++;
    298:         }
    299:       }
    300:     }
    301:     for (r = 0, k = 0, interfacesPerRank[rank] = 0; r < numProcs; r++) k += interfacesPerRank[r];

    303:     /* Get the degree of the vertex */
    304:     PetscCall(PetscMalloc3(k, &interfaces_lv, k, &interfaces_gv, numProcs + 1, &interfacesOffset));
    305:     interfacesOffset[0] = 0;
    306:     for (r = 0; r < numProcs; r++) {
    307:       interfacesOffset[r + 1] = interfacesOffset[r] + interfacesPerRank[r];
    308:       if (interfacesPerRank[r]) numNgbRanks++;
    309:       interfacesPerRank[r] = 0;
    310:     }

    312:     /* Get the local and global vertex numbers at interfaces */
    313:     PetscCall(DMPlexGetVertexNumbering(dm, &globalVertexNum));
    314:     PetscCall(ISGetIndices(globalVertexNum, &gV));
    315:     for (v = vStart; v < vEnd; v++) {
    316:       if (vertexNumber[v - vStart]) {
    317:         PetscCall(PetscSectionGetDof(rankSection, v, &shareCnt));
    318:         if (shareCnt) {
    319:           PetscCall(PetscSectionGetOffset(rankSection, v, &offset));
    320:           for (j = 0; j < shareCnt; j++) {
    321:             r = rankArray[offset + j];
    322:             if (r == rank) continue;
    323:             k                = interfacesOffset[r] + interfacesPerRank[r]++;
    324:             interfaces_lv[k] = vertexNumber[v - vStart];
    325:             gv               = gV[v - vStart];
    326:             interfaces_gv[k] = gv < 0 ? -gv : gv + 1;
    327:           }
    328:         }
    329:       }
    330:     }
    331:     PetscCall(PetscFree(interfacesPerRank));
    332:     PetscCall(PetscFree(rankArray));
    333:     PetscCall(ISRestoreIndices(globalVertexNum, &gV));
    334:     PetscCall(PetscSectionDestroy(&rankSection));
    335:   }
    336:   PetscCall(DMDestroy(&udm));
    337:   PetscCall(PetscFree(vertexNumber));

    339:   /* Send the data to ParMmg and remesh */
    340:   PetscCall(DMPlexMetricNoInsertion(dm, &noInsert));
    341:   PetscCall(DMPlexMetricNoSwapping(dm, &noSwap));
    342:   PetscCall(DMPlexMetricNoMovement(dm, &noMove));
    343:   PetscCall(DMPlexMetricNoSurf(dm, &noSurf));
    344:   PetscCall(DMPlexMetricGetVerbosity(dm, &verbosity));
    345:   PetscCall(DMPlexMetricGetNumIterations(dm, &numIter));
    346:   PetscCall(DMPlexMetricGetGradationFactor(dm, &gradationFactor));
    347:   PetscCall(DMPlexMetricGetHausdorffNumber(dm, &hausdorffNumber));
    348:   PetscCallMMG_NONSTANDARD(PMMG_Init_parMesh, PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_pMesh, PMMG_ARG_pMet, PMMG_ARG_dim, 3, PMMG_ARG_MPIComm, comm, PMMG_ARG_end);
    349:   PetscCallMMG_NONSTANDARD(PMMG_Set_meshSize, parmesh, numUsedVertices, numCellsNotShared, 0, numFaceTags, 0, 0);
    350:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_APImode, PMMG_APIDISTRIB_nodes);
    351:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_noinsert, noInsert);
    352:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_noswap, noSwap);
    353:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_nomove, noMove);
    354:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_nosurf, noSurf);
    355:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_verbose, verbosity);
    356:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_globalNum, 1);
    357:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_niter, numIter);
    358:   PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter, parmesh, PMMG_DPARAM_hgrad, gradationFactor);
    359:   PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter, parmesh, PMMG_DPARAM_hausd, hausdorffNumber);
    360:   PetscCallMMG_NONSTANDARD(PMMG_Set_vertices, parmesh, vertices, verTags);
    361:   PetscCallMMG_NONSTANDARD(PMMG_Set_tetrahedra, parmesh, cells, cellTags);
    362:   PetscCallMMG_NONSTANDARD(PMMG_Set_triangles, parmesh, bdFaces, faceTags);
    363:   PetscCallMMG_NONSTANDARD(PMMG_Set_metSize, parmesh, MMG5_Vertex, numUsedVertices, MMG5_Tensor);
    364:   PetscCallMMG_NONSTANDARD(PMMG_Set_tensorMets, parmesh, metric);
    365:   PetscCallMMG_NONSTANDARD(PMMG_Set_numberOfNodeCommunicators, parmesh, numNgbRanks);
    366:   for (r = 0, c = 0; r < numProcs; ++r) {
    367:     if (interfacesOffset[r + 1] > interfacesOffset[r]) {
    368:       PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicatorSize, parmesh, c, r, interfacesOffset[r + 1] - interfacesOffset[r]);
    369:       PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicator_nodes, parmesh, c++, &interfaces_lv[interfacesOffset[r]], &interfaces_gv[interfacesOffset[r]], 1);
    370:     }
    371:   }
    372:   PetscCallMMG(PMMG_parmmglib_distributed, parmesh);
    373:   PetscCall(PetscFree(cells));
    374:   PetscCall(PetscFree2(metric, vertices));
    375:   PetscCall(PetscFree2(bdFaces, faceTags));
    376:   PetscCall(PetscFree2(verTags, cellTags));
    377:   if (numProcs > 1) PetscCall(PetscFree3(interfaces_lv, interfaces_gv, interfacesOffset));

    379:   /* Retrieve mesh from Mmg */
    380:   numCornersNew = 4;
    381:   PetscCallMMG_NONSTANDARD(PMMG_Get_meshSize, parmesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0);
    382:   PetscCall(PetscMalloc4(dim * numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer));
    383:   PetscCall(PetscMalloc3((dim + 1) * numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells));
    384:   PetscCall(PetscMalloc4(dim * numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces));
    385:   PetscCallMMG_NONSTANDARD(PMMG_Get_vertices, parmesh, verticesNew, verTagsNew, corners, requiredVer);
    386:   PetscCallMMG_NONSTANDARD(PMMG_Get_tetrahedra, parmesh, cellsNew, cellTagsNew, requiredCells);
    387:   PetscCallMMG_NONSTANDARD(PMMG_Get_triangles, parmesh, facesNew, faceTagsNew, requiredFaces);
    388:   PetscCall(PetscMalloc2(numVerticesNew, &owners, numVerticesNew, &gv_new));
    389:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_globalNum, 1);
    390:   PetscCallMMG_NONSTANDARD(PMMG_Get_verticesGloNum, parmesh, gv_new, owners);
    391:   for (i = 0; i < dim * numFacesNew; ++i) facesNew[i] -= 1;
    392:   for (i = 0; i < (dim + 1) * numCellsNew; ++i) cellsNew[i] = gv_new[cellsNew[i] - 1] - 1;
    393:   for (i = 0, numVerticesNewLoc = 0; i < numVerticesNew; ++i) {
    394:     if (owners[i] == rank) numVerticesNewLoc++;
    395:   }
    396:   PetscCall(PetscMalloc2(numVerticesNewLoc * dim, &verticesNewLoc, numVerticesNew, &verticesNewSorted));
    397:   for (i = 0, c = 0; i < numVerticesNew; i++) {
    398:     if (owners[i] == rank) {
    399:       for (j = 0; j < dim; ++j) verticesNewLoc[dim * c + j] = verticesNew[dim * i + j];
    400:       c++;
    401:     }
    402:   }

    404:   /* Reorder for consistency with DMPlex */
    405:   for (i = 0; i < numCellsNew; ++i) PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4 * i]));

    407:   /* Create new plex */
    408:   PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNewLoc, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNewLoc, NULL, &verticesNewSorted, dmNew));
    409:   PetscCallMMG_NONSTANDARD(PMMG_Free_all, PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_end);
    410:   PetscCall(PetscFree4(verticesNew, verTagsNew, corners, requiredVer));

    412:   /* Get adapted mesh information */
    413:   PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd));
    414:   PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd));
    415:   PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd));

    417:   /* Rebuild boundary label */
    418:   PetscCall(DMCreateLabel(*dmNew, flg ? bdName : bdLabelName));
    419:   PetscCall(DMGetLabel(*dmNew, flg ? bdName : bdLabelName, &bdLabelNew));
    420:   for (i = 0; i < numFacesNew; i++) {
    421:     PetscBool       hasTag = PETSC_FALSE;
    422:     PetscInt        numCoveredPoints, numFaces = 0, facePoints[3];
    423:     const PetscInt *coveredPoints = NULL;

    425:     for (j = 0; j < dim; ++j) {
    426:       lv = facesNew[i * dim + j];
    427:       gv = gv_new[lv] - 1;
    428:       PetscCall(PetscFindInt(gv, numVerticesNew, verticesNewSorted, &lv));
    429:       facePoints[j] = lv + vStart;
    430:     }
    431:     PetscCall(DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints));
    432:     for (j = 0; j < numCoveredPoints; ++j) {
    433:       if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) {
    434:         numFaces++;
    435:         f = j;
    436:       }
    437:     }
    438:     PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " vertices cannot define more than 1 facet (%" PetscInt_FMT ")", dim, numFaces);
    439:     PetscCall(DMLabelHasStratum(bdLabel, faceTagsNew[i], &hasTag));
    440:     if (hasTag) PetscCall(DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i]));
    441:     PetscCall(DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints));
    442:   }
    443:   PetscCall(PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces));
    444:   PetscCall(PetscFree2(owners, gv_new));
    445:   PetscCall(PetscFree2(verticesNewLoc, verticesNewSorted));
    446:   if (flg) PetscCall(DMLabelDestroy(&bdLabel));

    448:   /* Rebuild cell labels */
    449:   PetscCall(DMCreateLabel(*dmNew, rgLabel ? rgLabelName : rgName));
    450:   PetscCall(DMGetLabel(*dmNew, rgLabel ? rgLabelName : rgName, &rgLabelNew));
    451:   for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(rgLabelNew, c, cellTagsNew[c - cStart]));
    452:   PetscCall(PetscFree3(cellsNew, cellTagsNew, requiredCells));
    453:   PetscFunctionReturn(PETSC_SUCCESS);
    454: }