xref: /petsc/src/dm/impls/plex/plexadapt.c (revision 19532ea364b07ecd049071e360c10aa6fe9ed42f) 
1 #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 #ifdef PETSC_HAVE_PRAGMATIC
3 #include <pragmatic/cpragmatic.h>
4 #endif
5 
6 static PetscErrorCode DMPlexLabelToVolumeConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscReal refRatio, PetscReal maxVolumes[])
7 {
8   PetscInt       dim, c;
9   PetscErrorCode ierr;
10 
11   PetscFunctionBegin;
12   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
13   refRatio = refRatio == PETSC_DEFAULT ? (PetscReal) ((PetscInt) 1 << dim) : refRatio;
14   for (c = cStart; c < cEnd; c++) {
15     PetscReal vol;
16     PetscInt  closureSize = 0, cl;
17     PetscInt *closure     = NULL;
18     PetscBool anyRefine   = PETSC_FALSE;
19     PetscBool anyCoarsen  = PETSC_FALSE;
20     PetscBool anyKeep     = PETSC_FALSE;
21 
22     ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr);
23     maxVolumes[c - cStart] = vol;
24     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
25     for (cl = 0; cl < closureSize*2; cl += 2) {
26       const PetscInt point = closure[cl];
27       PetscInt       refFlag;
28 
29       ierr = DMLabelGetValue(adaptLabel, point, &refFlag);CHKERRQ(ierr);
30       switch (refFlag) {
31       case DM_ADAPT_REFINE:
32         anyRefine  = PETSC_TRUE;break;
33       case DM_ADAPT_COARSEN:
34         anyCoarsen = PETSC_TRUE;break;
35       case DM_ADAPT_KEEP:
36         anyKeep    = PETSC_TRUE;break;
37       case DM_ADAPT_DETERMINE:
38         break;
39       default:
40         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMPlex does not support refinement flag %D\n", refFlag);break;
41       }
42       if (anyRefine) break;
43     }
44     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
45     if (anyRefine) {
46       maxVolumes[c - cStart] = vol / refRatio;
47     } else if (anyKeep) {
48       maxVolumes[c - cStart] = vol;
49     } else if (anyCoarsen) {
50       maxVolumes[c - cStart] = vol * refRatio;
51     }
52   }
53   PetscFunctionReturn(0);
54 }
55 
56 static PetscErrorCode DMPlexLabelToMetricConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscInt vStart, PetscInt vEnd, PetscReal refRatio, Vec *metricVec)
57 {
58   DM              udm, coordDM;
59   PetscSection    coordSection;
60   Vec             coordinates, mb, mx;
61   Mat             A;
62   PetscScalar    *metric, *eqns;
63   const PetscReal coarseRatio = refRatio == PETSC_DEFAULT ? PetscSqr(0.5) : 1/refRatio;
64   PetscInt        dim, Nv, Neq, c, v;
65   PetscErrorCode  ierr;
66 
67   PetscFunctionBegin;
68   ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr);
69   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
70   ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr);
71   ierr = DMGetDefaultSection(coordDM, &coordSection);CHKERRQ(ierr);
72   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
73   Nv   = vEnd - vStart;
74   ierr = VecCreateSeq(PETSC_COMM_SELF, Nv*PetscSqr(dim), metricVec);CHKERRQ(ierr);
75   ierr = VecGetArray(*metricVec, &metric);CHKERRQ(ierr);
76   Neq  = (dim*(dim+1))/2;
77   ierr = PetscMalloc1(PetscSqr(Neq), &eqns);CHKERRQ(ierr);
78   ierr = MatCreateSeqDense(PETSC_COMM_SELF, Neq, Neq, eqns, &A);CHKERRQ(ierr);
79   ierr = MatCreateVecs(A, &mx, &mb);CHKERRQ(ierr);
80   ierr = VecSet(mb, 1.0);CHKERRQ(ierr);
81   for (c = cStart; c < cEnd; ++c) {
82     const PetscScalar *sol;
83     PetscScalar       *cellCoords = NULL;
84     PetscReal          e[3], vol;
85     const PetscInt    *cone;
86     PetscInt           coneSize, cl, i, j, d, r;
87 
88     ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr);
89     /* Only works for simplices */
90     for (i = 0, r = 0; i < dim+1; ++i) {
91       for (j = 0; j < i; ++j, ++r) {
92         for (d = 0; d < dim; ++d) e[d] = PetscRealPart(cellCoords[i*dim+d] - cellCoords[j*dim+d]);
93         /* FORTRAN ORDERING */
94         switch (dim) {
95         case 2:
96           eqns[0*Neq+r] = PetscSqr(e[0]);
97           eqns[1*Neq+r] = 2.0*e[0]*e[1];
98           eqns[2*Neq+r] = PetscSqr(e[1]);
99           break;
100         case 3:
101           eqns[0*Neq+r] = PetscSqr(e[0]);
102           eqns[1*Neq+r] = 2.0*e[0]*e[1];
103           eqns[2*Neq+r] = 2.0*e[0]*e[2];
104           eqns[3*Neq+r] = PetscSqr(e[1]);
105           eqns[4*Neq+r] = 2.0*e[1]*e[2];
106           eqns[5*Neq+r] = PetscSqr(e[2]);
107           break;
108         }
109       }
110     }
111     ierr = MatSetUnfactored(A);CHKERRQ(ierr);
112     ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr);
113     ierr = MatLUFactor(A, NULL, NULL, NULL);CHKERRQ(ierr);
114     ierr = MatSolve(A, mb, mx);CHKERRQ(ierr);
115     ierr = VecGetArrayRead(mx, &sol);CHKERRQ(ierr);
116     ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr);
117     ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr);
118     ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr);
119     for (cl = 0; cl < coneSize; ++cl) {
120       const PetscInt v = cone[cl] - vStart;
121 
122       if (dim == 2) {
123         metric[v*4+0] += vol*coarseRatio*sol[0];
124         metric[v*4+1] += vol*coarseRatio*sol[1];
125         metric[v*4+2] += vol*coarseRatio*sol[1];
126         metric[v*4+3] += vol*coarseRatio*sol[2];
127       } else {
128         metric[v*9+0] += vol*coarseRatio*sol[0];
129         metric[v*9+1] += vol*coarseRatio*sol[1];
130         metric[v*9+3] += vol*coarseRatio*sol[1];
131         metric[v*9+2] += vol*coarseRatio*sol[2];
132         metric[v*9+6] += vol*coarseRatio*sol[2];
133         metric[v*9+4] += vol*coarseRatio*sol[3];
134         metric[v*9+5] += vol*coarseRatio*sol[4];
135         metric[v*9+7] += vol*coarseRatio*sol[4];
136         metric[v*9+8] += vol*coarseRatio*sol[5];
137       }
138     }
139     ierr = VecRestoreArrayRead(mx, &sol);CHKERRQ(ierr);
140   }
141   for (v = 0; v < Nv; ++v) {
142     const PetscInt *support;
143     PetscInt        supportSize, s;
144     PetscReal       vol, totVol = 0.0;
145 
146     ierr = DMPlexGetSupport(udm, v+vStart, &support);CHKERRQ(ierr);
147     ierr = DMPlexGetSupportSize(udm, v+vStart, &supportSize);CHKERRQ(ierr);
148     for (s = 0; s < supportSize; ++s) {ierr = DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL);CHKERRQ(ierr); totVol += vol;}
149     for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol;
150   }
151   ierr = PetscFree(eqns);CHKERRQ(ierr);
152   ierr = VecRestoreArray(*metricVec, &metric);CHKERRQ(ierr);
153   ierr = VecDestroy(&mx);CHKERRQ(ierr);
154   ierr = VecDestroy(&mb);CHKERRQ(ierr);
155   ierr = MatDestroy(&A);CHKERRQ(ierr);
156   ierr = DMDestroy(&udm);CHKERRQ(ierr);
157   PetscFunctionReturn(0);
158 }
159 
160 PetscErrorCode DMPlexRefine_Internal(DM dm, DMLabel adaptLabel, DM *dmRefined)
161 {
162   PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *);
163   PetscReal        refinementLimit;
164   PetscInt         dim, cStart, cEnd;
165   char             genname[1024], *name = NULL;
166   PetscBool        isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg, localized;
167   PetscErrorCode   ierr;
168 
169   PetscFunctionBegin;
170   ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
171   ierr = DMPlexGetRefinementLimit(dm, &refinementLimit);CHKERRQ(ierr);
172   ierr = DMPlexGetRefinementFunction(dm, &refinementFunc);CHKERRQ(ierr);
173   if (refinementLimit == 0.0 && !refinementFunc && !adaptLabel) PetscFunctionReturn(0);
174   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
175   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
176   ierr = PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);CHKERRQ(ierr);
177   if (flg) name = genname;
178   if (name) {
179     ierr = PetscStrcmp(name, "triangle", &isTriangle);CHKERRQ(ierr);
180     ierr = PetscStrcmp(name, "tetgen",   &isTetgen);CHKERRQ(ierr);
181     ierr = PetscStrcmp(name, "ctetgen",  &isCTetgen);CHKERRQ(ierr);
182   }
183   switch (dim) {
184   case 2:
185     if (!name || isTriangle) {
186 #if defined(PETSC_HAVE_TRIANGLE)
187       PetscReal *maxVolumes;
188       PetscInt  c;
189 
190       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
191       if (adaptLabel) {
192         ierr = DMPlexLabelToVolumeConstraint(dm, adaptLabel, cStart, cEnd, PETSC_DEFAULT, maxVolumes);CHKERRQ(ierr);
193       } else if (refinementFunc) {
194         for (c = cStart; c < cEnd; ++c) {
195           PetscReal vol, centroid[3];
196           PetscReal maxVol;
197 
198           ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);CHKERRQ(ierr);
199           ierr = (*refinementFunc)(centroid, &maxVol);CHKERRQ(ierr);
200           maxVolumes[c - cStart] = (double) maxVol;
201         }
202       } else {
203         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
204       }
205 #if !defined(PETSC_USE_REAL_DOUBLE)
206       {
207         double *mvols;
208         ierr = PetscMalloc1(cEnd - cStart,&mvols);CHKERRQ(ierr);
209         for (c = 0; c < cEnd-cStart; ++c) mvols[c] = (double)maxVolumes[c];
210         ierr = DMPlexRefine_Triangle(dm, mvols, dmRefined);CHKERRQ(ierr);
211         ierr = PetscFree(mvols);CHKERRQ(ierr);
212       }
213 #else
214       ierr = DMPlexRefine_Triangle(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
215 #endif
216       ierr = PetscFree(maxVolumes);CHKERRQ(ierr);
217 #else
218       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
219 #endif
220     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
221     break;
222   case 3:
223     if (!name || isCTetgen || isTetgen) {
224       PetscReal *maxVolumes;
225       PetscInt   c;
226 
227       ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr);
228       if (adaptLabel) {
229         ierr = DMPlexLabelToVolumeConstraint(dm, adaptLabel, cStart, cEnd, PETSC_DEFAULT, maxVolumes);CHKERRQ(ierr);
230       } else if (refinementFunc) {
231         for (c = cStart; c < cEnd; ++c) {
232           PetscReal vol, centroid[3];
233 
234           ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);CHKERRQ(ierr);
235           ierr = (*refinementFunc)(centroid, &maxVolumes[c-cStart]);CHKERRQ(ierr);
236         }
237       } else {
238         for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
239       }
240       if (!name) {
241 #if defined(PETSC_HAVE_CTETGEN)
242         ierr = DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
243 #elif defined(PETSC_HAVE_TETGEN)
244         ierr = DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
245 #else
246         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "External package CTetgen or Tetgen needed.\nPlease reconfigure with '--download-ctetgen' or '--with-clanguage=cxx --download-tetgen'.");
247 #endif
248       } else if (isCTetgen) {
249 #if defined(PETSC_HAVE_CTETGEN)
250         ierr = DMPlexRefine_CTetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
251 #else
252         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
253 #endif
254       } else {
255 #if defined(PETSC_HAVE_TETGEN)
256         ierr = DMPlexRefine_Tetgen(dm, maxVolumes, dmRefined);CHKERRQ(ierr);
257 #else
258         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --download-tetgen.");
259 #endif
260       }
261       ierr = PetscFree(maxVolumes);CHKERRQ(ierr);
262     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
263     break;
264   default:
265     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
266   }
267   ierr = DMCopyBoundary(dm, *dmRefined);CHKERRQ(ierr);
268   if (localized) {ierr = DMLocalizeCoordinates(*dmRefined);CHKERRQ(ierr);}
269   PetscFunctionReturn(0);
270 }
271 
272 PetscErrorCode DMPlexCoarsen_Internal(DM dm, DMLabel adaptLabel, DM *dmCoarsened)
273 {
274   Vec            metricVec;
275   PetscInt       cStart, cEnd, vStart, vEnd;
276   DMLabel        bdLabel = NULL;
277   char           bdLabelName[PETSC_MAX_PATH_LEN];
278   PetscBool      localized, flg;
279   PetscErrorCode ierr;
280 
281   PetscFunctionBegin;
282   ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
283   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
284   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
285   ierr = DMPlexLabelToMetricConstraint(dm, adaptLabel, cStart, cEnd, vStart, vEnd, PETSC_DEFAULT, &metricVec);CHKERRQ(ierr);
286   ierr = PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_bd_label", bdLabelName, PETSC_MAX_PATH_LEN-1, &flg);CHKERRQ(ierr);
287   if (flg) {ierr = DMGetLabel(dm, bdLabelName, &bdLabel);CHKERRQ(ierr);}
288   ierr = DMAdaptMetric_Plex(dm, metricVec, bdLabel, dmCoarsened);CHKERRQ(ierr);
289   ierr = VecDestroy(&metricVec);CHKERRQ(ierr);
290   if (localized) {ierr = DMLocalizeCoordinates(*dmCoarsened);CHKERRQ(ierr);}
291   PetscFunctionReturn(0);
292 }
293 
294 PetscErrorCode DMAdaptLabel_Plex(DM dm, DMLabel adaptLabel, DM *dmAdapted)
295 {
296   IS              flagIS;
297   const PetscInt *flags;
298   PetscInt        defFlag, minFlag, maxFlag, numFlags, f;
299   PetscErrorCode  ierr;
300 
301   PetscFunctionBegin;
302   ierr = DMLabelGetDefaultValue(adaptLabel, &defFlag);CHKERRQ(ierr);
303   minFlag = defFlag;
304   maxFlag = defFlag;
305   ierr = DMLabelGetValueIS(adaptLabel, &flagIS);CHKERRQ(ierr);
306   ierr = ISGetLocalSize(flagIS, &numFlags);CHKERRQ(ierr);
307   ierr = ISGetIndices(flagIS, &flags);CHKERRQ(ierr);
308   for (f = 0; f < numFlags; ++f) {
309     const PetscInt flag = flags[f];
310 
311     minFlag = PetscMin(minFlag, flag);
312     maxFlag = PetscMax(maxFlag, flag);
313   }
314   ierr = ISRestoreIndices(flagIS, &flags);CHKERRQ(ierr);
315   ierr = ISDestroy(&flagIS);CHKERRQ(ierr);
316   {
317     PetscInt minMaxFlag[2], minMaxFlagGlobal[2];
318 
319     minMaxFlag[0] =  minFlag;
320     minMaxFlag[1] = -maxFlag;
321     ierr = MPI_Allreduce(minMaxFlag, minMaxFlagGlobal, 2, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
322     minFlag =  minMaxFlagGlobal[0];
323     maxFlag = -minMaxFlagGlobal[1];
324   }
325   if (minFlag == maxFlag) {
326     switch (minFlag) {
327     case DM_ADAPT_DETERMINE:
328       *dmAdapted = NULL;break;
329     case DM_ADAPT_REFINE:
330       ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);
331       ierr = DMRefine(dm, MPI_COMM_NULL, dmAdapted);CHKERRQ(ierr);break;
332     case DM_ADAPT_COARSEN:
333       ierr = DMCoarsen(dm, MPI_COMM_NULL, dmAdapted);CHKERRQ(ierr);break;
334     default:
335       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP,"DMPlex does not support refinement flag %D\n", minFlag);break;
336     }
337   } else {
338     ierr = DMPlexSetRefinementUniform(dm, PETSC_FALSE);CHKERRQ(ierr);
339     ierr = DMPlexRefine_Internal(dm, adaptLabel, dmAdapted);CHKERRQ(ierr);
340   }
341   PetscFunctionReturn(0);
342 }
343 
344 /*
345   DMAdaptMetric_Plex - Generates a new mesh conforming to a metric field.
346 
347   Input Parameters:
348 + dm - The DM object
349 . vertexMetric - The metric to which the mesh is adapted, defined vertex-wise in a LOCAL vector
350 - bdLabel - Label for boundary tags which are preserved in dmNew, or NULL. Should not be named "_boundary_".
351 
352   Output Parameter:
353 . dmNew  - the new DM
354 
355   Level: advanced
356 
357 .seealso: DMCoarsen(), DMRefine()
358 */
359 PetscErrorCode DMAdaptMetric_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew)
360 {
361   MPI_Comm           comm;
362   const char        *bdName = "_boundary_";
363 #if 0
364   DM                 odm = dm;
365 #endif
366   DM                 udm, cdm;
367   DMLabel            bdLabelFull;
368   const char        *bdLabelName;
369   IS                 bdIS, globalVertexNum;
370   PetscSection       coordSection;
371   Vec                coordinates;
372   const PetscScalar *coords, *met;
373   const PetscInt    *bdFacesFull, *gV;
374   PetscInt          *bdFaces, *bdFaceIds, *l2gv;
375   PetscReal         *x, *y, *z, *metric;
376   PetscInt          *cells;
377   PetscInt           dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v;
378   PetscInt           off, maxConeSize, numBdFaces, f, bdSize;
379   PetscBool          flg;
380 #ifdef PETSC_HAVE_PRAGMATIC
381   DMLabel            bdLabelNew;
382   double            *coordsNew;
383   PetscInt          *bdTags;
384   PetscReal         *xNew[3] = {NULL, NULL, NULL};
385   PetscInt          *cellsNew;
386   PetscInt           d, numCellsNew, numVerticesNew;
387   PetscInt           numCornersNew, fStart, fEnd;
388 #endif
389   PetscMPIInt        numProcs;
390   PetscErrorCode     ierr;
391 
392   PetscFunctionBegin;
393   /* Check for FEM adjacency flags */
394   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
395   ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
396   if (bdLabel) {
397     ierr = DMLabelGetName(bdLabel, &bdLabelName);CHKERRQ(ierr);
398     ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr);
399     if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
400   }
401   /* Add overlap for Pragmatic */
402 #if 0
403   /* Check for overlap by looking for cell in the SF */
404   if (!overlapped) {
405     ierr = DMPlexDistributeOverlap(odm, 1, NULL, &dm);CHKERRQ(ierr);
406     if (!dm) {dm = odm; ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);}
407   }
408 #endif
409   /* Get mesh information */
410   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
411   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
412   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
413   ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr);
414   ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr);
415   numCells    = cEnd - cStart;
416   numVertices = vEnd - vStart;
417   ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr);
418   for (c = 0, coff = 0; c < numCells; ++c) {
419     const PetscInt *cone;
420     PetscInt        coneSize, cl;
421 
422     ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr);
423     ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr);
424     for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
425   }
426   ierr = PetscCalloc1(numVertices, &l2gv);CHKERRQ(ierr);
427   ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr);
428   ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr);
429   for (v = 0, numLocVertices = 0; v < numVertices; ++v) {
430     if (gV[v] >= 0) ++numLocVertices;
431     l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v];
432   }
433   ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr);
434   ierr = DMDestroy(&udm);CHKERRQ(ierr);
435   ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
436   ierr = DMGetDefaultSection(cdm, &coordSection);CHKERRQ(ierr);
437   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
438   ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr);
439   for (v = vStart; v < vEnd; ++v) {
440     ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
441     x[v-vStart] = PetscRealPart(coords[off+0]);
442     if (dim > 1) y[v-vStart] = PetscRealPart(coords[off+1]);
443     if (dim > 2) z[v-vStart] = PetscRealPart(coords[off+2]);
444   }
445   ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr);
446   /* Get boundary mesh */
447   ierr = DMLabelCreate(bdName, &bdLabelFull);CHKERRQ(ierr);
448   ierr = DMPlexMarkBoundaryFaces(dm, bdLabelFull);CHKERRQ(ierr);
449   ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr);
450   ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr);
451   ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr);
452   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
453     PetscInt *closure = NULL;
454     PetscInt  closureSize, cl;
455 
456     ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
457     for (cl = 0; cl < closureSize*2; cl += 2) {
458       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
459     }
460     ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
461   }
462   ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr);
463   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
464     PetscInt *closure = NULL;
465     PetscInt  closureSize, cl;
466 
467     ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
468     for (cl = 0; cl < closureSize*2; cl += 2) {
469       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart;
470     }
471     ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
472     if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);}
473     else         {bdFaceIds[f] = 1;}
474   }
475   ierr = ISDestroy(&bdIS);CHKERRQ(ierr);
476   ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr);
477   /* Get metric */
478   ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr);
479   for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]);
480   ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr);
481 #if 0
482   /* Destroy overlap mesh */
483   ierr = DMDestroy(&dm);CHKERRQ(ierr);
484 #endif
485   /* Create new mesh */
486 #ifdef PETSC_HAVE_PRAGMATIC
487   switch (dim) {
488   case 2:
489     pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);break;
490   case 3:
491     pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);break;
492   default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim);
493   }
494   pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
495   pragmatic_set_metric(metric);
496   pragmatic_adapt(((DM_Plex *) dm->data)->remeshBd ? 1 : 0);
497   ierr = PetscFree(l2gv);CHKERRQ(ierr);
498   /* Read out mesh */
499   pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew);
500   ierr = PetscMalloc1(numVerticesNew*dim, &coordsNew);CHKERRQ(ierr);
501   switch (dim) {
502   case 2:
503     numCornersNew = 3;
504     ierr = PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);CHKERRQ(ierr);
505     pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]);
506     break;
507   case 3:
508     numCornersNew = 4;
509     ierr = PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);CHKERRQ(ierr);
510     pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]);
511     break;
512   default:
513     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim);
514   }
515   for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = (double) xNew[d][v];}
516   ierr = PetscMalloc1(numCellsNew*(dim+1), &cellsNew);CHKERRQ(ierr);
517   pragmatic_get_elements(cellsNew);
518   ierr = DMPlexCreateFromCellListParallel(comm, dim, numCellsNew, numVerticesNew, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, dmNew);CHKERRQ(ierr);
519   /* Read out boundary label */
520   pragmatic_get_boundaryTags(&bdTags);
521   ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr);
522   ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr);
523   ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr);
524   ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr);
525   ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr);
526   for (c = cStart; c < cEnd; ++c) {
527     /* Only for simplicial meshes */
528     coff = (c-cStart)*(dim+1);
529     /* d is the local cell number of the vertex opposite to the face we are marking */
530     for (d = 0; d < dim+1; ++d) {
531       if (bdTags[coff+d]) {
532         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 */
533         const PetscInt *cone;
534 
535         /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
536         ierr = DMPlexGetCone(*dmNew, c, &cone);CHKERRQ(ierr);
537         ierr = DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);CHKERRQ(ierr);
538       }
539     }
540   }
541   /* Cleanup */
542   switch (dim) {
543   case 2: ierr = PetscFree2(xNew[0], xNew[1]);CHKERRQ(ierr);break;
544   case 3: ierr = PetscFree3(xNew[0], xNew[1], xNew[2]);CHKERRQ(ierr);break;
545   }
546   ierr = PetscFree(cellsNew);CHKERRQ(ierr);
547   ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr);
548   ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr);
549   ierr = PetscFree(coordsNew);CHKERRQ(ierr);
550   pragmatic_finalize();
551 #else
552   SETERRQ(comm, PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic.");
553 #endif
554   PetscFunctionReturn(0);
555 }
556