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: }