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