1 #include "../mmgcommon.h" /*I "petscdmplex.h" I*/ 2 #include <parmmg/libparmmg.h> 3 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"; 12 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; 41 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 } 57 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; 69 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; 75 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 } 80 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)); 92 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; 103 104 PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint)); 105 if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue; 106 numFaceTags++; 107 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; 118 119 PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint)); 120 if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue; 121 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 } 129 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 } 135 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)); 154 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)); 164 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 } 172 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 } 180 181 /* Count global number of ranks */ 182 for (p = 0; p < numProcs; ++p) { 183 if (numVerInterfaces[p]) numNgbRanks++; 184 } 185 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]; 197 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]; 203 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 } 220 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 } 237 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)); 253 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 } 294 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 } 319 320 /* Reorder for consistency with DMPlex */ 321 for (i = 0; i < numCellsNew; ++i) PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4 * i])); 322 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)); 327 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)); 332 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; 340 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)); 363 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)); 369 370 PetscFunctionReturn(PETSC_SUCCESS); 371 } 372