1 #include "../mmgcommon.h" /*I "petscdmplex.h" I*/ 2 #include <parmmg/libparmmg.h> 3 4 PETSC_EXTERN PetscErrorCode DMAdaptMetric_ParMmg_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DMLabel rgLabel, DM *dmNew) 5 { 6 MPI_Comm comm; 7 const char *bdName = "_boundary_"; 8 const char *rgName = "_regions_"; 9 DM udm, cdm; 10 DMLabel bdLabelNew, rgLabelNew; 11 const char *bdLabelName, *rgLabelName; 12 IS globalVertexNum; 13 PetscSection coordSection; 14 Vec coordinates; 15 PetscSF sf; 16 const PetscScalar *coords, *met; 17 PetscReal *vertices, *metric, *verticesNew, *verticesNewLoc, gradationFactor, hausdorffNumber; 18 PetscInt *cells, *cellsNew, *cellTags, *cellTagsNew, *verTags, *verTagsNew; 19 PetscInt *bdFaces, *faceTags, *facesNew, *faceTagsNew; 20 PetscInt *corners, *requiredCells, *requiredVer, *ridges, *requiredFaces; 21 PetscInt cStart, cEnd, c, numCells, fStart, fEnd, f, numFaceTags, vStart, vEnd, v, numVertices; 22 PetscInt dim, off, coff, maxConeSize, bdSize, i, j, k, Neq, verbosity, numIter; 23 PetscInt *numVerInterfaces, *ngbRanks, *verNgbRank, *interfaces_lv, *interfaces_gv, *intOffset; 24 PetscInt niranks, nrranks, numNgbRanks, numVerNgbRanksTotal, count, sliceSize, p, r, n, lv, gv; 25 PetscInt *gv_new, *owners, *verticesNewSorted, pStart, pEnd; 26 PetscInt numCellsNew, numVerticesNew, numCornersNew, numFacesNew, numVerticesNewLoc; 27 const PetscInt *gV, *ioffset, *irootloc, *roffset, *rmine, *rremote; 28 PetscBool flg = PETSC_FALSE, noInsert, noSwap, noMove, noSurf, isotropic, uniform; 29 const PetscMPIInt *iranks, *rranks; 30 PetscMPIInt numProcs, rank; 31 PMMG_pParMesh parmesh = NULL; 32 33 PetscFunctionBegin; 34 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 35 PetscCallMPI(MPI_Comm_size(comm, &numProcs)); 36 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 37 if (bdLabel) { 38 PetscCall(PetscObjectGetName((PetscObject) bdLabel, &bdLabelName)); 39 PetscCall(PetscStrcmp(bdLabelName, bdName, &flg)); 40 PetscCheck(!flg,comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 41 } 42 if (rgLabel) { 43 PetscCall(PetscObjectGetName((PetscObject) rgLabel, &rgLabelName)); 44 PetscCall(PetscStrcmp(rgLabelName, rgName, &flg)); 45 PetscCheck(!flg,comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for element tags", rgLabelName); 46 } 47 48 /* Get mesh information */ 49 PetscCall(DMGetDimension(dm, &dim)); 50 PetscCheck(dim == 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "ParMmg only works in 3D."); 51 Neq = (dim*(dim+1))/2; 52 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 53 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 54 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 55 PetscCall(DMPlexUninterpolate(dm, &udm)); 56 PetscCall(DMPlexGetMaxSizes(udm, &maxConeSize, NULL)); 57 numCells = cEnd - cStart; 58 numVertices = vEnd - vStart; 59 60 /* Get cell offsets */ 61 PetscCall(PetscMalloc1(numCells*maxConeSize, &cells)); 62 for (c = 0, coff = 0; c < numCells; ++c) { 63 const PetscInt *cone; 64 PetscInt coneSize, cl; 65 66 PetscCall(DMPlexGetConeSize(udm, c, &coneSize)); 67 PetscCall(DMPlexGetCone(udm, c, &cone)); 68 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart+1; 69 } 70 71 /* Get vertex coordinate array */ 72 PetscCall(DMGetCoordinateDM(dm, &cdm)); 73 PetscCall(DMGetLocalSection(cdm, &coordSection)); 74 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 75 PetscCall(VecGetArrayRead(coordinates, &coords)); 76 PetscCall(PetscMalloc2(numVertices*Neq, &metric, dim*numVertices, &vertices)); 77 for (v = 0; v < vEnd-vStart; ++v) { 78 PetscCall(PetscSectionGetOffset(coordSection, v+vStart, &off)); 79 for (i = 0; i < dim; ++i) vertices[dim*v+i] = PetscRealPart(coords[off+i]); 80 } 81 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 82 83 /* Get face tags */ 84 if (!bdLabel) { 85 flg = PETSC_TRUE; 86 PetscCall(DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabel)); 87 PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabel)); 88 } 89 PetscCall(DMLabelGetBounds(bdLabel, &pStart, &pEnd)); 90 for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) { 91 PetscBool hasPoint; 92 PetscInt *closure = NULL, closureSize, cl; 93 94 PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint)); 95 if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue; 96 numFaceTags++; 97 98 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 99 for (cl = 0; cl < closureSize*2; cl += 2) { 100 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 101 } 102 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 103 } 104 PetscCall(PetscMalloc2(bdSize, &bdFaces, numFaceTags, &faceTags)); 105 for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) { 106 PetscBool hasPoint; 107 PetscInt *closure = NULL, closureSize, cl; 108 109 PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint)); 110 if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue; 111 112 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 113 for (cl = 0; cl < closureSize*2; cl += 2) { 114 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart+1; 115 } 116 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 117 PetscCall(DMLabelGetValue(bdLabel, f, &faceTags[numFaceTags++])); 118 } 119 120 /* Get cell tags */ 121 PetscCall(PetscCalloc2(numVertices, &verTags, numCells, &cellTags)); 122 if (rgLabel) { 123 for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelGetValue(rgLabel, c, &cellTags[c])); 124 } 125 126 /* Get metric */ 127 PetscCall(VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view")); 128 PetscCall(VecGetArrayRead(vertexMetric, &met)); 129 PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic)); 130 PetscCall(DMPlexMetricIsUniform(dm, &uniform)); 131 for (v = 0; v < (vEnd-vStart); ++v) { 132 for (i = 0, k = 0; i < dim; ++i) { 133 for (j = i; j < dim; ++j, ++k) { 134 if (isotropic) { 135 if (i == j) { 136 if (uniform) metric[Neq*v+k] = PetscRealPart(met[0]); 137 else metric[Neq*v+k] = PetscRealPart(met[v]); 138 } else metric[Neq*v+k] = 0.0; 139 } else metric[Neq*v+k] = PetscRealPart(met[dim*dim*v+dim*i+j]); 140 } 141 } 142 } 143 PetscCall(VecRestoreArrayRead(vertexMetric, &met)); 144 145 /* Build ParMMG communicators: the list of vertices between two partitions */ 146 niranks = nrranks = 0; 147 numNgbRanks = 0; 148 if (numProcs > 1) { 149 PetscCall(DMGetPointSF(dm, &sf)); 150 PetscCall(PetscSFSetUp(sf)); 151 PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, &ioffset, &irootloc)); 152 PetscCall(PetscSFGetRootRanks(sf, &nrranks, &rranks, &roffset, &rmine, &rremote)); 153 PetscCall(PetscCalloc1(numProcs, &numVerInterfaces)); 154 155 /* Count number of roots associated with each leaf */ 156 for (r = 0; r < niranks; ++r) { 157 for (i=ioffset[r], count=0; i<ioffset[r+1]; ++i) { 158 if (irootloc[i] >= vStart && irootloc[i] < vEnd) count++; 159 } 160 numVerInterfaces[iranks[r]] += count; 161 } 162 163 /* Count number of leaves associated with each root */ 164 for (r = 0; r < nrranks; ++r) { 165 for (i=roffset[r], count=0; i<roffset[r+1]; ++i) { 166 if (rmine[i] >= vStart && rmine[i] < vEnd) count++; 167 } 168 numVerInterfaces[rranks[r]] += count; 169 } 170 171 /* Count global number of ranks */ 172 for (p = 0; p < numProcs; ++p) { 173 if (numVerInterfaces[p]) numNgbRanks++; 174 } 175 176 /* Provide numbers of vertex interfaces */ 177 PetscCall(PetscMalloc2(numNgbRanks, &ngbRanks, numNgbRanks, &verNgbRank)); 178 for (p = 0, n = 0; p < numProcs; ++p) { 179 if (numVerInterfaces[p]) { 180 ngbRanks[n] = p; 181 verNgbRank[n] = numVerInterfaces[p]; 182 n++; 183 } 184 } 185 numVerNgbRanksTotal = 0; 186 for (p = 0; p < numNgbRanks; ++p) numVerNgbRanksTotal += verNgbRank[p]; 187 188 /* For each neighbor, fill in interface arrays */ 189 PetscCall(PetscMalloc3(numVerNgbRanksTotal, &interfaces_lv, numVerNgbRanksTotal, &interfaces_gv, numNgbRanks+1, &intOffset)); 190 intOffset[0] = 0; 191 for (p = 0, r = 0, i = 0; p < numNgbRanks; ++p) { 192 intOffset[p+1] = intOffset[p]; 193 194 /* Leaf case */ 195 if (iranks && iranks[i] == ngbRanks[p]) { 196 197 /* Add the right slice of irootloc at the right place */ 198 sliceSize = ioffset[i+1]-ioffset[i]; 199 for (j = 0, count = 0; j < sliceSize; ++j) { 200 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]); 201 v = irootloc[ioffset[i]+j]; 202 if (v >= vStart && v < vEnd) { 203 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); 204 interfaces_lv[intOffset[p+1]+count] = v-vStart; 205 count++; 206 } 207 } 208 intOffset[p+1] += count; 209 i++; 210 } 211 212 /* Root case */ 213 if (rranks && rranks[r] == ngbRanks[p]) { 214 215 /* Add the right slice of rmine at the right place */ 216 sliceSize = roffset[r+1]-roffset[r]; 217 for (j = 0, count = 0; j < sliceSize; ++j) { 218 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]); 219 v = rmine[roffset[r]+j]; 220 if (v >= vStart && v < vEnd) { 221 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); 222 interfaces_lv[intOffset[p+1]+count] = v-vStart; 223 count++; 224 } 225 } 226 intOffset[p+1] += count; 227 r++; 228 } 229 230 /* Check validity of offsets */ 231 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]); 232 } 233 PetscCall(DMPlexGetVertexNumbering(udm, &globalVertexNum)); 234 PetscCall(ISGetIndices(globalVertexNum, &gV)); 235 for (i = 0; i < numVerNgbRanksTotal; ++i) { 236 v = gV[interfaces_lv[i]]; 237 interfaces_gv[i] = v < 0 ? -v-1 : v; 238 interfaces_lv[i] += 1; 239 interfaces_gv[i] += 1; 240 } 241 PetscCall(ISRestoreIndices(globalVertexNum, &gV)); 242 PetscCall(PetscFree(numVerInterfaces)); 243 } 244 PetscCall(DMDestroy(&udm)); 245 246 /* Send the data to ParMmg and remesh */ 247 PetscCall(DMPlexMetricNoInsertion(dm, &noInsert)); 248 PetscCall(DMPlexMetricNoSwapping(dm, &noSwap)); 249 PetscCall(DMPlexMetricNoMovement(dm, &noMove)); 250 PetscCall(DMPlexMetricNoSurf(dm, &noSurf)); 251 PetscCall(DMPlexMetricGetVerbosity(dm, &verbosity)); 252 PetscCall(DMPlexMetricGetNumIterations(dm, &numIter)); 253 PetscCall(DMPlexMetricGetGradationFactor(dm, &gradationFactor)); 254 PetscCall(DMPlexMetricGetHausdorffNumber(dm, &hausdorffNumber)); 255 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)); 256 PetscCallMMG_NONSTANDARD(PMMG_Set_meshSize(parmesh, numVertices, numCells, 0, numFaceTags, 0, 0)); 257 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter(parmesh, PMMG_IPARAM_APImode, PMMG_APIDISTRIB_nodes)); 258 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter(parmesh, PMMG_IPARAM_noinsert, noInsert)); 259 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter(parmesh, PMMG_IPARAM_noswap, noSwap)); 260 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter(parmesh, PMMG_IPARAM_nomove, noMove)); 261 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter(parmesh, PMMG_IPARAM_nosurf, noSurf)); 262 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter(parmesh, PMMG_IPARAM_verbose, verbosity)); 263 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter(parmesh, PMMG_IPARAM_globalNum, 1)); 264 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter(parmesh, PMMG_IPARAM_niter, numIter)); 265 PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter(parmesh, PMMG_DPARAM_hgrad, gradationFactor)); 266 PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter(parmesh, PMMG_DPARAM_hausd, hausdorffNumber)); 267 PetscCallMMG_NONSTANDARD(PMMG_Set_vertices(parmesh, vertices, verTags)); 268 PetscCallMMG_NONSTANDARD(PMMG_Set_tetrahedra(parmesh, cells, cellTags)); 269 PetscCallMMG_NONSTANDARD(PMMG_Set_triangles(parmesh, bdFaces, faceTags)); 270 PetscCallMMG_NONSTANDARD(PMMG_Set_metSize(parmesh, MMG5_Vertex, numVertices, MMG5_Tensor)); 271 PetscCallMMG_NONSTANDARD(PMMG_Set_tensorMets(parmesh, metric)); 272 PetscCallMMG_NONSTANDARD(PMMG_Set_numberOfNodeCommunicators(parmesh, numNgbRanks)); 273 for (c = 0; c < numNgbRanks; ++c) { 274 PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicatorSize(parmesh, c, ngbRanks[c], intOffset[c+1]-intOffset[c])); 275 PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicator_nodes(parmesh, c, &interfaces_lv[intOffset[c]], &interfaces_gv[intOffset[c]], 1)); 276 } 277 PetscCallMMG(PMMG_parmmglib_distributed(parmesh)); 278 PetscCall(PetscFree(cells)); 279 PetscCall(PetscFree2(metric, vertices)); 280 PetscCall(PetscFree2(bdFaces, faceTags)); 281 PetscCall(PetscFree2(verTags, cellTags)); 282 if (numProcs > 1) { 283 PetscCall(PetscFree2(ngbRanks, verNgbRank)); 284 PetscCall(PetscFree3(interfaces_lv, interfaces_gv, intOffset)); 285 } 286 287 /* Retrieve mesh from Mmg */ 288 numCornersNew = 4; 289 PetscCallMMG_NONSTANDARD(PMMG_Get_meshSize(parmesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0)); 290 PetscCall(PetscMalloc4(dim*numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer)); 291 PetscCall(PetscMalloc3((dim+1)*numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells)); 292 PetscCall(PetscMalloc4(dim*numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces)); 293 PetscCallMMG_NONSTANDARD(PMMG_Get_vertices(parmesh, verticesNew, verTagsNew, corners, requiredVer)); 294 PetscCallMMG_NONSTANDARD(PMMG_Get_tetrahedra(parmesh, cellsNew, cellTagsNew, requiredCells)); 295 PetscCallMMG_NONSTANDARD(PMMG_Get_triangles(parmesh, facesNew, faceTagsNew, requiredFaces)); 296 PetscCall(PetscMalloc2(numVerticesNew, &owners, numVerticesNew, &gv_new)); 297 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter(parmesh, PMMG_IPARAM_globalNum, 1)); 298 PetscCallMMG_NONSTANDARD(PMMG_Get_verticesGloNum(parmesh, gv_new, owners)); 299 for (i = 0; i < dim*numFacesNew; ++i) facesNew[i] -= 1; 300 for (i = 0; i < (dim+1)*numCellsNew; ++i) cellsNew[i] = gv_new[cellsNew[i]-1]-1; 301 for (i = 0, numVerticesNewLoc = 0; i < numVerticesNew; ++i) { 302 if (owners[i] == rank) numVerticesNewLoc++; 303 } 304 PetscCall(PetscMalloc2(numVerticesNewLoc*dim, &verticesNewLoc, numVerticesNew, &verticesNewSorted)); 305 for (i = 0, c = 0; i < numVerticesNew; i++) { 306 if (owners[i] == rank) { 307 for (j=0; j<dim; ++j) verticesNewLoc[dim*c+j] = verticesNew[dim*i+j]; 308 c++; 309 } 310 } 311 312 /* Reorder for consistency with DMPlex */ 313 for (i = 0; i < numCellsNew; ++i) PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4*i])); 314 315 /* Create new plex */ 316 PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNewLoc, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNewLoc, NULL, &verticesNewSorted, dmNew)); 317 PetscCallMMG_NONSTANDARD(PMMG_Free_all(PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_end)); 318 PetscCall(PetscFree4(verticesNew, verTagsNew, corners, requiredVer)); 319 320 /* Get adapted mesh information */ 321 PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd)); 322 PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd)); 323 PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd)); 324 325 /* Rebuild boundary label */ 326 PetscCall(DMCreateLabel(*dmNew, flg ? bdName : bdLabelName)); 327 PetscCall(DMGetLabel(*dmNew, flg ? bdName : bdLabelName, &bdLabelNew)); 328 for (i = 0; i < numFacesNew; i++) { 329 PetscBool hasTag = PETSC_FALSE; 330 PetscInt numCoveredPoints, numFaces = 0, facePoints[3]; 331 const PetscInt *coveredPoints = NULL; 332 333 for (j = 0; j < dim; ++j) { 334 lv = facesNew[i*dim+j]; 335 gv = gv_new[lv]-1; 336 PetscCall(PetscFindInt(gv, numVerticesNew, verticesNewSorted, &lv)); 337 facePoints[j] = lv+vStart; 338 } 339 PetscCall(DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints)); 340 for (j = 0; j < numCoveredPoints; ++j) { 341 if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) { 342 numFaces++; 343 f = j; 344 } 345 } 346 PetscCheck(numFaces == 1,comm, PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " vertices cannot define more than 1 facet (%" PetscInt_FMT ")", dim, numFaces); 347 PetscCall(DMLabelHasStratum(bdLabel, faceTagsNew[i], &hasTag)); 348 if (hasTag) PetscCall(DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i])); 349 PetscCall(DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints)); 350 } 351 PetscCall(PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces)); 352 PetscCall(PetscFree2(owners, gv_new)); 353 PetscCall(PetscFree2(verticesNewLoc, verticesNewSorted)); 354 if (flg) PetscCall(DMLabelDestroy(&bdLabel)); 355 356 /* Rebuild cell labels */ 357 PetscCall(DMCreateLabel(*dmNew, rgLabel ? rgLabelName : rgName)); 358 PetscCall(DMGetLabel(*dmNew, rgLabel ? rgLabelName : rgName, &rgLabelNew)); 359 for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(rgLabelNew, c, cellTagsNew[c-cStart])); 360 PetscCall(PetscFree3(cellsNew, cellTagsNew, requiredCells)); 361 362 PetscFunctionReturn(0); 363 } 364