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, DM *dmNew) 5 { 6 MPI_Comm comm; 7 const char *bdName = "_boundary_"; 8 DM udm, cdm; 9 DMLabel bdLabelFull; 10 const char *bdLabelName; 11 IS bdIS; 12 PetscSection coordSection; 13 Vec coordinates; 14 const PetscScalar *coords, *met; 15 const PetscInt *bdFacesFull; 16 PetscInt *bdFaces, *bdFaceIds; 17 PetscReal *vertices, *metric, *verticesNew, gradationFactor; 18 PetscInt *cells, *cellsNew; 19 PetscInt *facesNew, *faceTagsNew; 20 PetscInt *verTags, *cellTags, *verTagsNew, *cellTagsNew; 21 PetscInt *corners, *requiredCells, *requiredVer, *ridges, *requiredFaces; 22 PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, v; 23 PetscInt off, maxConeSize, numBdFaces, f, bdSize, fStart, fEnd, i, j, k, Neq; 24 PetscInt numCellsNew, numVerticesNew, numCornersNew, numFacesNew; 25 PetscInt verbosity; 26 PetscBool flg, noInsert, noSwap, noMove; 27 DMLabel bdLabelNew; 28 MMG5_pMesh mmg_mesh = NULL; 29 MMG5_pSol mmg_metric = NULL; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 34 if (bdLabel) { 35 ierr = PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);CHKERRQ(ierr); 36 ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr); 37 if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 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 = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 45 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 46 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); 47 numCells = cEnd - cStart; 48 numVertices = vEnd - vStart; 49 50 /* Get cell offsets */ 51 ierr = PetscMalloc1(numCells*maxConeSize, &cells);CHKERRQ(ierr); 52 for (c = 0, coff = 0; c < numCells; ++c) { 53 const PetscInt *cone; 54 PetscInt coneSize, cl; 55 56 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 57 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 58 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart + 1; 59 } 60 61 /* Get vertex coordinate array */ 62 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 63 ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr); 64 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 65 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 66 ierr = PetscMalloc2(numVertices*Neq, &metric, dim*numVertices, &vertices);CHKERRQ(ierr); 67 for (v = 0; v < vEnd-vStart; ++v) { 68 ierr = PetscSectionGetOffset(coordSection, v+vStart, &off);CHKERRQ(ierr); 69 for (i = 0; i < dim; ++i) vertices[dim*v+i] = PetscRealPart(coords[off+i]); 70 } 71 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 72 ierr = DMDestroy(&udm);CHKERRQ(ierr); 73 74 /* Get boundary mesh */ 75 ierr = DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);CHKERRQ(ierr); 76 ierr = DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);CHKERRQ(ierr); 77 ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr); 78 ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr); 79 ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr); 80 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 81 PetscInt *closure = NULL; 82 PetscInt closureSize, cl; 83 84 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 85 for (cl = 0; cl < closureSize*2; cl += 2) { 86 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 87 } 88 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 89 } 90 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 91 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 92 PetscInt *closure = NULL; 93 PetscInt closureSize, cl; 94 95 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 96 for (cl = 0; cl < closureSize*2; cl += 2) { 97 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart + 1; 98 } 99 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 100 if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);} 101 else {bdFaceIds[f] = 1;} 102 } 103 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 104 ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr); 105 106 /* Get metric */ 107 ierr = VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");CHKERRQ(ierr); 108 ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr); 109 for (v = 0; v < (vEnd-vStart); ++v) { 110 for (i = 0, k = 0; i < dim; ++i) { 111 for (j = i; j < dim; ++j) { 112 metric[Neq*v+k] = PetscRealPart(met[dim*dim*v+dim*i+j]); 113 k++; 114 } 115 } 116 } 117 ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr); 118 119 /* Send mesh to Mmg and remesh */ 120 ierr = DMPlexMetricGetVerbosity(dm, &verbosity);CHKERRQ(ierr); 121 ierr = DMPlexMetricGetGradationFactor(dm, &gradationFactor);CHKERRQ(ierr); 122 ierr = DMPlexMetricNoInsertion(dm, &noInsert);CHKERRQ(ierr); 123 ierr = DMPlexMetricNoSwapping(dm, &noSwap);CHKERRQ(ierr); 124 ierr = DMPlexMetricNoMovement(dm, &noMove);CHKERRQ(ierr); 125 ierr = PetscCalloc2(numVertices, &verTags, numCells, &cellTags);CHKERRQ(ierr); 126 switch (dim) { 127 case 2: 128 ierr = MMG2D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 129 ierr = MMG2D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_noinsert, noInsert); 130 ierr = MMG2D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_noswap, noSwap); 131 ierr = MMG2D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_nomove, noMove); 132 ierr = MMG2D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_verbose, verbosity); 133 ierr = MMG2D_Set_dparameter(mmg_mesh, mmg_metric, MMG2D_DPARAM_hgrad, gradationFactor); 134 ierr = MMG2D_Set_meshSize(mmg_mesh, numVertices, numCells, 0, numBdFaces); 135 ierr = MMG2D_Set_vertices(mmg_mesh, vertices, verTags); 136 ierr = MMG2D_Set_triangles(mmg_mesh, cells, cellTags); 137 ierr = MMG2D_Set_edges(mmg_mesh, bdFaces, bdFaceIds); 138 ierr = MMG2D_Set_solSize(mmg_mesh, mmg_metric, MMG5_Vertex, numVertices, MMG5_Tensor); 139 ierr = MMG2D_Set_tensorSols(mmg_metric, metric); 140 ierr = MMG2D_mmg2dlib(mmg_mesh, mmg_metric); 141 break; 142 case 3: 143 ierr = MMG3D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 144 ierr = MMG3D_Set_iparameter(mmg_mesh, mmg_metric, MMG3D_IPARAM_noinsert, noInsert); 145 ierr = MMG3D_Set_iparameter(mmg_mesh, mmg_metric, MMG3D_IPARAM_noswap, noSwap); 146 ierr = MMG3D_Set_iparameter(mmg_mesh, mmg_metric, MMG3D_IPARAM_nomove, noMove); 147 ierr = MMG3D_Set_iparameter(mmg_mesh, mmg_metric, MMG3D_IPARAM_verbose, verbosity); 148 ierr = MMG3D_Set_dparameter(mmg_mesh, mmg_metric, MMG3D_DPARAM_hgrad, gradationFactor); 149 ierr = MMG3D_Set_meshSize(mmg_mesh, numVertices, numCells, 0, numBdFaces, 0, 0); 150 ierr = MMG3D_Set_vertices(mmg_mesh, vertices, verTags); 151 ierr = MMG3D_Set_tetrahedra(mmg_mesh, cells, cellTags); 152 ierr = MMG3D_Set_triangles(mmg_mesh, bdFaces, bdFaceIds); 153 ierr = MMG3D_Set_solSize(mmg_mesh, mmg_metric, MMG5_Vertex, numVertices, MMG5_Tensor); 154 ierr = MMG3D_Set_tensorSols(mmg_metric, metric); 155 ierr = MMG3D_mmg3dlib(mmg_mesh, mmg_metric); 156 break; 157 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Mmg adaptation defined for dimension %D", dim); 158 } 159 160 /* Retrieve mesh from Mmg and create new Plex*/ 161 switch (dim) { 162 case 2: 163 numCornersNew = 3; 164 ierr = MMG2D_Get_meshSize(mmg_mesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew); 165 ierr = PetscMalloc4(2*numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer);CHKERRQ(ierr); 166 ierr = PetscMalloc3(3*numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells);CHKERRQ(ierr); 167 ierr = PetscMalloc4(2*numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces);CHKERRQ(ierr); 168 ierr = MMG2D_Get_vertices(mmg_mesh, verticesNew, verTagsNew, corners, requiredVer); 169 ierr = MMG2D_Get_triangles(mmg_mesh, cellsNew, cellTagsNew, requiredCells); 170 ierr = MMG2D_Get_edges(mmg_mesh, facesNew, faceTagsNew, ridges, requiredFaces); 171 break; 172 case 3: 173 numCornersNew = 4; 174 ierr = MMG3D_Get_meshSize(mmg_mesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0); 175 ierr = PetscMalloc4(3*numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer);CHKERRQ(ierr); 176 ierr = PetscMalloc3(4*numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells);CHKERRQ(ierr); 177 ierr = PetscMalloc4(3*numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces);CHKERRQ(ierr); 178 ierr = MMG3D_Get_vertices(mmg_mesh, verticesNew, verTagsNew, corners, requiredVer); 179 ierr = MMG3D_Get_tetrahedra(mmg_mesh, cellsNew, cellTagsNew, requiredCells); 180 ierr = MMG3D_Get_triangles(mmg_mesh, facesNew, faceTagsNew, requiredFaces); 181 break; 182 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Mmg adaptation defined for dimension %D", dim); 183 } 184 for (i = 0; i < (dim+1)*numCellsNew; i++) cellsNew[i] -= 1; 185 for (i = 0; i < dim*numFacesNew; i++) facesNew[i] -= 1; 186 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNew, NULL, NULL, dmNew);CHKERRQ(ierr); 187 switch (dim) { 188 case 2: 189 ierr = MMG2D_Free_all(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 190 break; 191 case 3: 192 ierr = MMG3D_Free_all(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 193 break; 194 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Mmg adaptation defined for dimension %D", dim); 195 } 196 197 /* Rebuild boundary labels */ 198 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 199 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 200 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 201 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr); 202 for (i = 0; i < numFacesNew; i++) { 203 PetscInt numCoveredPoints, numFaces = 0, facePoints[3]; 204 const PetscInt *coveredPoints = NULL; 205 206 for (j = 0; j < dim; ++j) facePoints[j] = facesNew[i*dim+j]+vStart; 207 ierr = DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr); 208 for (j = 0; j < numCoveredPoints; ++j) { 209 if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) { 210 numFaces++; 211 f = j; 212 } 213 } 214 if (numFaces != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "%d vertices cannot define more than 1 facet (%d)", dim, numFaces); 215 ierr = DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i]);CHKERRQ(ierr); 216 ierr = DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr); 217 } 218 219 /* Clean up */ 220 ierr = PetscFree(cells);CHKERRQ(ierr); 221 ierr = PetscFree2(metric, vertices);CHKERRQ(ierr); 222 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 223 ierr = PetscFree2(verTags, cellTags);CHKERRQ(ierr); 224 ierr = PetscFree4(verticesNew, verTagsNew, corners, requiredVer);CHKERRQ(ierr); 225 ierr = PetscFree3(cellsNew, cellTagsNew, requiredCells);CHKERRQ(ierr); 226 ierr = PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229