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