Actual source code: mmgadapt.c
1: #include "../mmgcommon.h" /*I "petscdmplex.h" I*/
2: #include <mmg/libmmg.h>
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";
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;
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: }
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;
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;
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: }
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));
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;
96: PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint));
97: if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;
98: numFaceTags++;
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;
111: PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint));
112: if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;
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));
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: }
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));
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));
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);
224: /* Reorder for consistency with DMPlex */
225: for (i = 0; i < numCellsNew; ++i) PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4 * i]));
226: break;
228: default:
229: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Mmg adaptation defined for dimension %" PetscInt_FMT, dim);
230: }
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));
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));
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;
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));
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: }