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

 42:   PetscFunctionBegin;
 43:   PetscCall(PetscCitationsRegister(ParMmgCitation, &ParMmgCite));
 44:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
 45:   PetscCallMPI(MPI_Comm_size(comm, &numProcs));
 46:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 47:   if (bdLabel) {
 48:     PetscCall(PetscObjectGetName((PetscObject)bdLabel, &bdLabelName));
 49:     PetscCall(PetscStrcmp(bdLabelName, bdName, &flg));
 50:     PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
 51:   }
 52:   if (rgLabel) {
 53:     PetscCall(PetscObjectGetName((PetscObject)rgLabel, &rgLabelName));
 54:     PetscCall(PetscStrcmp(rgLabelName, rgName, &flg));
 55:     PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for element tags", rgLabelName);
 56:   }

 58:   /* Get mesh information */
 59:   PetscCall(DMGetDimension(dm, &dim));
 60:   PetscCheck(dim == 3, comm, PETSC_ERR_ARG_OUTOFRANGE, "ParMmg only works in 3D.");
 61:   Neq = (dim * (dim + 1)) / 2;
 62:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 63:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
 64:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
 65:   PetscCall(DMPlexUninterpolate(dm, &udm));
 66:   PetscCall(DMPlexGetMaxSizes(udm, &maxConeSize, NULL));
 67:   numCells    = cEnd - cStart;
 68:   numVertices = vEnd - vStart;

 70:   /* Get cell offsets */
 71:   PetscCall(PetscMalloc1(numCells * maxConeSize, &cells));
 72:   for (c = 0, coff = 0; c < numCells; ++c) {
 73:     const PetscInt *cone;
 74:     PetscInt        coneSize, cl;

 76:     PetscCall(DMPlexGetConeSize(udm, c, &coneSize));
 77:     PetscCall(DMPlexGetCone(udm, c, &cone));
 78:     for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart + 1;
 79:   }

 81:   /* Get vertex coordinate array */
 82:   PetscCall(DMGetCoordinateDM(dm, &cdm));
 83:   PetscCall(DMGetLocalSection(cdm, &coordSection));
 84:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
 85:   PetscCall(VecGetArrayRead(coordinates, &coords));
 86:   PetscCall(PetscMalloc2(numVertices * Neq, &metric, dim * numVertices, &vertices));
 87:   for (v = 0; v < vEnd - vStart; ++v) {
 88:     PetscCall(PetscSectionGetOffset(coordSection, v + vStart, &off));
 89:     for (i = 0; i < dim; ++i) vertices[dim * v + i] = PetscRealPart(coords[off + i]);
 90:   }
 91:   PetscCall(VecRestoreArrayRead(coordinates, &coords));

 93:   /* Get face tags */
 94:   if (!bdLabel) {
 95:     flg = PETSC_TRUE;
 96:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabel));
 97:     PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabel));
 98:   }
 99:   PetscCall(DMLabelGetBounds(bdLabel, &pStart, &pEnd));
100:   for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
101:     PetscBool hasPoint;
102:     PetscInt *closure = NULL, closureSize, cl;

104:     PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint));
105:     if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;
106:     numFaceTags++;

108:     PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
109:     for (cl = 0; cl < closureSize * 2; cl += 2) {
110:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
111:     }
112:     PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
113:   }
114:   PetscCall(PetscMalloc2(bdSize, &bdFaces, numFaceTags, &faceTags));
115:   for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
116:     PetscBool hasPoint;
117:     PetscInt *closure = NULL, closureSize, cl;

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

122:     PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
123:     for (cl = 0; cl < closureSize * 2; cl += 2) {
124:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart + 1;
125:     }
126:     PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
127:     PetscCall(DMLabelGetValue(bdLabel, f, &faceTags[numFaceTags++]));
128:   }

130:   /* Get cell tags */
131:   PetscCall(PetscCalloc2(numVertices, &verTags, numCells, &cellTags));
132:   if (rgLabel) {
133:     for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelGetValue(rgLabel, c, &cellTags[c]));
134:   }

