xref: /petsc/src/dm/impls/plex/adaptors/mmg/mmgadapt.c (revision 970231d20df44f79b27787157e39d441e79f434b) 
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 
DMAdaptMetric_Mmg_Plex(DM dm,Vec vertexMetric,DMLabel bdLabel,DMLabel rgLabel,DM * dmNew)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   PetscFunctionReturn(PETSC_SUCCESS);
280 }
281