xref: /petsc/src/dm/impls/plex/plexadapt.c (revision 57c00dc39306b6bef142b767add45be987161876) 
1 #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 #if defined(PETSC_HAVE_PRAGMATIC)
3 #include <pragmatic/cpragmatic.h>
4 #endif
5 #if defined(PETSC_HAVE_MMG)
6 #include <mmg/libmmg.h>
7 #endif
8 #if defined(PETSC_HAVE_PARMMG)
9 #include <parmmg/libparmmg.h>
10 #endif
11 
12 static PetscErrorCode DMPlexLabelToVolumeConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscReal refRatio, PetscReal maxVolumes[])
13 {
14   PetscInt       dim, c;
15   PetscErrorCode ierr;
16 
17   PetscFunctionBegin;
18   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
19   refRatio = refRatio == PETSC_DEFAULT ? (PetscReal) ((PetscInt) 1 << dim) : refRatio;
20   for (c = cStart; c < cEnd; c++) {
21     PetscReal vol;
22     PetscInt  closureSize = 0, cl;
23     PetscInt *closure     = NULL;
24     PetscBool anyRefine   = PETSC_FALSE;
25     PetscBool anyCoarsen  = PETSC_FALSE;
26     PetscBool anyKeep     = PETSC_FALSE;
27 
28     ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr);
29     maxVolumes[c - cStart] = vol;
30     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
31     for (cl = 0; cl < closureSize*2; cl += 2) {
32       const PetscInt point = closure[cl];
33       PetscInt       refFlag;
34 
35       ierr = DMLabelGetValue(adaptLabel, point, &refFlag);CHKERRQ(ierr);
36       switch (refFlag) {
37       case DM_ADAPT_REFINE:
38         anyRefine  = PETSC_TRUE;break;
39       case DM_ADAPT_COARSEN:
40         anyCoarsen = PETSC_TRUE;break;
41       case DM_ADAPT_KEEP:
42         anyKeep    = PETSC_TRUE;break;
43       case DM_ADAPT_DETERMINE:
44         break;
45       default:
46         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMPlex does not support refinement flag %D\n", refFlag);
47       }
48       if (anyRefine) break;
49     }
50     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
51     if (anyRefine) {
52       maxVolumes[c - cStart] = vol / refRatio;
53     } else if (anyKeep) {
54       maxVolumes[c - cStart] = vol;
55     } else if (anyCoarsen) {
56       maxVolumes[c - cStart] = vol * refRatio;
57     }
58   }
59   PetscFunctionReturn(0);
60 }
61 
62 static PetscErrorCode DMPlexLabelToMetricConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscInt vStart, PetscInt vEnd, PetscReal refRatio, Vec *metricVec)
63 {
64   DM              udm, coordDM;
65   PetscSection    coordSection;
66   Vec             coordinates, mb, mx;
67   Mat             A;
68   PetscScalar    *metric, *eqns;
69   const PetscReal coarseRatio = refRatio == PETSC_DEFAULT ? PetscSqr(0.5) : 1/refRatio;
70   PetscInt        dim, Nv, Neq, c, v;
71   PetscErrorCode  ierr;
72 
73   PetscFunctionBegin;
74   ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr);
75   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
76   ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr);
77   ierr = DMGetLocalSection(coordDM, &coordSection);CHKERRQ(ierr);
78   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
79   Nv   = vEnd - vStart;
80   ierr = VecCreateSeq(PETSC_COMM_SELF, Nv*PetscSqr(dim), metricVec);CHKERRQ(ierr);
81   ierr = VecGetArray(*metricVec, &metric);CHKERRQ(ierr);
82   Neq  = (dim*(dim+1))/2;
83   ierr = PetscMalloc1(PetscSqr(Neq), &eqns);CHKERRQ(ierr);
84   ierr = MatCreateSeqDense(PETSC_COMM_SELF, Neq, Neq, eqns, &A);CHKERRQ(ierr);
85   ierr = MatCreateVecs(A, &mx, &mb);CHKERRQ(ierr);
86   ierr = VecSet(mb, 1.0);CHKERRQ(ierr);
87   for (c = cStart; c < cEnd; ++c) {
88     const PetscScalar *sol;
89     PetscScalar       *cellCoords = NULL;
90     PetscReal          e[3], vol;
91     const PetscInt    *cone;
92     PetscInt           coneSize, cl, i, j, d, r;
93 
94     ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr);
95     /* Only works for simplices */
96     for (i = 0, r = 0; i < dim+1; ++i) {
97       for (j = 0; j < i; ++j, ++r) {
98         for (d = 0; d < dim; ++d) e[d] = PetscRealPart(cellCoords[i*dim+d] - cellCoords[j*dim+d]);
99         /* FORTRAN ORDERING */
100         switch (dim) {
101         case 2:
102           eqns[0*Neq+r] = PetscSqr(e[0]);
103           eqns[1*Neq+r] = 2.0*e[0]*e[1];
104           eqns[2*Neq+r] = PetscSqr(e[1]);
105           break;
106         case 3:
107           eqns[0*Neq+r] = PetscSqr(e[0]);
108           eqns[1*Neq+r] = 2.0*e[0]*e[1];
109           eqns[2*Neq+r] = 2.0*e[0]*e[2];
110           eqns[3*Neq+r] = PetscSqr(e[1]);
111           eqns[4*Neq+r] = 2.0*e[1]*e[2];
112           eqns[5*Neq+r] = PetscSqr(e[2]);
113           break;
114         }
115       }
116     }
117     ierr = MatSetUnfactored(A);CHKERRQ(ierr);
118     ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr);
119     ierr = MatLUFactor(A, NULL, NULL, NULL);CHKERRQ(ierr);
120     ierr = MatSolve(A, mb, mx);CHKERRQ(ierr);
121     ierr = VecGetArrayRead(mx, &sol);CHKERRQ(ierr);
122     ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr);
123     ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr);
124     ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr);
125     for (cl = 0; cl < coneSize; ++cl) {
126       const PetscInt v = cone[cl] - vStart;
127 
128       if (dim == 2) {
129         metric[v*4+0] += vol*coarseRatio*sol[0];
130         metric[v*4+1] += vol*coarseRatio*sol[1];
131         metric[v*4+2] += vol*coarseRatio*sol[1];
132         metric[v*4+3] += vol*coarseRatio*sol[2];
133       } else {
134         metric[v*9+0] += vol*coarseRatio*sol[0];
135         metric[v*9+1] += vol*coarseRatio*sol[1];
136         metric[v*9+3] += vol*coarseRatio*sol[1];
137         metric[v*9+2] += vol*coarseRatio*sol[2];
138         metric[v*9+6] += vol*coarseRatio*sol[2];
139         metric[v*9+4] += vol*coarseRatio*sol[3];
140         metric[v*9+5] += vol*coarseRatio*sol[4];
141         metric[v*9+7] += vol*coarseRatio*sol[4];
142         metric[v*9+8] += vol*coarseRatio*sol[5];
143       }
144     }
145     ierr = VecRestoreArrayRead(mx, &sol);CHKERRQ(ierr);
146   }
147   for (v = 0; v < Nv; ++v) {
148     const PetscInt *support;
149     PetscInt        supportSize, s;
150     PetscReal       vol, totVol = 0.0;
151 
152     ierr = DMPlexGetSupport(udm, v+vStart, &support);CHKERRQ(ierr);
153     ierr = DMPlexGetSupportSize(udm, v+vStart, &supportSize);CHKERRQ(ierr);
154     for (s = 0; s < supportSize; ++s) {ierr = DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL);CHKERRQ(ierr); totVol += vol;}
155     for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol;
156   }
157   ierr = PetscFree(eqns);CHKERRQ(ierr);
158   ierr = VecRestoreArray(*metricVec, &metric);CHKERRQ(ierr);
159   ierr = VecDestroy(&mx);CHKERRQ(ierr);
160   ierr = VecDestroy(&mb);CHKERRQ(ierr);
161   ierr = MatDestroy(&A);CHKERRQ(ierr);
162   ierr = DMDestroy(&udm);CHKERRQ(ierr);
163   PetscFunctionReturn(0);
164 }
165 
166 /*
167    Contains the list of registered DMPlexGenerators routines
168 */
169 extern PlexGeneratorFunctionList DMPlexGenerateList;
170 
171 PetscErrorCode DMPlexRefine_Internal(DM dm, DMLabel adaptLabel, DM *dmRefined)
172 {
173   PlexGeneratorFunctionList fl;
174   PetscErrorCode          (*refine)(DM,PetscReal*,DM*);
175   PetscErrorCode          (*adapt)(DM,DMLabel,DM*);
176   PetscErrorCode          (*refinementFunc)(const PetscReal [], PetscReal *);
177   char                      genname[PETSC_MAX_PATH_LEN], *name = NULL;
178   PetscReal                 refinementLimit;
179   PetscReal                *maxVolumes;
180   PetscInt                  dim, cStart, cEnd, c;
181   PetscBool                 flg, flg2, localized;
182   PetscErrorCode            ierr;
183 
184   PetscFunctionBegin;
185   ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
186   ierr = DMPlexGetRefinementLimit(dm, &refinementLimit);CHKERRQ(ierr);
187   ierr = DMPlexGetRefinementFunction(dm, &refinementFunc);CHKERRQ(ierr);
188   if (refinementLimit == 0.0 && !refinementFunc && !adaptLabel) PetscFunctionReturn(0);
189   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
190   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
191   ierr = PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_adaptor", genname, sizeof(genname), &flg);CHKERRQ(ierr);
192   if (flg) name = genname;
193   else {
194     ierr = PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, sizeof(genname), &flg2);CHKERRQ(ierr);
195     if (flg2) name = genname;
196   }
197 
198   fl = DMPlexGenerateList;
199   if (name) {
200     while (fl) {
201       ierr = PetscStrcmp(fl->name,name,&flg);CHKERRQ(ierr);
202       if (flg) {
203         refine = fl->refine;
204         adapt  = fl->adaptlabel;
205         goto gotit;
206       }
207       fl = fl->next;
208     }
209     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Grid refiner %s not registered",name);
210   } else {
211     while (fl) {
212       if (fl->dim < 0 || dim-1 == fl->dim) {
213         refine = fl->refine;
214         adapt  = fl->adaptlabel;
215         goto gotit;
216       }
217       fl = fl->next;
218     }
219     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No grid refiner of dimension %D registered",dim);
220   }
221 
222   gotit:
223   switch (dim) {
224     case 1:
225     case 2:
226     case 3:
227       if (adapt) {
228         ierr = (*adapt)(dm, adaptLabel, dmRefined);CHKERRQ(ierr);
229       } else {
230         ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
231         if (adaptLabel) {
232           ierr = DMPlexLabelToVolumeConstraint(dm, adaptLabel, cStart, cEnd, PETSC_DEFAULT, maxVolumes);CHKERRQ(ierr);
233         } else if (refinementFunc) {
234           for (c = cStart; c < cEnd; ++c) {
235             PetscReal vol, centroid[3];
236 
237             ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);CHKERRQ(ierr);
238             ierr = (*refinementFunc)(centroid, &maxVolumes[c-cStart]);CHKERRQ(ierr);
239           }
240         } else {
241           for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
242         }
243         ierr = (*refine)(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
244         ierr = PetscFree(maxVolumes);CHKERRQ(ierr);
245       }
246       break;
247     default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %D is not supported.", dim);
248   }
249   ((DM_Plex *) (*dmRefined)->data)->useHashLocation = ((DM_Plex *) dm->data)->useHashLocation;
250   ierr = DMCopyDisc(dm, *dmRefined);CHKERRQ(ierr);
251   if (localized) {ierr = DMLocalizeCoordinates(*dmRefined);CHKERRQ(ierr);}
252   PetscFunctionReturn(0);
253 }
254 
255 PetscErrorCode DMPlexCoarsen_Internal(DM dm, DMLabel adaptLabel, DM *dmCoarsened)
256 {
257   Vec            metricVec;
258   PetscInt       cStart, cEnd, vStart, vEnd;
259   DMLabel        bdLabel = NULL;
260   char           bdLabelName[PETSC_MAX_PATH_LEN];
261   PetscBool      localized, flg;
262   PetscErrorCode ierr;
263 
264   PetscFunctionBegin;
265   ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
266   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
267   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
268   ierr = DMPlexLabelToMetricConstraint(dm, adaptLabel, cStart, cEnd, vStart, vEnd, PETSC_DEFAULT, &metricVec);CHKERRQ(ierr);
269   ierr = PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_bd_label", bdLabelName, sizeof(bdLabelName), &flg);CHKERRQ(ierr);
270   if (flg) {ierr = DMGetLabel(dm, bdLabelName, &bdLabel);CHKERRQ(ierr);}
271   ierr = DMAdaptMetric_Plex(dm, metricVec, bdLabel, dmCoarsened);CHKERRQ(ierr);
272   ierr = VecDestroy(&metricVec);CHKERRQ(ierr);
273   ((DM_Plex *) (*dmCoarsened)->data)->useHashLocation = ((DM_Plex *) dm->data)->useHashLocation;
274   ierr = DMCopyDisc(dm, *dmCoarsened);CHKERRQ(ierr);
275   if (localized) {ierr = DMLocalizeCoordinates(*dmCoarsened);CHKERRQ(ierr);}
276   PetscFunctionReturn(0);
277 }
278 
279 PetscErrorCode DMAdaptLabel_Plex(DM dm, DMLabel adaptLabel, DM *dmAdapted)
280 {
281   IS              flagIS;
282   const PetscInt *flags;
283   PetscInt        defFlag, minFlag, maxFlag, numFlags, f;
284   PetscErrorCode  ierr;
285 
286   PetscFunctionBegin;
287   ierr = DMLabelGetDefaultValue(adaptLabel, &defFlag);CHKERRQ(ierr);
288   minFlag = defFlag;
289   maxFlag = defFlag;
290   ierr = DMLabelGetValueIS(adaptLabel, &flagIS);CHKERRQ(ierr);
291   ierr = ISGetLocalSize(flagIS, &numFlags);CHKERRQ(ierr);
292   ierr = ISGetIndices(flagIS, &flags);CHKERRQ(ierr);
293   for (f = 0; f < numFlags; ++f) {
294     const PetscInt flag = flags[f];
295 
296     minFlag = PetscMin(minFlag, flag);
297     maxFlag = PetscMax(maxFlag, flag);
298   }
299   ierr = ISRestoreIndices(flagIS, &flags);CHKERRQ(ierr);
300   ierr = ISDestroy(&flagIS);CHKERRQ(ierr);
301   {
302     PetscInt minMaxFlag[2], minMaxFlagGlobal[2];
303 
304     minMaxFlag[0] =  minFlag;
305     minMaxFlag[1] = -maxFlag;
306     ierr = MPI_Allreduce(minMaxFlag, minMaxFlagGlobal, 2, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRMPI(ierr);
307     minFlag =  minMaxFlagGlobal[0];
308     maxFlag = -minMaxFlagGlobal[1];
309   }
310   if (minFlag == maxFlag) {
311     switch (minFlag) {
312     case DM_ADAPT_DETERMINE:
313       *dmAdapted = NULL;break;
314     case DM_ADAPT_REFINE:
315       ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);
316       ierr = DMRefine(dm, MPI_COMM_NULL, dmAdapted);CHKERRQ(ierr);break;
317     case DM_ADAPT_COARSEN:
318       ierr = DMCoarsen(dm, MPI_COMM_NULL, dmAdapted);CHKERRQ(ierr);break;
319     default:
320       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP,"DMPlex does not support refinement flag %D\n", minFlag);
321     }
322   } else {
323     ierr = DMPlexSetRefinementUniform(dm, PETSC_FALSE);CHKERRQ(ierr);
324     ierr = DMPlexRefine_Internal(dm, adaptLabel, dmAdapted);CHKERRQ(ierr);
325   }
326   PetscFunctionReturn(0);
327 }
328 
329 PETSC_EXTERN PetscErrorCode DMAdaptMetric_Pragmatic_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew)
330 {
331 #if defined(PETSC_HAVE_PRAGMATIC)
332   MPI_Comm           comm;
333   const char        *bdName = "_boundary_";
334 #if 0
335   DM                 odm = dm;
336 #endif
337   DM                 udm, cdm;
338   DMLabel            bdLabelFull;
339   const char        *bdLabelName;
340   IS                 bdIS, globalVertexNum;
341   PetscSection       coordSection;
342   Vec                coordinates;
343   const PetscScalar *coords, *met;
344   const PetscInt    *bdFacesFull, *gV;
345   PetscInt          *bdFaces, *bdFaceIds, *l2gv;
346   PetscReal         *x, *y, *z, *metric;
347   PetscInt          *cells;
348   PetscInt           dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v;
349   PetscInt           off, maxConeSize, numBdFaces, f, bdSize;
350   PetscBool          flg;
351   DMLabel            bdLabelNew;
352   PetscReal         *coordsNew;
353   PetscInt          *bdTags;
354   PetscReal         *xNew[3] = {NULL, NULL, NULL};
355   PetscInt          *cellsNew;
356   PetscInt           d, numCellsNew, numVerticesNew;
357   PetscInt           numCornersNew, fStart, fEnd;
358   PetscMPIInt        numProcs;
359   PetscErrorCode     ierr;
360 
361   PetscFunctionBegin;
362 
363   /* Check for FEM adjacency flags */
364   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
365   ierr = MPI_Comm_size(comm, &numProcs);CHKERRMPI(ierr);
366   if (bdLabel) {
367     ierr = PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);CHKERRQ(ierr);
368     ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr);
369     if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
370   }
371 #if 0
372   /* Check for overlap by looking for cell in the SF */
373   if (!overlapped) {
374     ierr = DMPlexDistributeOverlap(odm, 1, NULL, &dm);CHKERRQ(ierr);
375     if (!dm) {dm = odm; ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);}
376   }
377 #endif
378 
379   /* Get mesh information */
380   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
381   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
382   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
383   ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr);
384   ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr);
385   numCells = cEnd - cStart;
386   if (numCells == 0) {
387     PetscMPIInt rank;
388 
389     ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
390     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform mesh adaptation because process %d does not own any cells.", rank);
391   }
392   numVertices = vEnd - vStart;
393   ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr);
394 
395   /* Get cell offsets */
396   for (c = 0, coff = 0; c < numCells; ++c) {
397     const PetscInt *cone;
398     PetscInt        coneSize, cl;
399 
400     ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr);
401     ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr);
402     for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
403   }
404 
405   /* Get local-to-global vertex map */
406   ierr = PetscCalloc1(numVertices, &l2gv);CHKERRQ(ierr);
407   ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr);
408   ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr);
409   for (v = 0, numLocVertices = 0; v < numVertices; ++v) {
410     if (gV[v] >= 0) ++numLocVertices;
411     l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v];
412   }
413   ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr);
414   ierr = DMDestroy(&udm);CHKERRQ(ierr);
415 
416   /* Get vertex coordinate arrays */
417   ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
418   ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr);
419   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
420   ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr);
421   for (v = vStart; v < vEnd; ++v) {
422     ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
423     x[v-vStart] = PetscRealPart(coords[off+0]);
424     if (dim > 1) y[v-vStart] = PetscRealPart(coords[off+1]);
425     if (dim > 2) z[v-vStart] = PetscRealPart(coords[off+2]);
426   }
427   ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr);
428 
429   /* Get boundary mesh */
430   ierr = DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);CHKERRQ(ierr);
431   ierr = DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);CHKERRQ(ierr);
432   ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr);
433   ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr);
434   ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr);
435   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
436     PetscInt *closure = NULL;
437     PetscInt  closureSize, cl;
438 
439     ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
440     for (cl = 0; cl < closureSize*2; cl += 2) {
441       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
442     }
443     ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
444   }
445   ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr);
446   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
447     PetscInt *closure = NULL;
448     PetscInt  closureSize, cl;
449 
450     ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
451     for (cl = 0; cl < closureSize*2; cl += 2) {
452       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart;
453     }
454     ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
455     if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);}
456     else         {bdFaceIds[f] = 1;}
457   }
458   ierr = ISDestroy(&bdIS);CHKERRQ(ierr);
459   ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr);
460 
461   /* Get metric */
462   ierr = VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");CHKERRQ(ierr);
463   ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr);
464   for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]);
465   ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr);
466 
467 #if 0
468   /* Destroy overlap mesh */
469   ierr = DMDestroy(&dm);CHKERRQ(ierr);
470 #endif
471   /* Send to Pragmatic and remesh */
472   switch (dim) {
473   case 2:
474     pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);
475     break;
476   case 3:
477     pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);
478     break;
479   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim);
480   }
481   pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
482   pragmatic_set_metric(metric);
483   pragmatic_adapt(((DM_Plex *) dm->data)->remeshBd ? 1 : 0);
484   ierr = PetscFree(l2gv);CHKERRQ(ierr);
485 
486   /* Retrieve mesh from Pragmatic and create new plex */
487   pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew);
488   ierr = PetscMalloc1(numVerticesNew*dim, &coordsNew);CHKERRQ(ierr);
489   switch (dim) {
490   case 2:
491     numCornersNew = 3;
492     ierr = PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);CHKERRQ(ierr);
493     pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]);
494     break;
495   case 3:
496     numCornersNew = 4;
497     ierr = PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);CHKERRQ(ierr);
498     pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]);
499     break;
500   default:
501     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim);
502   }
503   for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];}
504   ierr = PetscMalloc1(numCellsNew*(dim+1), &cellsNew);CHKERRQ(ierr);
505   pragmatic_get_elements(cellsNew);
506   ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew);CHKERRQ(ierr);
507 
508   /* Rebuild boundary label */
509   pragmatic_get_boundaryTags(&bdTags);
510   ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr);
511   ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr);
512   ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr);
513   ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr);
514   ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr);
515   for (c = cStart; c < cEnd; ++c) {
516 
517     /* Only for simplicial meshes */
518     coff = (c-cStart)*(dim+1);
519 
520     /* d is the local cell number of the vertex opposite to the face we are marking */
521     for (d = 0; d < dim+1; ++d) {
522       if (bdTags[coff+d]) {
523         const PetscInt  perm[4][4] = {{-1, -1, -1, -1}, {-1, -1, -1, -1}, {1, 2, 0, -1}, {3, 2, 1, 0}}; /* perm[d] = face opposite */
524         const PetscInt *cone;
525 
526         /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
527         ierr = DMPlexGetCone(*dmNew, c, &cone);CHKERRQ(ierr);
528         ierr = DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);CHKERRQ(ierr);
529       }
530     }
531   }
532 
533   /* Clean up */
534   switch (dim) {
535   case 2: ierr = PetscFree2(xNew[0], xNew[1]);CHKERRQ(ierr);
536   break;
537   case 3: ierr = PetscFree3(xNew[0], xNew[1], xNew[2]);CHKERRQ(ierr);
538   break;
539   }
540   ierr = PetscFree(cellsNew);CHKERRQ(ierr);
541   ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr);
542   ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr);
543   ierr = PetscFree(coordsNew);CHKERRQ(ierr);
544   pragmatic_finalize();
545   PetscFunctionReturn(0);
546 #else
547   PetscFunctionBegin;
548   SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic.");
549 #endif
550 }
551 
552 PETSC_EXTERN PetscErrorCode DMAdaptMetric_Mmg_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew)
553 {
554 #if defined(PETSC_HAVE_MMG)
555   MPI_Comm           comm;
556   const char        *bdName = "_boundary_";
557   DM                 udm, cdm;
558   DMLabel            bdLabelFull;
559   const char        *bdLabelName;
560   IS                 bdIS;
561   PetscSection       coordSection;
562   Vec                coordinates;
563   const PetscScalar *coords, *met;
564   const PetscInt    *bdFacesFull;
565   PetscInt          *bdFaces, *bdFaceIds;
566   PetscReal         *vertices, *metric, *verticesNew;
567   PetscInt          *cells, *cellsNew;
568   PetscInt          *facesNew, *faceTagsNew;
569   PetscInt          *verTags, *cellTags, *verTagsNew, *cellTagsNew;
570   PetscInt          *corners, *requiredCells, *requiredVer, *ridges, *requiredFaces;
571   PetscInt           dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, v;
572   PetscInt           off, maxConeSize, numBdFaces, f, bdSize, fStart, fEnd, i, j, k, Neq;
573   PetscInt           numCellsNew, numVerticesNew, numCornersNew, numFacesNew;
574   PetscBool          flg;
575   DMLabel            bdLabelNew;
576   MMG5_pMesh         mmg_mesh = NULL;
577   MMG5_pSol          mmg_metric = NULL;
578   PetscErrorCode     ierr;
579 
580   PetscFunctionBegin;
581   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
582   if (bdLabel) {
583     ierr = PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);CHKERRQ(ierr);
584     ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr);
585     if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
586   }
587 
588   /* Get mesh information */
589   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
590   Neq  = (dim*(dim+1))/2;
591   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
592   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
593   ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr);
594   ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr);
595   numCells    = cEnd - cStart;
596   numVertices = vEnd - vStart;
597 
598   /* Get cell offsets */
599   ierr = PetscMalloc1(numCells*maxConeSize, &cells);CHKERRQ(ierr);
600   for (c = 0, coff = 0; c < numCells; ++c) {
601     const PetscInt *cone;
602     PetscInt        coneSize, cl;
603 
604     ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr);
605     ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr);
606     for (cl = 0; cl < coneSize; ++cl) { cells[coff++] = cone[cl] - vStart + 1; }
607   }
608 
609   /* Get vertex coordinate array */
610   ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
611   ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr);
612   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
613   ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr);
614   ierr = PetscMalloc2(numVertices*Neq, &metric, dim*numVertices, &vertices);CHKERRQ(ierr);
615   for (v = 0; v < vEnd-vStart; ++v) {
616     ierr = PetscSectionGetOffset(coordSection, v+vStart, &off);CHKERRQ(ierr);
617     for (i = 0; i < dim; ++i) { vertices[dim*v+i] = PetscRealPart(coords[off+i]); }
618   }
619   ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr);
620   ierr = DMDestroy(&udm);CHKERRQ(ierr);
621 
622   /* Get boundary mesh */
623   ierr = DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);CHKERRQ(ierr);
624   ierr = DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);CHKERRQ(ierr);
625   ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr);
626   ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr);
627   ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr);
628   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
629     PetscInt *closure = NULL;
630     PetscInt  closureSize, cl;
631 
632     ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
633     for (cl = 0; cl < closureSize*2; cl += 2) {
634       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) { ++bdSize; }
635     }
636     ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
637   }
638   ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr);
639   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
640     PetscInt *closure = NULL;
641     PetscInt  closureSize, cl;
642 
643     ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
644     for (cl = 0; cl < closureSize*2; cl += 2) {
645       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) { bdFaces[bdSize++] = closure[cl] - vStart + 1; }
646     }
647     ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
648     if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);}
649     else         {bdFaceIds[f] = 1;}
650   }
651   ierr = ISDestroy(&bdIS);CHKERRQ(ierr);
652   ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr);
653 
654   /* Get metric */
655   ierr = VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");CHKERRQ(ierr);
656   ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr);
657   for (v = 0; v < (vEnd-vStart); ++v) {
658     for (i = 0, k = 0; i < dim; ++i) {
659       for (j = i; j < dim; ++j) {
660         metric[Neq*v+k] = PetscRealPart(met[dim*dim*v+dim*i+j]);
661         k++;
662       }
663     }
664   }
665   ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr);
666 
667   /* Send mesh to Mmg and remesh */
668   ierr = DMPlexMetricGetVerbosity(dm, &verbosity);CHKERRQ(ierr);
669   ierr = PetscCalloc2(numVertices, &verTags, numCells, &cellTags);CHKERRQ(ierr);
670   switch (dim) {
671   case 2:
672     ierr = MMG2D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end);
673     ierr = MMG2D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_verbose, 10);
674     ierr = MMG2D_Set_meshSize(mmg_mesh, numVertices, numCells, 0, numBdFaces);
675     ierr = MMG2D_Set_vertices(mmg_mesh, vertices, verTags);
676     ierr = MMG2D_Set_triangles(mmg_mesh, cells, cellTags);
677     ierr = MMG2D_Set_edges(mmg_mesh, bdFaces, bdFaceIds);
678     ierr = MMG2D_Set_solSize(mmg_mesh, mmg_metric, MMG5_Vertex, numVertices, MMG5_Tensor);
679     ierr = MMG2D_Set_tensorSols(mmg_metric, metric);
680     ierr = MMG2D_mmg2dlib(mmg_mesh, mmg_metric);
681     break;
682   case 3:
683     ierr = MMG3D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end);
684     ierr = MMG3D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_verbose, ctx->verbosity);
685     ierr = MMG3D_Set_meshSize(mmg_mesh, numVertices, numCells, 0, numBdFaces, 0, 0);
686     ierr = MMG3D_Set_vertices(mmg_mesh, vertices, verTags);
687     ierr = MMG3D_Set_tetrahedra(mmg_mesh, cells, cellTags);
688     ierr = MMG3D_Set_triangles(mmg_mesh, bdFaces, bdFaceIds);
689     ierr = MMG3D_Set_solSize(mmg_mesh, mmg_metric, MMG5_Vertex, numVertices, MMG5_Tensor);
690     ierr = MMG3D_Set_tensorSols(mmg_metric, metric);
691     ierr = MMG3D_mmg3dlib(mmg_mesh, mmg_metric);
692     break;
693   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Mmg adaptation defined for dimension %D", dim);
694   }
695 
696   /* Retrieve mesh from Mmg and create new Plex*/
697   switch (dim) {
698   case 2:
699     numCornersNew = 3;
700     ierr = MMG2D_Get_meshSize(mmg_mesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew);
701     ierr = PetscMalloc4(2*numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer);CHKERRQ(ierr);
702     ierr = PetscMalloc3(3*numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells);CHKERRQ(ierr);
703     ierr = PetscMalloc4(2*numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces);CHKERRQ(ierr);
704     ierr = MMG2D_Get_vertices(mmg_mesh, verticesNew, verTagsNew, corners, requiredVer);
705     ierr = MMG2D_Get_triangles(mmg_mesh, cellsNew, cellTagsNew, requiredCells);
706     ierr = MMG2D_Get_edges(mmg_mesh, facesNew, faceTagsNew, ridges, requiredFaces);
707     break;
708   case 3:
709     numCornersNew = 4;
710     ierr = MMG3D_Get_meshSize(mmg_mesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0);
711     ierr = PetscMalloc4(3*numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer);CHKERRQ(ierr);
712     ierr = PetscMalloc3(4*numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells);CHKERRQ(ierr);
713     ierr = PetscMalloc4(3*numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces);CHKERRQ(ierr);
714     ierr = MMG3D_Get_vertices(mmg_mesh, verticesNew, verTagsNew, corners, requiredVer);
715     ierr = MMG3D_Get_tetrahedra(mmg_mesh, cellsNew, cellTagsNew, requiredCells);
716     ierr = MMG3D_Get_triangles(mmg_mesh, facesNew, faceTagsNew, requiredFaces);
717     break;
718   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Mmg adaptation defined for dimension %D", dim);
719   }
720   for (i = 0; i < (dim+1)*numCellsNew; i++) { cellsNew[i] -= 1; }
721   for (i = 0; i < dim*numFacesNew; i++) { facesNew[i] -= 1; }
722   ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNew, NULL, NULL, dmNew);CHKERRQ(ierr);
723   switch (dim) {
724   case 2:
725     ierr = MMG2D_Free_all(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end);
726     break;
727   case 3:
728     ierr = MMG3D_Free_all(MMG5_ARG_start,MMG5_ARG_ppMesh,&mmg_mesh,MMG5_ARG_ppMet,&mmg_metric,MMG5_ARG_end);
729     break;
730   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No MMG adaptation defined for dimension %D", dim);
731   }
732 
733   /* Rebuild boundary labels */
734   ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr);
735   ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr);
736   ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr);
737   ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr);
738   for (i = 0; i < numFacesNew; i++) {
739     PetscInt        numCoveredPoints, numFaces = 0, facePoints[3];
740     const PetscInt *coveredPoints = NULL;
741 
742     for (j = 0; j < dim; ++j) { facePoints[j] = facesNew[i*dim+j]+vStart; }
743     ierr = DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr);
744     for (j = 0; j < numCoveredPoints; ++j) {
745       if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) {
746         numFaces++;
747         f = j;
748       }
749     }
750     if (numFaces != 1) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "%d vertices cannot define more than 1 facet (%d)", dim, numFaces);
751     ierr = DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i]);CHKERRQ(ierr);
752     ierr = DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr);
753   }
754 
755   /* Clean up */
756   ierr = PetscFree(cells);CHKERRQ(ierr);
757   ierr = PetscFree2(metric, vertices);CHKERRQ(ierr);
758   ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr);
759   ierr = PetscFree2(verTags, cellTags);CHKERRQ(ierr);
760   ierr = PetscFree4(verticesNew, verTagsNew, corners, requiredVer);CHKERRQ(ierr);
761   ierr = PetscFree3(cellsNew, cellTagsNew, requiredCells);CHKERRQ(ierr);
762   ierr = PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces);CHKERRQ(ierr);
763 
764   PetscFunctionReturn(0);
765 #else
766   PetscFunctionBegin;
767   SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-mmg.");
768   PetscFunctionReturn(0);
769 #endif
770 }
771 
772 PETSC_EXTERN PetscErrorCode DMAdaptMetric_ParMmg_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew)
773 {
774 #if defined(PETSC_HAVE_PARMMG)
775   MPI_Comm           comm, tmpComm;
776   const char        *bdName = "_boundary_";
777   DM                 udm, cdm;
778   DMLabel            bdLabelFull;
779   const char        *bdLabelName;
780   IS                 bdIS, globalVertexNum;
781   PetscSection       coordSection;
782   Vec                coordinates;
783   PetscSF            sf;
784   const PetscScalar *coords, *met;
785   const PetscInt    *bdFacesFull;
786   PetscInt          *bdFaces, *bdFaceIds;
787   PetscReal         *vertices, *metric, *verticesNew;
788   PetscInt          *cells;
789   const PetscInt    *gV, *ioffset, *irootloc, *roffset, *rmine, *rremote;
790   PetscInt          *numVerInterfaces, *ngbRanks, *verNgbRank, *interfaces_lv, *interfaces_gv, *intOffset;
791   PetscInt           niranks, nrranks, numNgbRanks, numVerNgbRanksTotal, count, sliceSize;
792   PetscInt          *verTags, *cellTags;
793   PetscInt           dim, Neq, cStart, cEnd, numCells, coff, vStart, vEnd, numVertices;
794   PetscInt           maxConeSize, numBdFaces, bdSize, fStart, fEnd;
795   PetscBool          flg;
796   DMLabel            bdLabelNew;
797   PetscReal         *verticesNewLoc;
798   PetscInt          *verTagsNew, *cellsNew, *cellTagsNew, *corners;
799   PetscInt          *requiredCells, *requiredVer, *facesNew, *faceTagsNew, *ridges, *requiredFaces;
800   PetscInt          *gv_new, *owners, *verticesNewSorted;
801   PetscInt           numCellsNew, numVerticesNew, numCornersNew, numFacesNew, numVerticesNewLoc;
802   PetscInt           i, j, k, p, r, n, v, c, f, off, lv, gv;
803   const PetscMPIInt *iranks, *rranks;
804   PetscMPIInt        numProcs, rank;
805   PMMG_pParMesh      parmesh = NULL;
806   PetscErrorCode     ierr;
807 
808   PetscFunctionBegin;
809   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
810   ierr = MPI_Comm_size(comm, &numProcs);CHKERRMPI(ierr);
811   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
812   MPI_Comm_dup(comm, &tmpComm);
813   if (bdLabel) {
814     ierr = PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);CHKERRQ(ierr);
815     ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr);
816     if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
817   }
818 
819   /* Get mesh information */
820   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
821   if (dim != 3) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "ParMmg only works in 3D.\n");
822   Neq  = (dim*(dim+1))/2;
823   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
824   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
825   ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr);
826   ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr);
827   numCells    = cEnd - cStart;
828   numVertices = vEnd - vStart;
829 
830   /* Get cell offsets */
831   ierr = PetscMalloc1(numCells*maxConeSize, &cells);CHKERRQ(ierr);
832   for (c = 0, coff = 0; c < numCells; ++c) {
833     const PetscInt *cone;
834     PetscInt        coneSize, cl;
835 
836     ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr);
837     ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr);
838     for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart+1;
839   }
840 
841   /* Get vertex coordinate array */
842   ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
843   ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr);
844   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
845   ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr);
846   ierr = PetscMalloc2(numVertices*Neq, &metric, dim*numVertices, &vertices);CHKERRQ(ierr);
847   for (v = 0; v < vEnd-vStart; ++v) {
848     ierr = PetscSectionGetOffset(coordSection, v+vStart, &off);CHKERRQ(ierr);
849     for (i = 0; i < dim; ++i) { vertices[dim*v+i]=PetscRealPart(coords[off+i]); }
850   }
851   ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr);
852 
853   /* Get boundary mesh */
854   ierr = DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);CHKERRQ(ierr);
855   ierr = DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);CHKERRQ(ierr);
856   ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr);
857   ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr);
858   ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr);
859   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
860     PetscInt *closure = NULL;
861     PetscInt  closureSize, cl;
862 
863     ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
864     for (cl = 0; cl < closureSize*2; cl += 2) {
865       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
866     }
867     ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
868   }
869   ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr);
870   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
871     PetscInt *closure = NULL;
872     PetscInt  closureSize, cl;
873 
874     ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
875     for (cl = 0; cl < closureSize*2; cl += 2) {
876       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart+1;
877     }
878     ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
879     if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);}
880     else         {bdFaceIds[f] = 1;}
881   }
882   ierr = ISDestroy(&bdIS);CHKERRQ(ierr);
883   ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr);
884 
885   /* Get metric */
886   ierr = VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");CHKERRQ(ierr);
887   ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr);
888   for (v = 0; v < (vEnd-vStart); ++v) {
889     for (i = 0, k = 0; i < dim; ++i) {
890       for (j = i; j < dim; ++j) {
891         metric[Neq*v+k] = PetscRealPart(met[dim*dim*v+dim*i+j]);
892         k++;
893       }
894     }
895   }
896   ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr);
897 
898   /* Build ParMMG communicators: the list of vertices between two partitions  */
899   niranks = nrranks = 0;
900   numNgbRanks = 0;
901   if (numProcs > 1) {
902     ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
903     ierr = PetscSFSetUp(sf);CHKERRQ(ierr);
904     ierr = PetscSFGetLeafRanks(sf, &niranks, &iranks, &ioffset, &irootloc);CHKERRQ(ierr);
905     ierr = PetscSFGetRootRanks(sf, &nrranks, &rranks, &roffset, &rmine, &rremote);CHKERRQ(ierr);
906     ierr = PetscCalloc1(numProcs, &numVerInterfaces);CHKERRQ(ierr);
907 
908     /* Counting */
909     for (r = 0; r < niranks; ++r) {
910       for (i=ioffset[r], count=0; i<ioffset[r+1]; ++i) {
911         if (irootloc[i] >= vStart && irootloc[i] < vEnd) {count++; }
912       }
913       numVerInterfaces[iranks[r]] += count;
914     }
915     for (r = 0; r < nrranks; ++r) {
916       for (i=roffset[r], count=0; i<roffset[r+1]; ++i) {
917         if (rmine[i] >= vStart && rmine[i] < vEnd) {count++; }
918       }
919       numVerInterfaces[rranks[r]] += count;
920     }
921     for (p = 0; p < numProcs; ++p) {
922       if (numVerInterfaces[p]) { numNgbRanks++; }
923     }
924     ierr = PetscMalloc2(numNgbRanks, &ngbRanks, numNgbRanks, &verNgbRank);CHKERRQ(ierr);
925     for (p = 0, n = 0; p < numProcs; ++p) {
926       if (numVerInterfaces[p]) {
927         ngbRanks[n] = p;
928         verNgbRank[n] = numVerInterfaces[p];
929         n++;
930       }
931     }
932     numVerNgbRanksTotal = 0;
933     for (p = 0; p < numNgbRanks; ++p) { numVerNgbRanksTotal += verNgbRank[p]; }
934 
935     /* For each neighbor, fill in interface arrays */
936     ierr = PetscMalloc3(numVerNgbRanksTotal, &interfaces_lv, numVerNgbRanksTotal, &interfaces_gv,  numNgbRanks+1, &intOffset);CHKERRQ(ierr);
937     intOffset[0] = 0;
938     for (p = 0, r = 0, i = 0; p < numNgbRanks; ++p) {
939       intOffset[p+1] = intOffset[p];
940       if (iranks && iranks[i] == ngbRanks[p]) {
941 
942         /* Add the right slice of irootloc at the right place */
943         sliceSize = ioffset[i+1]-ioffset[i];
944         for (j = 0, count = 0; j < sliceSize; ++j) {
945           if (ioffset[i]+j >= ioffset[niranks]) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Offset out of range");
946           v = irootloc[ioffset[i]+j];
947           if (v >= vStart && v < vEnd) {
948             if (intOffset[p+1]+count >= numVerNgbRanksTotal) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Offset out of range");
949             interfaces_lv[intOffset[p+1]+count] = v-vStart;
950             count++;
951           }
952         }
953         intOffset[p+1] += count;
954         i++;
955       }
956       if (rranks && rranks[r] == ngbRanks[p]) {
957 
958         /* Add the right slice of rmine at the right place */
959         sliceSize = roffset[r+1]-roffset[r];
960         for (j = 0, count = 0; j < sliceSize; ++j) {
961           if (roffset[r]+j >= roffset[nrranks]) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Offset out of range");
962           v = rmine[roffset[r]+j];
963           if (v >= vStart && v < vEnd) {
964             if (intOffset[p+1]+count >= numVerNgbRanksTotal) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Offset out of range");
965             interfaces_lv[intOffset[p+1]+count] = v-vStart;
966             count++;
967           }
968         }
969         intOffset[p+1] += count;
970         r++;
971       }
972       if (intOffset[p+1] != intOffset[p] + verNgbRank[p]) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unequal offsets");
973     }
974     ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr);
975     ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr);
976     for (i = 0; i < numVerNgbRanksTotal; ++i) {
977       v = gV[interfaces_lv[i]];
978       interfaces_gv[i] = v < 0 ? -v-1 : v;
979       interfaces_lv[i] += 1;
980       interfaces_gv[i] += 1;
981     }
982     ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr);
983     ierr = PetscFree(numVerInterfaces);CHKERRQ(ierr);
984   }
985   ierr = DMDestroy(&udm);CHKERRQ(ierr);
986 
987   /* Send the data to ParMmg and remesh */
988   ierr = PetscCalloc2(numVertices, &verTags, numCells, &cellTags);CHKERRQ(ierr);
989   ierr = PMMG_Init_parMesh(PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_pMesh, PMMG_ARG_pMet, PMMG_ARG_dim, 3, PMMG_ARG_MPIComm, tmpComm, PMMG_ARG_end);
990   ierr = PMMG_Set_meshSize(parmesh, numVertices, numCells, 0, numBdFaces, 0, 0);
991   ierr = PMMG_Set_iparameter(parmesh, PMMG_IPARAM_APImode, PMMG_APIDISTRIB_nodes);
992   ierr = PMMG_Set_iparameter(parmesh, PMMG_IPARAM_verbose, 10);
993   ierr = PMMG_Set_iparameter(parmesh, PMMG_IPARAM_globalNum, 1);
994   ierr = PMMG_Set_vertices(parmesh, vertices, verTags);
995   for (i=0; i<numCells; i++) ierr = PMMG_Set_tetrahedron(parmesh, cells[4*i+0], cells[4*i+1], cells[4*i+2], cells[4*i+3], 0, i+1);
996   ierr = PMMG_Set_triangles(parmesh, bdFaces, bdFaceIds);
997   ierr = PMMG_Set_metSize(parmesh, MMG5_Vertex, numVertices, MMG5_Tensor);
998   for (i=0; i<numVertices; i++) {
999     PMMG_Set_tensorMet(parmesh, metric[6*i], metric[6*i+1], metric[6*i+2], metric[6*i+3], metric[6*i+4], metric[6*i+5], i+1);
1000   }
1001   ierr = PMMG_Set_numberOfNodeCommunicators(parmesh, numNgbRanks);
1002   for (c = 0; c < numNgbRanks; ++c) {
1003     ierr = PMMG_Set_ithNodeCommunicatorSize(parmesh, c, ngbRanks[c], intOffset[c+1]-intOffset[c]);
1004     ierr = PMMG_Set_ithNodeCommunicator_nodes(parmesh, c, &interfaces_lv[intOffset[c]], &interfaces_gv[intOffset[c]], 1);
1005   }
1006   ierr = PMMG_Set_iparameter(parmesh, PMMG_IPARAM_niter, 3);
1007   ierr = PMMG_parmmglib_distributed(parmesh);
1008   ierr = PetscFree(cells);CHKERRQ(ierr);
1009   ierr = PetscFree2(metric, vertices);CHKERRQ(ierr);
1010   ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr);
1011   ierr = PetscFree2(verTags, cellTags);CHKERRQ(ierr);
1012   if (numProcs > 1) {
1013     ierr = PetscFree2(ngbRanks, verNgbRank);CHKERRQ(ierr);
1014     ierr = PetscFree3(interfaces_lv, interfaces_gv, intOffset);CHKERRQ(ierr);
1015   }
1016 
1017   /* Retrieve mesh from Mmg and create new Plex */
1018   numCornersNew = 4;
1019   ierr = PMMG_Get_meshSize(parmesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0);
1020   ierr = PetscMalloc4(dim*numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer);CHKERRQ(ierr);
1021   ierr = PetscMalloc3((dim+1)*numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells);CHKERRQ(ierr);
1022   ierr = PetscMalloc4(dim*numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces);CHKERRQ(ierr);
1023   ierr = PMMG_Get_vertices(parmesh, verticesNew, verTagsNew, corners, requiredVer);
1024   ierr = PMMG_Get_tetrahedra(parmesh, cellsNew, cellTagsNew, requiredCells);
1025   ierr = PMMG_Get_triangles(parmesh, facesNew, faceTagsNew, requiredFaces);
1026   ierr = PetscMalloc2(numVerticesNew, &owners, numVerticesNew, &gv_new);
1027   ierr = PMMG_Set_iparameter(parmesh, PMMG_IPARAM_globalNum, 1);
1028   ierr = PMMG_Get_verticesGloNum(parmesh, gv_new, owners);
1029   for (i = 0; i < dim*numFacesNew; ++i) { facesNew[i] -= 1; }
1030   for (i = 0; i < (dim+1)*numCellsNew; ++i) {
1031     cellsNew[i] = gv_new[cellsNew[i]-1]-1;
1032   }
1033   numVerticesNewLoc = 0;
1034   for (i = 0; i < numVerticesNew; ++i) {
1035     if (owners[i] == rank) { numVerticesNewLoc++; }
1036   }
1037   ierr = PetscMalloc2(numVerticesNewLoc*dim, &verticesNewLoc, numVerticesNew, &verticesNewSorted);CHKERRQ(ierr);
1038 
1039   for (i = 0, c = 0; i < numVerticesNew; i++) {
1040     if (owners[i] == rank) {
1041       for (j=0; j<dim; ++j) { verticesNewLoc[dim*c+j] = verticesNew[dim*i+j]; }
1042         c++;
1043     }
1044   }
1045   ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNewLoc, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNewLoc, NULL, &verticesNewSorted, dmNew);CHKERRQ(ierr);
1046   ierr = PMMG_Free_all(PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_end);
1047 
1048   /* Rebuild boundary label*/
1049   ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr);
1050   ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr);
1051   ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr);
1052   ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr);
1053   ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr);
1054   for (i = 0; i < numFacesNew; i++) {
1055     PetscInt        numCoveredPoints, numFaces = 0, facePoints[3];
1056     const PetscInt *coveredPoints = NULL;
1057 
1058     for (j = 0; j < dim; ++j) {
1059       lv = facesNew[i*dim+j];
1060       gv = gv_new[lv]-1;
1061       ierr = PetscFindInt(gv, numVerticesNew, verticesNewSorted, &lv);CHKERRQ(ierr);
1062       facePoints[j] = lv+vStart;
1063     }
1064     ierr = DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr);
1065     for (j = 0; j < numCoveredPoints; ++j) {
1066       if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) {
1067         numFaces++;
1068         f = j;
1069       }
1070     }
1071     if (numFaces != 1) { SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "%d vertices cannot define more than 1 facet (%d)", dim, numFaces); }
1072     ierr = DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i]);CHKERRQ(ierr);
1073     ierr = DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);CHKERRQ(ierr);
1074   }
1075 
1076   /* Clean up */
1077   ierr = PetscFree4(verticesNew, verTagsNew, corners, requiredVer);CHKERRQ(ierr);
1078   ierr = PetscFree3(cellsNew, cellTagsNew, requiredCells);CHKERRQ(ierr);
1079   ierr = PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces);CHKERRQ(ierr);
1080   ierr = PetscFree2(owners, gv_new);CHKERRQ(ierr);
1081   ierr = PetscFree2(verticesNewLoc, verticesNewSorted);CHKERRQ(ierr);
1082   PetscFunctionReturn(0);
1083 #else
1084   PetscFunctionBegin;
1085   SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-parmmg.");
1086   PetscFunctionReturn(0);
1087 #endif
1088 }
1089 
1090 /*
1091   DMAdaptMetric_Plex - Generates a new mesh conforming to a metric field.
1092 
1093   Input Parameters:
1094 + dm - The DM object
1095 . vertexMetric - The metric to which the mesh is adapted, defined vertex-wise in a LOCAL vector
1096 - bdLabel - Label for boundary tags which are preserved in dmNew, or NULL. Should not be named "_boundary_".
1097 
1098   Output Parameter:
1099 . dmNew  - the new DM
1100 
1101   Level: advanced
1102 
1103 .seealso: DMCoarsen(), DMRefine()
1104 */
1105 PetscErrorCode DMAdaptMetric_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew)
1106 {
1107   PetscInt remesher = 2;
1108 
1109   PetscFunctionBegin;
1110   switch (remesher) {
1111   case 0:
1112     DMAdaptMetricPragmatic_Plex(dm, vertexMetric, bdLabel, dmNew);
1113     break;
1114   case 1:
1115     DMAdaptMetricMMG_Plex(dm, vertexMetric, bdLabel, dmNew);
1116     break;
1117   case 2:
1118     DMAdaptMetricParMMG_Plex(dm, vertexMetric, bdLabel, dmNew);
1119     break;
1120   }
1121   PetscFunctionReturn(0);
1122 }
1123