136:   /* Get metric */
137:   PetscCall(VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view"));
138:   PetscCall(VecGetArrayRead(vertexMetric, &met));
139:   PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
140:   PetscCall(DMPlexMetricIsUniform(dm, &uniform));
141:   for (v = 0; v < (vEnd - vStart); ++v) {
142:     for (i = 0, k = 0; i < dim; ++i) {
143:       for (j = i; j < dim; ++j, ++k) {
144:         if (isotropic) {
145:           if (i == j) {
146:             if (uniform) metric[Neq * v + k] = PetscRealPart(met[0]);
147:             else metric[Neq * v + k] = PetscRealPart(met[v]);
148:           } else metric[Neq * v + k] = 0.0;
149:         } else metric[Neq * v + k] = PetscRealPart(met[dim * dim * v + dim * i + j]);
150:       }
151:     }
152:   }
153:   PetscCall(VecRestoreArrayRead(vertexMetric, &met));

155:   /* Build ParMMG communicators: the list of vertices between two partitions  */
156:   niranks = nrranks = 0;
157:   numNgbRanks       = 0;
158:   if (numProcs > 1) {
159:     PetscCall(DMGetPointSF(dm, &sf));
160:     PetscCall(PetscSFSetUp(sf));
161:     PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, &ioffset, &irootloc));
162:     PetscCall(PetscSFGetRootRanks(sf, &nrranks, &rranks, &roffset, &rmine, &rremote));
163:     PetscCall(PetscCalloc1(numProcs, &numVerInterfaces));

165:     /* Count number of roots associated with each leaf */
166:     for (r = 0; r < niranks; ++r) {
167:       for (i = ioffset[r], count = 0; i < ioffset[r + 1]; ++i) {
168:         if (irootloc[i] >= vStart && irootloc[i] < vEnd) count++;
169:       }
170:       numVerInterfaces[iranks[r]] += count;
171:     }

173:     /* Count number of leaves associated with each root */
174:     for (r = 0; r < nrranks; ++r) {
175:       for (i = roffset[r], count = 0; i < roffset[r + 1]; ++i) {
176:         if (rmine[i] >= vStart && rmine[i] < vEnd) count++;
177:       }
178:       numVerInterfaces[rranks[r]] += count;
179:     }

181:     /* Count global number of ranks */
182:     for (p = 0; p < numProcs; ++p) {
183:       if (numVerInterfaces[p]) numNgbRanks++;
184:     }

186:     /* Provide numbers of vertex interfaces */
187:     PetscCall(PetscMalloc2(numNgbRanks, &ngbRanks, numNgbRanks, &verNgbRank));
188:     for (p = 0, n = 0; p < numProcs; ++p) {
189:       if (numVerInterfaces[p]) {
190:         ngbRanks[n]   = p;
191:         verNgbRank[n] = numVerInterfaces[p];
192:         n++;
193:       }
194:     }
195:     numVerNgbRanksTotal = 0;
196:     for (p = 0; p < numNgbRanks; ++p) numVerNgbRanksTotal += verNgbRank[p];

198:     /* For each neighbor, fill in interface arrays */
199:     PetscCall(PetscMalloc3(numVerNgbRanksTotal, &interfaces_lv, numVerNgbRanksTotal, &interfaces_gv, numNgbRanks + 1, &intOffset));
200:     intOffset[0] = 0;
201:     for (p = 0, r = 0, i = 0; p < numNgbRanks; ++p) {
202:       intOffset[p + 1] = intOffset[p];

204:       /* Leaf case */
205:       if (iranks && iranks[i] == ngbRanks[p]) {
206:         /* Add the right slice of irootloc at the right place */
207:         sliceSize = ioffset[i + 1] - ioffset[i];
208:         for (j = 0, count = 0; j < sliceSize; ++j) {
209:           PetscCheck(ioffset[i] + j < ioffset[niranks], comm, PETSC_ERR_ARG_OUTOFRANGE, "Leaf index %" PetscInt_FMT " out of range (expected < %" PetscInt_FMT ")", ioffset[i] + j, ioffset[niranks]);
210:           v = irootloc[ioffset[i] + j];
211:           if (v >= vStart && v < vEnd) {
212:             PetscCheck(intOffset[p + 1] + count < numVerNgbRanksTotal, comm, PETSC_ERR_ARG_OUTOFRANGE, "Leaf interface index %" PetscInt_FMT " out of range (expected < %" PetscInt_FMT ")", intOffset[p + 1] + count, numVerNgbRanksTotal);
213:             interfaces_lv[intOffset[p + 1] + count] = v - vStart;
214:             count++;
215:           }
216:         }
217:         intOffset[p + 1] += count;
218:         i++;
219:       }

221:       /* Root case */
222:       if (rranks && rranks[r] == ngbRanks[p]) {
223:         /* Add the right slice of rmine at the right place */
224:         sliceSize = roffset[r + 1] - roffset[r];
225:         for (j = 0, count = 0; j < sliceSize; ++j) {
226:           PetscCheck(roffset[r] + j < roffset[nrranks], comm, PETSC_ERR_ARG_OUTOFRANGE, "Root index %" PetscInt_FMT " out of range (expected < %" PetscInt_FMT ")", roffset[r] + j, roffset[nrranks]);
227:           v = rmine[roffset[r] + j];
228:           if (v >= vStart && v < vEnd) {
229:             PetscCheck(intOffset[p + 1] + count < numVerNgbRanksTotal, comm, PETSC_ERR_ARG_OUTOFRANGE, "Root interface index %" PetscInt_FMT " out of range (expected < %" PetscInt_FMT ")", intOffset[p + 1] + count, numVerNgbRanksTotal);
230:             interfaces_lv[intOffset[p + 1] + count] = v - vStart;
231:             count++;
232:           }
233:         }
234:         intOffset[p + 1] += count;
235:         r++;
236:       }

238:       /* Check validity of offsets */
239:       PetscCheck(intOffset[p + 1] == intOffset[p] + verNgbRank[p], comm, PETSC_ERR_ARG_OUTOFRANGE, "Missing offsets (expected %" PetscInt_FMT ", got %" PetscInt_FMT ")", intOffset[p] + verNgbRank[p], intOffset[p + 1]);
240:     }
241:     PetscCall(DMPlexGetVertexNumbering(udm, &globalVertexNum));
242:     PetscCall(ISGetIndices(globalVertexNum, &gV));
243:     for (i = 0; i < numVerNgbRanksTotal; ++i) {
244:       v                = gV[interfaces_lv[i]];
245:       interfaces_gv[i] = v < 0 ? -v - 1 : v;
246:       interfaces_lv[i] += 1;
247:       interfaces_gv[i] += 1;
248:     }
249:     PetscCall(ISRestoreIndices(globalVertexNum, &gV));
250:     PetscCall(PetscFree(numVerInterfaces));
251:   }
252:   PetscCall(DMDestroy(&udm));

