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