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