254:   /* Send the data to ParMmg and remesh */
255:   PetscCall(DMPlexMetricNoInsertion(dm, &noInsert));
256:   PetscCall(DMPlexMetricNoSwapping(dm, &noSwap));
257:   PetscCall(DMPlexMetricNoMovement(dm, &noMove));
258:   PetscCall(DMPlexMetricNoSurf(dm, &noSurf));
259:   PetscCall(DMPlexMetricGetVerbosity(dm, &verbosity));
260:   PetscCall(DMPlexMetricGetNumIterations(dm, &numIter));
261:   PetscCall(DMPlexMetricGetGradationFactor(dm, &gradationFactor));
262:   PetscCall(DMPlexMetricGetHausdorffNumber(dm, &hausdorffNumber));
263:   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);
264:   PetscCallMMG_NONSTANDARD(PMMG_Set_meshSize, parmesh, numVertices, numCells, 0, numFaceTags, 0, 0);
265:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_APImode, PMMG_APIDISTRIB_nodes);
266:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_noinsert, noInsert);
267:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_noswap, noSwap);
268:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_nomove, noMove);
269:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_nosurf, noSurf);
270:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_verbose, verbosity);
271:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_globalNum, 1);
272:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_niter, numIter);
273:   PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter, parmesh, PMMG_DPARAM_hgrad, gradationFactor);
274:   PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter, parmesh, PMMG_DPARAM_hausd, hausdorffNumber);
275:   PetscCallMMG_NONSTANDARD(PMMG_Set_vertices, parmesh, vertices, verTags);
276:   PetscCallMMG_NONSTANDARD(PMMG_Set_tetrahedra, parmesh, cells, cellTags);
277:   PetscCallMMG_NONSTANDARD(PMMG_Set_triangles, parmesh, bdFaces, faceTags);
278:   PetscCallMMG_NONSTANDARD(PMMG_Set_metSize, parmesh, MMG5_Vertex, numVertices, MMG5_Tensor);
279:   PetscCallMMG_NONSTANDARD(PMMG_Set_tensorMets, parmesh, metric);
280:   PetscCallMMG_NONSTANDARD(PMMG_Set_numberOfNodeCommunicators, parmesh, numNgbRanks);
281:   for (c = 0; c < numNgbRanks; ++c) {
282:     PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicatorSize, parmesh, c, ngbRanks[c], intOffset[c + 1] - intOffset[c]);
283:     PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicator_nodes, parmesh, c, &interfaces_lv[intOffset[c]], &interfaces_gv[intOffset[c]], 1);
284:   }
285:   PetscCallMMG(PMMG_parmmglib_distributed, parmesh);
286:   PetscCall(PetscFree(cells));
287:   PetscCall(PetscFree2(metric, vertices));
288:   PetscCall(PetscFree2(bdFaces, faceTags));
289:   PetscCall(PetscFree2(verTags, cellTags));
290:   if (numProcs > 1) {
291:     PetscCall(PetscFree2(ngbRanks, verNgbRank));
292:     PetscCall(PetscFree3(interfaces_lv, interfaces_gv, intOffset));
293:   }

