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