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