295:   /* Retrieve mesh from Mmg */
296:   numCornersNew = 4;
297:   PetscCallMMG_NONSTANDARD(PMMG_Get_meshSize, parmesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0);
298:   PetscCall(PetscMalloc4(dim * numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer));
299:   PetscCall(PetscMalloc3((dim + 1) * numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells));
300:   PetscCall(PetscMalloc4(dim * numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces));
301:   PetscCallMMG_NONSTANDARD(PMMG_Get_vertices, parmesh, verticesNew, verTagsNew, corners, requiredVer);
302:   PetscCallMMG_NONSTANDARD(PMMG_Get_tetrahedra, parmesh, cellsNew, cellTagsNew, requiredCells);
303:   PetscCallMMG_NONSTANDARD(PMMG_Get_triangles, parmesh, facesNew, faceTagsNew, requiredFaces);
304:   PetscCall(PetscMalloc2(numVerticesNew, &owners, numVerticesNew, &gv_new));
305:   PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_globalNum, 1);
306:   PetscCallMMG_NONSTANDARD(PMMG_Get_verticesGloNum, parmesh, gv_new, owners);
307:   for (i = 0; i < dim * numFacesNew; ++i) facesNew[i] -= 1;
308:   for (i = 0; i < (dim + 1) * numCellsNew; ++i) cellsNew[i] = gv_new[cellsNew[i] - 1] - 1;
309:   for (i = 0, numVerticesNewLoc = 0; i < numVerticesNew; ++i) {
310:     if (owners[i] == rank) numVerticesNewLoc++;
311:   }
312:   PetscCall(PetscMalloc2(numVerticesNewLoc * dim, &verticesNewLoc, numVerticesNew, &verticesNewSorted));
313:   for (i = 0, c = 0; i < numVerticesNew; i++) {
314:     if (owners[i] == rank) {
315:       for (j = 0; j < dim; ++j) verticesNewLoc[dim * c + j] = verticesNew[dim * i + j];
316:       c++;
317:     }
318:   }

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

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

328:   /* Get adapted mesh information */
329:   PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd));
330:   PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd));
331:   PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd));

333:   /* Rebuild boundary label */
334:   PetscCall(DMCreateLabel(*dmNew, flg ? bdName : bdLabelName));
335:   PetscCall(DMGetLabel(*dmNew, flg ? bdName : bdLabelName, &bdLabelNew));
336:   for (i = 0; i < numFacesNew; i++) {
337:     PetscBool       hasTag = PETSC_FALSE;
338:     PetscInt        numCoveredPoints, numFaces = 0, facePoints[3];
339:     const PetscInt *coveredPoints = NULL;

341:     for (j = 0; j < dim; ++j) {
342:       lv = facesNew[i * dim + j];
343:       gv = gv_new[lv] - 1;
344:       PetscCall(PetscFindInt(gv, numVerticesNew, verticesNewSorted, &lv));
345:       facePoints[j] = lv + vStart;
346:     }
347:     PetscCall(DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints));
348:     for (j = 0; j < numCoveredPoints; ++j) {
349:       if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) {
350:         numFaces++;
351:         f = j;
352:       }
353:     }
354:     PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " vertices cannot define more than 1 facet (%" PetscInt_FMT ")", dim, numFaces);
355:     PetscCall(DMLabelHasStratum(bdLabel, faceTagsNew[i], &hasTag));
356:     if (hasTag) PetscCall(DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i]));
357:     PetscCall(DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints));
358:   }
359:   PetscCall(PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces));
360:   PetscCall(PetscFree2(owners, gv_new));
361:   PetscCall(PetscFree2(verticesNewLoc, verticesNewSorted));
362:   if (flg) PetscCall(DMLabelDestroy(&bdLabel));

364:   /* Rebuild cell labels */
365:   PetscCall(DMCreateLabel(*dmNew, rgLabel ? rgLabelName : rgName));
366:   PetscCall(DMGetLabel(*dmNew, rgLabel ? rgLabelName : rgName, &rgLabelNew));
367:   for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(rgLabelNew, c, cellTagsNew[c - cStart]));
368:   PetscCall(PetscFree3(cellsNew, cellTagsNew, requiredCells));

370:   PetscFunctionReturn(PETSC_SUCCESS);
371: }