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 /* 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 PetscInt niranks, nrranks, numNgbRanks, r, lv, gv; 38 PetscInt *gv_new, *owners, *verticesNewSorted, pStart, pEnd; 39 PetscInt numCellsNew, numVerticesNew, numCornersNew, numFacesNew, numVerticesNewLoc; 40 const PetscInt *gV, *ioffset, *irootloc, *roffset, *rmine, *rremote; 41 PetscBool flg = PETSC_FALSE, noInsert, noSwap, noMove, noSurf, isotropic, uniform; 42 const PetscMPIInt *iranks, *rranks; 43 PetscMPIInt numProcs, rank; 44 PMMG_pParMesh parmesh = NULL; 45 46 // DEVELOPER NOTE: ParMmg wants to know the rank of every process which is sharing a given point and 47 // for this information to be conveyed to every process that is sharing that point. 48 49 PetscFunctionBegin; 50 PetscCall(PetscCitationsRegister(ParMmgCitation, &ParMmgCite)); 51 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 52 PetscCallMPI(MPI_Comm_size(comm, &numProcs)); 53 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 54 if (bdLabel) { 55 PetscCall(PetscObjectGetName((PetscObject)bdLabel, &bdLabelName)); 56 PetscCall(PetscStrcmp(bdLabelName, bdName, &flg)); 57 PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 58 } 59 if (rgLabel) { 60 PetscCall(PetscObjectGetName((PetscObject)rgLabel, &rgLabelName)); 61 PetscCall(PetscStrcmp(rgLabelName, rgName, &flg)); 62 PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for element tags", rgLabelName); 63 } 64 65 /* Get mesh information */ 66 PetscCall(DMGetDimension(dm, &dim)); 67 PetscCheck(dim == 3, comm, PETSC_ERR_ARG_OUTOFRANGE, "ParMmg only works in 3D."); 68 Neq = (dim * (dim + 1)) / 2; 69 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 70 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 71 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 72 PetscCall(DMPlexUninterpolate(dm, &udm)); 73 PetscCall(DMPlexGetMaxSizes(udm, &maxConeSize, NULL)); 74 numCells = cEnd - cStart; 75 numVertices = vEnd - vStart; 76 77 /* Get parallel data; work out which cells are owned and which are leaves */ 78 PetscCall(PetscCalloc1(numCells, &cIsLeaf)); 79 numCellsNotShared = numCells; 80 niranks = nrranks = 0; 81 if (numProcs > 1) { 82 PetscCall(DMGetPointSF(dm, &sf)); 83 PetscCall(PetscSFSetUp(sf)); 84 PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, &ioffset, &irootloc)); 85 PetscCall(PetscSFGetRootRanks(sf, &nrranks, &rranks, &roffset, &rmine, &rremote)); 86 for (r = 0; r < nrranks; ++r) { 87 for (i = roffset[r]; i < roffset[r + 1]; ++i) { 88 if (rmine[i] >= cStart && rmine[i] < cEnd) { 89 cIsLeaf[rmine[i] - cStart] = 1; 90 numCellsNotShared--; 91 } 92 } 93 } 94 } 95 96 /* 97 Create a vertex numbering for ParMmg starting at 1. Vertices not included in any 98 owned cell remain 0 and will be removed. Using this numbering, create cells. 99 */ 100 numUsedVertices = 0; 101 PetscCall(PetscCalloc1(numVertices, &vertexNumber)); 102 PetscCall(PetscMalloc1(numCellsNotShared * maxConeSize, &cells)); 103 for (c = 0, coff = 0; c < numCells; ++c) { 104 const PetscInt *cone; 105 PetscInt coneSize, cl; 106 107 if (!cIsLeaf[c]) { 108 PetscCall(DMPlexGetConeSize(udm, cStart + c, &coneSize)); 109 PetscCall(DMPlexGetCone(udm, cStart + c, &cone)); 110 for (cl = 0; cl < coneSize; ++cl) { 111 if (!vertexNumber[cone[cl] - vStart]) vertexNumber[cone[cl] - vStart] = ++numUsedVertices; 112 cells[coff++] = vertexNumber[cone[cl] - vStart]; 113 } 114 } 115 } 116 117 /* Get array of vertex coordinates */ 118 PetscCall(DMGetCoordinateDM(dm, &cdm)); 119 PetscCall(DMGetLocalSection(cdm, &coordSection)); 120 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 121 PetscCall(VecGetArrayRead(coordinates, &coords)); 122 PetscCall(PetscMalloc2(numUsedVertices * Neq, &metric, dim * numUsedVertices, &vertices)); 123 for (v = 0; v < vEnd - vStart; ++v) { 124 PetscCall(PetscSectionGetOffset(coordSection, v + vStart, &off)); 125 if (vertexNumber[v]) { 126 for (i = 0; i < dim; ++i) vertices[dim * (vertexNumber[v] - 1) + i] = PetscRealPart(coords[off + i]); 127 } 128 } 129 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 130 131 /* Get face tags */ 132 if (!bdLabel) { 133 flg = PETSC_TRUE; 134 PetscCall(DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabel)); 135 PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabel)); 136 } 137 PetscCall(DMLabelGetBounds(bdLabel, &pStart, &pEnd)); 138 PetscCall(PetscCalloc1(pEnd - pStart, &fIsIncluded)); 139 for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) { 140 PetscBool hasPoint; 141 PetscInt *closure = NULL, closureSize, cl; 142 143 PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint)); 144 if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue; 145 146 /* Only faces adjacent to an owned (non-leaf) cell are included */ 147 PetscInt nnbrs; 148 const PetscInt *nbrs; 149 PetscCall(DMPlexGetSupportSize(dm, f, &nnbrs)); 150 PetscCall(DMPlexGetSupport(dm, f, &nbrs)); 151 for (c = 0; c < nnbrs; ++c) fIsIncluded[f - pStart] = fIsIncluded[f - pStart] || !cIsLeaf[nbrs[c]]; 152 if (!fIsIncluded[f - pStart]) continue; 153 154 numFaceTags++; 155 156 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 157 for (cl = 0; cl < closureSize * 2; cl += 2) { 158 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 159 } 160 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 161 } 162 PetscCall(PetscMalloc2(bdSize, &bdFaces, numFaceTags, &faceTags)); 163 for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) { 164 PetscInt *closure = NULL, closureSize, cl; 165 166 if (!fIsIncluded[f - pStart]) continue; 167 168 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 169 for (cl = 0; cl < closureSize * 2; cl += 2) { 170 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = vertexNumber[closure[cl] - vStart]; 171 } 172 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 173 PetscCall(DMLabelGetValue(bdLabel, f, &faceTags[numFaceTags++])); 174 } 175 PetscCall(PetscFree(fIsIncluded)); 176 177 /* Get cell tags */ 178 PetscCall(PetscCalloc2(numUsedVertices, &verTags, numCellsNotShared, &cellTags)); 179 if (rgLabel) { 180 for (c = cStart, coff = 0; c < cEnd; ++c) { 181 if (!cIsLeaf[c - cStart]) { PetscCall(DMLabelGetValue(rgLabel, c, &cellTags[coff++])); } 182 } 183 } 184 PetscCall(PetscFree(cIsLeaf)); 185 186 /* Get metric, using only the upper triangular part */ 187 PetscCall(VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view")); 188 PetscCall(VecGetArrayRead(vertexMetric, &met)); 189 PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic)); 190 PetscCall(DMPlexMetricIsUniform(dm, &uniform)); 191 for (v = 0; v < (vEnd - vStart); ++v) { 192 PetscInt vv = vertexNumber[v]; 193 if (!vv--) continue; 194 for (i = 0, k = 0; i < dim; ++i) { 195 for (j = i; j < dim; ++j, ++k) { 196 if (isotropic) { 197 if (i == j) { 198 if (uniform) metric[Neq * vv + k] = PetscRealPart(met[0]); 199 else metric[Neq * vv + k] = PetscRealPart(met[v]); 200 } else metric[Neq * vv + k] = 0.0; 201 } else metric[Neq * vv + k] = PetscRealPart(met[dim * dim * v + dim * i + j]); 202 } 203 } 204 } 205 PetscCall(VecRestoreArrayRead(vertexMetric, &met)); 206 207 /* Build ParMmg communicators: the list of vertices between two partitions */ 208 numNgbRanks = 0; 209 if (numProcs > 1) { 210 DM rankdm; 211 PetscSection rankSection, rankGlobalSection; 212 PetscSF rankSF; 213 const PetscInt *degree; 214 PetscInt *rankOfUsedVertices, *rankOfUsedMultiRootLeaves, *usedCopies; 215 PetscInt *rankArray, *rankGlobalArray, *interfacesPerRank; 216 PetscInt offset, mrl, rootDegreeCnt, s, shareCnt, gv; 217 218 PetscCall(PetscSFComputeDegreeBegin(sf, °ree)); 219 PetscCall(PetscSFComputeDegreeEnd(sf, °ree)); 220 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 221 for (i = 0, rootDegreeCnt = 0; i < pEnd - pStart; ++i) rootDegreeCnt += degree[i]; 222 223 /* rankOfUsedVertices, point-array: rank+1 if vertex and in use */ 224 PetscCall(PetscCalloc1(pEnd - pStart, &rankOfUsedVertices)); 225 for (i = 0; i < pEnd - pStart; ++i) rankOfUsedVertices[i] = -1; 226 for (i = vStart; i < vEnd; ++i) { 227 if (vertexNumber[i - vStart]) rankOfUsedVertices[i] = rank; 228 } 229 230 /* rankOfUsedMultiRootLeaves, multiroot-array: rank if vertex and in use, else -1 */ 231 PetscCall(PetscMalloc1(rootDegreeCnt, &rankOfUsedMultiRootLeaves)); 232 PetscCall(PetscSFGatherBegin(sf, MPIU_INT, rankOfUsedVertices, rankOfUsedMultiRootLeaves)); 233 PetscCall(PetscSFGatherEnd(sf, MPIU_INT, rankOfUsedVertices, rankOfUsedMultiRootLeaves)); 234 PetscCall(PetscFree(rankOfUsedVertices)); 235 236 /* usedCopies, point-array: if vertex, shared by how many processes */ 237 PetscCall(PetscCalloc1(pEnd - pStart, &usedCopies)); 238 for (i = 0, mrl = 0; i < vStart - pStart; i++) mrl += degree[i]; 239 for (i = vStart - pStart; i < vEnd - pStart; ++i) { 240 for (j = 0; j < degree[i]; j++, mrl++) { 241 if (rankOfUsedMultiRootLeaves[mrl] != -1) usedCopies[i]++; 242 } 243 if (vertexNumber[i - vStart + pStart]) usedCopies[i]++; 244 } 245 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, usedCopies, usedCopies, MPI_REPLACE)); 246 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, usedCopies, usedCopies, MPI_REPLACE)); 247 248 /* Create a section to store ranks of vertices shared by more than one process */ 249 PetscCall(PetscSectionCreate(comm, &rankSection)); 250 PetscCall(PetscSectionSetNumFields(rankSection, 1)); 251 PetscCall(PetscSectionSetChart(rankSection, pStart, pEnd)); 252 for (i = vStart - pStart; i < vEnd - pStart; ++i) { 253 if (usedCopies[i] > 1) { PetscCall(PetscSectionSetDof(rankSection, i + pStart, usedCopies[i])); } 254 } 255 PetscCall(PetscSectionSetUp(rankSection)); 256 PetscCall(PetscSectionCreateGlobalSection(rankSection, sf, PETSC_FALSE, PETSC_FALSE, PETSC_TRUE, &rankGlobalSection)); 257 258 PetscCall(PetscSectionGetStorageSize(rankGlobalSection, &s)); 259 PetscCall(PetscMalloc1(s, &rankGlobalArray)); 260 for (i = 0, mrl = 0; i < vStart - pStart; i++) mrl += degree[i]; 261 for (i = vStart - pStart, k = 0; i < vEnd - pStart; ++i) { 262 if (usedCopies[i] > 1 && degree[i]) { 263 PetscCall(PetscSectionGetOffset(rankSection, k, &offset)); 264 if (vertexNumber[i - vStart + pStart]) rankGlobalArray[k++] = rank; 265 for (j = 0; j < degree[i]; j++, mrl++) { 266 if (rankOfUsedMultiRootLeaves[mrl] != -1) { rankGlobalArray[k++] = rankOfUsedMultiRootLeaves[mrl]; } 267 } 268 } else mrl += degree[i]; 269 } 270 PetscCall(PetscFree(rankOfUsedMultiRootLeaves)); 271 PetscCall(PetscFree(usedCopies)); 272 PetscCall(PetscSectionDestroy(&rankGlobalSection)); 273 274 /* 275 Broadcast the array of ranks. 276 (We want all processes to know all the ranks that are looking at each point. 277 Above, we tell the roots. Here, the roots tell the leaves.) 278 */ 279 PetscCall(DMClone(dm, &rankdm)); 280 PetscCall(DMSetLocalSection(rankdm, rankSection)); 281 PetscCall(DMGetSectionSF(rankdm, &rankSF)); 282 PetscCall(PetscSectionGetStorageSize(rankSection, &s)); 283 PetscCall(PetscMalloc1(s, &rankArray)); 284 PetscCall(PetscSFBcastBegin(rankSF, MPI_INT, rankGlobalArray, rankArray, MPI_REPLACE)); 285 PetscCall(PetscSFBcastEnd(rankSF, MPI_INT, rankGlobalArray, rankArray, MPI_REPLACE)); 286 PetscCall(PetscFree(rankGlobalArray)); 287 PetscCall(DMDestroy(&rankdm)); 288 289 /* Count the number of interfaces per rank, not including those on the root */ 290 PetscCall(PetscCalloc1(numProcs, &interfacesPerRank)); 291 for (v = vStart; v < vEnd; v++) { 292 if (vertexNumber[v - vStart]) { 293 PetscCall(PetscSectionGetDof(rankSection, v, &shareCnt)); 294 if (shareCnt) { 295 PetscCall(PetscSectionGetOffset(rankSection, v, &offset)); 296 for (j = 0; j < shareCnt; j++) { interfacesPerRank[rankArray[offset + j]]++; } 297 } 298 } 299 } 300 for (r = 0, k = 0, interfacesPerRank[rank] = 0; r < numProcs; r++) k += interfacesPerRank[r]; 301 302 /* Get the degree of the vertex */ 303 PetscCall(PetscMalloc3(k, &interfaces_lv, k, &interfaces_gv, numProcs + 1, &interfacesOffset)); 304 interfacesOffset[0] = 0; 305 for (r = 0; r < numProcs; r++) { 306 interfacesOffset[r + 1] = interfacesOffset[r] + interfacesPerRank[r]; 307 if (interfacesPerRank[r]) numNgbRanks++; 308 interfacesPerRank[r] = 0; 309 } 310 311 /* Get the local and global vertex numbers at interfaces */ 312 PetscCall(DMPlexGetVertexNumbering(dm, &globalVertexNum)); 313 PetscCall(ISGetIndices(globalVertexNum, &gV)); 314 for (v = vStart; v < vEnd; v++) { 315 if (vertexNumber[v - vStart]) { 316 PetscCall(PetscSectionGetDof(rankSection, v, &shareCnt)); 317 if (shareCnt) { 318 PetscCall(PetscSectionGetOffset(rankSection, v, &offset)); 319 for (j = 0; j < shareCnt; j++) { 320 r = rankArray[offset + j]; 321 if (r == rank) continue; 322 k = interfacesOffset[r] + interfacesPerRank[r]++; 323 interfaces_lv[k] = vertexNumber[v - vStart]; 324 gv = gV[v - vStart]; 325 interfaces_gv[k] = gv < 0 ? -gv : gv + 1; 326 } 327 } 328 } 329 } 330 PetscCall(PetscFree(interfacesPerRank)); 331 PetscCall(PetscFree(rankArray)); 332 PetscCall(ISRestoreIndices(globalVertexNum, &gV)); 333 PetscCall(PetscSectionDestroy(&rankSection)); 334 } 335 PetscCall(DMDestroy(&udm)); 336 PetscCall(PetscFree(vertexNumber)); 337 338 /* Send the data to ParMmg and remesh */ 339 PetscCall(DMPlexMetricNoInsertion(dm, &noInsert)); 340 PetscCall(DMPlexMetricNoSwapping(dm, &noSwap)); 341 PetscCall(DMPlexMetricNoMovement(dm, &noMove)); 342 PetscCall(DMPlexMetricNoSurf(dm, &noSurf)); 343 PetscCall(DMPlexMetricGetVerbosity(dm, &verbosity)); 344 PetscCall(DMPlexMetricGetNumIterations(dm, &numIter)); 345 PetscCall(DMPlexMetricGetGradationFactor(dm, &gradationFactor)); 346 PetscCall(DMPlexMetricGetHausdorffNumber(dm, &hausdorffNumber)); 347 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); 348 PetscCallMMG_NONSTANDARD(PMMG_Set_meshSize, parmesh, numUsedVertices, numCellsNotShared, 0, numFaceTags, 0, 0); 349 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_APImode, PMMG_APIDISTRIB_nodes); 350 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_noinsert, noInsert); 351 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_noswap, noSwap); 352 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_nomove, noMove); 353 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_nosurf, noSurf); 354 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_verbose, verbosity); 355 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_globalNum, 1); 356 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_niter, numIter); 357 PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter, parmesh, PMMG_DPARAM_hgrad, gradationFactor); 358 PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter, parmesh, PMMG_DPARAM_hausd, hausdorffNumber); 359 PetscCallMMG_NONSTANDARD(PMMG_Set_vertices, parmesh, vertices, verTags); 360 PetscCallMMG_NONSTANDARD(PMMG_Set_tetrahedra, parmesh, cells, cellTags); 361 PetscCallMMG_NONSTANDARD(PMMG_Set_triangles, parmesh, bdFaces, faceTags); 362 PetscCallMMG_NONSTANDARD(PMMG_Set_metSize, parmesh, MMG5_Vertex, numUsedVertices, MMG5_Tensor); 363 PetscCallMMG_NONSTANDARD(PMMG_Set_tensorMets, parmesh, metric); 364 PetscCallMMG_NONSTANDARD(PMMG_Set_numberOfNodeCommunicators, parmesh, numNgbRanks); 365 for (r = 0, c = 0; r < numProcs; ++r) { 366 if (interfacesOffset[r + 1] > interfacesOffset[r]) { 367 PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicatorSize, parmesh, c, r, interfacesOffset[r + 1] - interfacesOffset[r]); 368 PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicator_nodes, parmesh, c++, &interfaces_lv[interfacesOffset[r]], &interfaces_gv[interfacesOffset[r]], 1); 369 } 370 } 371 PetscCallMMG(PMMG_parmmglib_distributed, parmesh); 372 PetscCall(PetscFree(cells)); 373 PetscCall(PetscFree2(metric, vertices)); 374 PetscCall(PetscFree2(bdFaces, faceTags)); 375 PetscCall(PetscFree2(verTags, cellTags)); 376 if (numProcs > 1) { PetscCall(PetscFree3(interfaces_lv, interfaces_gv, interfacesOffset)); } 377 378 /* Retrieve mesh from Mmg */ 379 numCornersNew = 4; 380 PetscCallMMG_NONSTANDARD(PMMG_Get_meshSize, parmesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0); 381 PetscCall(PetscMalloc4(dim * numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer)); 382 PetscCall(PetscMalloc3((dim + 1) * numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells)); 383 PetscCall(PetscMalloc4(dim * numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces)); 384 PetscCallMMG_NONSTANDARD(PMMG_Get_vertices, parmesh, verticesNew, verTagsNew, corners, requiredVer); 385 PetscCallMMG_NONSTANDARD(PMMG_Get_tetrahedra, parmesh, cellsNew, cellTagsNew, requiredCells); 386 PetscCallMMG_NONSTANDARD(PMMG_Get_triangles, parmesh, facesNew, faceTagsNew, requiredFaces); 387 PetscCall(PetscMalloc2(numVerticesNew, &owners, numVerticesNew, &gv_new)); 388 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_globalNum, 1); 389 PetscCallMMG_NONSTANDARD(PMMG_Get_verticesGloNum, parmesh, gv_new, owners); 390 for (i = 0; i < dim * numFacesNew; ++i) facesNew[i] -= 1; 391 for (i = 0; i < (dim + 1) * numCellsNew; ++i) cellsNew[i] = gv_new[cellsNew[i] - 1] - 1; 392 for (i = 0, numVerticesNewLoc = 0; i < numVerticesNew; ++i) { 393 if (owners[i] == rank) numVerticesNewLoc++; 394 } 395 PetscCall(PetscMalloc2(numVerticesNewLoc * dim, &verticesNewLoc, numVerticesNew, &verticesNewSorted)); 396 for (i = 0, c = 0; i < numVerticesNew; i++) { 397 if (owners[i] == rank) { 398 for (j = 0; j < dim; ++j) verticesNewLoc[dim * c + j] = verticesNew[dim * i + j]; 399 c++; 400 } 401 } 402 403 /* Reorder for consistency with DMPlex */ 404 for (i = 0; i < numCellsNew; ++i) PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4 * i])); 405 406 /* Create new plex */ 407 PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNewLoc, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNewLoc, NULL, &verticesNewSorted, dmNew)); 408 PetscCallMMG_NONSTANDARD(PMMG_Free_all, PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_end); 409 PetscCall(PetscFree4(verticesNew, verTagsNew, corners, requiredVer)); 410 411 /* Get adapted mesh information */ 412 PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd)); 413 PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd)); 414 PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd)); 415 416 /* Rebuild boundary label */ 417 PetscCall(DMCreateLabel(*dmNew, flg ? bdName : bdLabelName)); 418 PetscCall(DMGetLabel(*dmNew, flg ? bdName : bdLabelName, &bdLabelNew)); 419 for (i = 0; i < numFacesNew; i++) { 420 PetscBool hasTag = PETSC_FALSE; 421 PetscInt numCoveredPoints, numFaces = 0, facePoints[3]; 422 const PetscInt *coveredPoints = NULL; 423 424 for (j = 0; j < dim; ++j) { 425 lv = facesNew[i * dim + j]; 426 gv = gv_new[lv] - 1; 427 PetscCall(PetscFindInt(gv, numVerticesNew, verticesNewSorted, &lv)); 428 facePoints[j] = lv + vStart; 429 } 430 PetscCall(DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints)); 431 for (j = 0; j < numCoveredPoints; ++j) { 432 if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) { 433 numFaces++; 434 f = j; 435 } 436 } 437 PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " vertices cannot define more than 1 facet (%" PetscInt_FMT ")", dim, numFaces); 438 PetscCall(DMLabelHasStratum(bdLabel, faceTagsNew[i], &hasTag)); 439 if (hasTag) PetscCall(DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i])); 440 PetscCall(DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints)); 441 } 442 PetscCall(PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces)); 443 PetscCall(PetscFree2(owners, gv_new)); 444 PetscCall(PetscFree2(verticesNewLoc, verticesNewSorted)); 445 if (flg) PetscCall(DMLabelDestroy(&bdLabel)); 446 447 /* Rebuild cell labels */ 448 PetscCall(DMCreateLabel(*dmNew, rgLabel ? rgLabelName : rgName)); 449 PetscCall(DMGetLabel(*dmNew, rgLabel ? rgLabelName : rgName, &rgLabelNew)); 450 for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(rgLabelNew, c, cellTagsNew[c - cStart])); 451 PetscCall(PetscFree3(cellsNew, cellTagsNew, requiredCells)); 452 PetscFunctionReturn(PETSC_SUCCESS); 453 } 454