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