1 #include "../mmgcommon.h" /*I "petscdmplex.h" I*/ 2 #include <mmg/libmmg.h> 3 4 PetscBool MmgCite = PETSC_FALSE; 5 const char MmgCitation[] = "@article{DAPOGNY2014358,\n" 6 " title = {Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems},\n" 7 " journal = {Journal of Computational Physics},\n" 8 " author = {C. Dapogny and C. Dobrzynski and P. Frey},\n" 9 " volume = {262},\n" 10 " pages = {358--378},\n" 11 " doi = {10.1016/j.jcp.2014.01.005},\n" 12 " year = {2014}\n}\n"; 13 14 PETSC_EXTERN PetscErrorCode DMAdaptMetric_Mmg_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DMLabel rgLabel, DM *dmNew) 15 { 16 MPI_Comm comm; 17 const char *bdName = "_boundary_"; 18 const char *rgName = "_regions_"; 19 DM udm, cdm; 20 DMLabel bdLabelNew, rgLabelNew; 21 const char *bdLabelName, *rgLabelName; 22 PetscSection coordSection; 23 Vec coordinates; 24 const PetscScalar *coords, *met; 25 PetscReal *vertices, *metric, *verticesNew, gradationFactor, hausdorffNumber; 26 PetscInt *cells, *cellsNew, *cellTags, *cellTagsNew, *verTags, *verTagsNew; 27 PetscInt *bdFaces, *faceTags, *facesNew, *faceTagsNew; 28 int *corners, *requiredCells, *requiredVer, *ridges, *requiredFaces; 29 PetscInt cStart, cEnd, c, numCells, fStart, fEnd, numFaceTags, f, vStart, vEnd, v, numVertices; 30 PetscInt dim, off, coff, maxConeSize, bdSize, i, j, k, Neq, verbosity, pStart, pEnd; 31 PetscInt numCellsNew, numVerticesNew, numCornersNew, numFacesNew; 32 PetscBool flg = PETSC_FALSE, noInsert, noSwap, noMove, noSurf, isotropic, uniform; 33 MMG5_pMesh mmg_mesh = NULL; 34 MMG5_pSol mmg_metric = NULL; 35 36 PetscFunctionBegin; 37 PetscCall(PetscCitationsRegister(MmgCitation, &MmgCite)); 38 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 39 if (bdLabel) { 40 PetscCall(PetscObjectGetName((PetscObject)bdLabel, &bdLabelName)); 41 PetscCall(PetscStrcmp(bdLabelName, bdName, &flg)); 42 PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 43 } 44 if (rgLabel) { 45 PetscCall(PetscObjectGetName((PetscObject)rgLabel, &rgLabelName)); 46 PetscCall(PetscStrcmp(rgLabelName, rgName, &flg)); 47 PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for element tags", rgLabelName); 48 } 49 50 /* Get mesh information */ 51 PetscCall(DMGetDimension(dm, &dim)); 52 Neq = (dim * (dim + 1)) / 2; 53 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 54 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 55 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 56 PetscCall(DMPlexUninterpolate(dm, &udm)); 57 PetscCall(DMPlexGetMaxSizes(udm, &maxConeSize, NULL)); 58 numCells = cEnd - cStart; 59 numVertices = vEnd - vStart; 60 61 /* Get cell offsets */ 62 PetscCall(PetscMalloc1(numCells * maxConeSize, &cells)); 63 for (c = 0, coff = 0; c < numCells; ++c) { 64 const PetscInt *cone; 65 PetscInt coneSize, cl; 66 67 PetscCall(DMPlexGetConeSize(udm, c, &coneSize)); 68 PetscCall(DMPlexGetCone(udm, c, &cone)); 69 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart + 1; 70 } 71 72 /* Get vertex coordinate array */ 73 PetscCall(DMGetCoordinateDM(dm, &cdm)); 74 PetscCall(DMGetLocalSection(cdm, &coordSection)); 75 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 76 PetscCall(VecGetArrayRead(coordinates, &coords)); 77 PetscCall(PetscMalloc2(numVertices * Neq, &metric, dim * numVertices, &vertices)); 78 for (v = 0; v < vEnd - vStart; ++v) { 79 PetscCall(PetscSectionGetOffset(coordSection, v + vStart, &off)); 80 for (i = 0; i < dim; ++i) vertices[dim * v + i] = PetscRealPart(coords[off + i]); 81 } 82 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 83 PetscCall(DMDestroy(&udm)); 84 85 /* Get face tags */ 86 if (!bdLabel) { 87 flg = PETSC_TRUE; 88 PetscCall(DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabel)); 89 PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabel)); 90 } 91 PetscCall(DMLabelGetBounds(bdLabel, &pStart, &pEnd)); 92 for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) { 93 PetscBool hasPoint; 94 PetscInt *closure = NULL, closureSize, cl; 95 96 PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint)); 97 if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue; 98 numFaceTags++; 99 100 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 101 for (cl = 0; cl < closureSize * 2; cl += 2) { 102 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 103 } 104 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 105 } 106 PetscCall(PetscMalloc2(bdSize, &bdFaces, numFaceTags, &faceTags)); 107 for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) { 108 PetscBool hasPoint; 109 PetscInt *closure = NULL, closureSize, cl; 110 111 PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint)); 112 if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue; 113 114 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 115 for (cl = 0; cl < closureSize * 2; cl += 2) { 116 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart + 1; 117 } 118 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure)); 119 PetscCall(DMLabelGetValue(bdLabel, f, &faceTags[numFaceTags++])); 120 } 121 if (flg) PetscCall(DMLabelDestroy(&bdLabel)); 122 123 /* Get cell tags */ 124 PetscCall(PetscCalloc2(numVertices, &verTags, numCells, &cellTags)); 125 if (rgLabel) { 126 for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelGetValue(rgLabel, c, &cellTags[c])); 127 } 128 129 /* Get metric */ 130 PetscCall(VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view")); 131 PetscCall(VecGetArrayRead(vertexMetric, &met)); 132 PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic)); 133 PetscCall(DMPlexMetricIsUniform(dm, &uniform)); 134 for (v = 0; v < (vEnd - vStart); ++v) { 135 for (i = 0, k = 0; i < dim; ++i) { 136 for (j = i; j < dim; ++j) { 137 if (isotropic) { 138 if (i == j) { 139 if (uniform) metric[Neq * v + k] = PetscRealPart(met[0]); 140 else metric[Neq * v + k] = PetscRealPart(met[v]); 141 } else metric[Neq * v + k] = 0.0; 142 } else { 143 metric[Neq * v + k] = PetscRealPart(met[dim * dim * v + dim * i + j]); 144 } 145 k++; 146 } 147 } 148 } 149 PetscCall(VecRestoreArrayRead(vertexMetric, &met)); 150 151 /* Send mesh to Mmg and remesh */ 152 PetscCall(DMPlexMetricGetVerbosity(dm, &verbosity)); 153 PetscCall(DMPlexMetricGetGradationFactor(dm, &gradationFactor)); 154 PetscCall(DMPlexMetricGetHausdorffNumber(dm, &hausdorffNumber)); 155 PetscCall(DMPlexMetricNoInsertion(dm, &noInsert)); 156 PetscCall(DMPlexMetricNoSwapping(dm, &noSwap)); 157 PetscCall(DMPlexMetricNoMovement(dm, &noMove)); 158 PetscCall(DMPlexMetricNoSurf(dm, &noSurf)); 159 switch (dim) { 160 case 2: 161 PetscCallMMG_NONSTANDARD(MMG2D_Init_mesh, MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 162 PetscCallMMG_NONSTANDARD(MMG2D_Set_iparameter, mmg_mesh, mmg_metric, MMG2D_IPARAM_noinsert, noInsert); 163 PetscCallMMG_NONSTANDARD(MMG2D_Set_iparameter, mmg_mesh, mmg_metric, MMG2D_IPARAM_noswap, noSwap); 164 PetscCallMMG_NONSTANDARD(MMG2D_Set_iparameter, mmg_mesh, mmg_metric, MMG2D_IPARAM_nomove, noMove); 165 PetscCallMMG_NONSTANDARD(MMG2D_Set_iparameter, mmg_mesh, mmg_metric, MMG2D_IPARAM_nosurf, noSurf); 166 PetscCallMMG_NONSTANDARD(MMG2D_Set_iparameter, mmg_mesh, mmg_metric, MMG2D_IPARAM_verbose, verbosity); 167 PetscCallMMG_NONSTANDARD(MMG2D_Set_dparameter, mmg_mesh, mmg_metric, MMG2D_DPARAM_hgrad, gradationFactor); 168 PetscCallMMG_NONSTANDARD(MMG2D_Set_dparameter, mmg_mesh, mmg_metric, MMG2D_DPARAM_hausd, hausdorffNumber); 169 PetscCallMMG_NONSTANDARD(MMG2D_Set_meshSize, mmg_mesh, numVertices, numCells, 0, numFaceTags); 170 PetscCallMMG_NONSTANDARD(MMG2D_Set_vertices, mmg_mesh, vertices, verTags); 171 PetscCallMMG_NONSTANDARD(MMG2D_Set_triangles, mmg_mesh, cells, cellTags); 172 PetscCallMMG_NONSTANDARD(MMG2D_Set_edges, mmg_mesh, bdFaces, faceTags); 173 PetscCallMMG_NONSTANDARD(MMG2D_Set_solSize, mmg_mesh, mmg_metric, MMG5_Vertex, numVertices, MMG5_Tensor); 174 PetscCallMMG_NONSTANDARD(MMG2D_Set_tensorSols, mmg_metric, metric); 175 PetscCallMMG(MMG2D_mmg2dlib, mmg_mesh, mmg_metric); 176 break; 177 case 3: 178 PetscCallMMG_NONSTANDARD(MMG3D_Init_mesh, MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 179 PetscCallMMG_NONSTANDARD(MMG3D_Set_iparameter, mmg_mesh, mmg_metric, MMG3D_IPARAM_noinsert, noInsert); 180 PetscCallMMG_NONSTANDARD(MMG3D_Set_iparameter, mmg_mesh, mmg_metric, MMG3D_IPARAM_noswap, noSwap); 181 PetscCallMMG_NONSTANDARD(MMG3D_Set_iparameter, mmg_mesh, mmg_metric, MMG3D_IPARAM_nomove, noMove); 182 PetscCallMMG_NONSTANDARD(MMG3D_Set_iparameter, mmg_mesh, mmg_metric, MMG3D_IPARAM_nosurf, noSurf); 183 PetscCallMMG_NONSTANDARD(MMG3D_Set_iparameter, mmg_mesh, mmg_metric, MMG3D_IPARAM_verbose, verbosity); 184 PetscCallMMG_NONSTANDARD(MMG3D_Set_dparameter, mmg_mesh, mmg_metric, MMG3D_DPARAM_hgrad, gradationFactor); 185 PetscCallMMG_NONSTANDARD(MMG3D_Set_dparameter, mmg_mesh, mmg_metric, MMG3D_DPARAM_hausd, hausdorffNumber); 186 PetscCallMMG_NONSTANDARD(MMG3D_Set_meshSize, mmg_mesh, numVertices, numCells, 0, numFaceTags, 0, 0); 187 PetscCallMMG_NONSTANDARD(MMG3D_Set_vertices, mmg_mesh, vertices, verTags); 188 PetscCallMMG_NONSTANDARD(MMG3D_Set_tetrahedra, mmg_mesh, cells, cellTags); 189 PetscCallMMG_NONSTANDARD(MMG3D_Set_triangles, mmg_mesh, bdFaces, faceTags); 190 PetscCallMMG_NONSTANDARD(MMG3D_Set_solSize, mmg_mesh, mmg_metric, MMG5_Vertex, numVertices, MMG5_Tensor); 191 PetscCallMMG_NONSTANDARD(MMG3D_Set_tensorSols, mmg_metric, metric); 192 PetscCallMMG(MMG3D_mmg3dlib, mmg_mesh, mmg_metric); 193 break; 194 default: 195 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Mmg adaptation defined for dimension %" PetscInt_FMT, dim); 196 } 197 PetscCall(PetscFree(cells)); 198 PetscCall(PetscFree2(metric, vertices)); 199 PetscCall(PetscFree2(bdFaces, faceTags)); 200 PetscCall(PetscFree2(verTags, cellTags)); 201 202 /* Retrieve mesh from Mmg */ 203 switch (dim) { 204 case 2: 205 numCornersNew = 3; 206 PetscCallMMG_NONSTANDARD(MMG2D_Get_meshSize, mmg_mesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew); 207 PetscCall(PetscMalloc4(2 * numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer)); 208 PetscCall(PetscMalloc3(3 * numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells)); 209 PetscCall(PetscMalloc4(2 * numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces)); 210 PetscCallMMG_NONSTANDARD(MMG2D_Get_vertices, mmg_mesh, verticesNew, verTagsNew, corners, requiredVer); 211 PetscCallMMG_NONSTANDARD(MMG2D_Get_triangles, mmg_mesh, cellsNew, cellTagsNew, requiredCells); 212 PetscCallMMG_NONSTANDARD(MMG2D_Get_edges, mmg_mesh, facesNew, faceTagsNew, ridges, requiredFaces); 213 break; 214 case 3: 215 numCornersNew = 4; 216 PetscCallMMG_NONSTANDARD(MMG3D_Get_meshSize, mmg_mesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0); 217 PetscCall(PetscMalloc4(3 * numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer)); 218 PetscCall(PetscMalloc3(4 * numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells)); 219 PetscCall(PetscMalloc4(3 * numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces)); 220 PetscCallMMG_NONSTANDARD(MMG3D_Get_vertices, mmg_mesh, verticesNew, verTagsNew, corners, requiredVer); 221 PetscCallMMG_NONSTANDARD(MMG3D_Get_tetrahedra, mmg_mesh, cellsNew, cellTagsNew, requiredCells); 222 PetscCallMMG_NONSTANDARD(MMG3D_Get_triangles, mmg_mesh, facesNew, faceTagsNew, requiredFaces); 223 224 /* Reorder for consistency with DMPlex */ 225 for (i = 0; i < numCellsNew; ++i) PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4 * i])); 226 break; 227 228 default: 229 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Mmg adaptation defined for dimension %" PetscInt_FMT, dim); 230 } 231 232 /* Create new Plex */ 233 for (i = 0; i < (dim + 1) * numCellsNew; i++) cellsNew[i] -= 1; 234 for (i = 0; i < dim * numFacesNew; i++) facesNew[i] -= 1; 235 PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNew, NULL, NULL, dmNew)); 236 switch (dim) { 237 case 2: 238 PetscCallMMG_NONSTANDARD(MMG2D_Free_all, MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 239 break; 240 case 3: 241 PetscCallMMG_NONSTANDARD(MMG3D_Free_all, MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 242 break; 243 default: 244 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Mmg adaptation defined for dimension %" PetscInt_FMT, dim); 245 } 246 PetscCall(PetscFree4(verticesNew, verTagsNew, corners, requiredVer)); 247 248 /* Get adapted mesh information */ 249 PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd)); 250 PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd)); 251 PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd)); 252 253 /* Rebuild boundary labels */ 254 PetscCall(DMCreateLabel(*dmNew, flg ? bdName : bdLabelName)); 255 PetscCall(DMGetLabel(*dmNew, flg ? bdName : bdLabelName, &bdLabelNew)); 256 for (i = 0; i < numFacesNew; i++) { 257 PetscInt numCoveredPoints, numFaces = 0, facePoints[3]; 258 const PetscInt *coveredPoints = NULL; 259 260 for (j = 0; j < dim; ++j) facePoints[j] = facesNew[i * dim + j] + vStart; 261 PetscCall(DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints)); 262 for (j = 0; j < numCoveredPoints; ++j) { 263 if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) { 264 numFaces++; 265 f = j; 266 } 267 } 268 PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " vertices cannot define more than 1 facet (%" PetscInt_FMT ")", dim, numFaces); 269 PetscCall(DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i])); 270 PetscCall(DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints)); 271 } 272 PetscCall(PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces)); 273 274 /* Rebuild cell labels */ 275 PetscCall(DMCreateLabel(*dmNew, rgLabel ? rgLabelName : rgName)); 276 PetscCall(DMGetLabel(*dmNew, rgLabel ? rgLabelName : rgName, &rgLabelNew)); 277 for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(rgLabelNew, c, cellTagsNew[c - cStart])); 278 PetscCall(PetscFree3(cellsNew, cellTagsNew, requiredCells)); 279 280 PetscFunctionReturn(PETSC_SUCCESS); 281 } 282