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