Actual source code: plexadapt.c
petsc-3.9.4 2018-09-11
1: #include <petsc/private/dmpleximpl.h>
2: #ifdef PETSC_HAVE_PRAGMATIC
3: #include <pragmatic/cpragmatic.h>
4: #endif
6: static PetscErrorCode DMPlexLabelToVolumeConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscReal refRatio, PetscReal maxVolumes[])
7: {
8: PetscInt dim, c;
12: DMGetDimension(dm, &dim);
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;
22: DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);
23: maxVolumes[c - cStart] = vol;
24: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
25: for (cl = 0; cl < closureSize*2; cl += 2) {
26: const PetscInt point = closure[cl];
27: PetscInt refFlag;
29: DMLabelGetValue(adaptLabel, point, &refFlag);
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: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
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: return(0);
54: }
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;
68: DMPlexUninterpolate(dm, &udm);
69: DMGetDimension(dm, &dim);
70: DMGetCoordinateDM(dm, &coordDM);
71: DMGetDefaultSection(coordDM, &coordSection);
72: DMGetCoordinatesLocal(dm, &coordinates);
73: Nv = vEnd - vStart;
74: VecCreateSeq(PETSC_COMM_SELF, Nv*PetscSqr(dim), metricVec);
75: VecGetArray(*metricVec, &metric);
76: Neq = (dim*(dim+1))/2;
77: PetscMalloc1(PetscSqr(Neq), &eqns);
78: MatCreateSeqDense(PETSC_COMM_SELF, Neq, Neq, eqns, &A);
79: MatCreateVecs(A, &mx, &mb);
80: VecSet(mb, 1.0);
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;
88: DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);
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: MatSetUnfactored(A);
112: DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);
113: MatLUFactor(A, NULL, NULL, NULL);
114: MatSolve(A, mb, mx);
115: VecGetArrayRead(mx, &sol);
116: DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);
117: DMPlexGetCone(udm, c, &cone);
118: DMPlexGetConeSize(udm, c, &coneSize);
119: for (cl = 0; cl < coneSize; ++cl) {
120: const PetscInt v = cone[cl] - vStart;
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: VecRestoreArrayRead(mx, &sol);
140: }
141: for (v = 0; v < Nv; ++v) {
142: const PetscInt *support;
143: PetscInt supportSize, s;
144: PetscReal vol, totVol = 0.0;
146: DMPlexGetSupport(udm, v+vStart, &support);
147: DMPlexGetSupportSize(udm, v+vStart, &supportSize);
148: for (s = 0; s < supportSize; ++s) {DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL); totVol += vol;}
149: for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol;
150: }
151: PetscFree(eqns);
152: VecRestoreArray(*metricVec, &metric);
153: VecDestroy(&mx);
154: VecDestroy(&mb);
155: MatDestroy(&A);
156: DMDestroy(&udm);
157: return(0);
158: }
160: /*
161: Contains the list of registered DMPlexGenerators routines
162: */
163: extern PetscFunctionList DMPlexGenerateList;
165: struct _n_PetscFunctionList {
166: PetscErrorCode (*generate)(DM, PetscBool, DM*);
167: PetscErrorCode (*refine)(DM,PetscReal*, DM*);
168: char *name; /* string to identify routine */
169: PetscInt dim;
170: PetscFunctionList next; /* next pointer */
171: };
173: PetscErrorCode DMPlexRefine_Internal(DM dm, DMLabel adaptLabel, DM *dmRefined)
174: {
175: PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *);
176: PetscReal refinementLimit;
177: PetscInt dim, cStart, cEnd;
178: char genname[1024], *name = NULL;
179: PetscBool flg, localized;
180: PetscErrorCode ierr;
181: PetscErrorCode (*refine)(DM,PetscReal*,DM*);
182: PetscFunctionList fl;
183: PetscReal *maxVolumes;
184: PetscInt c;
187: DMGetCoordinatesLocalized(dm, &localized);
188: DMPlexGetRefinementLimit(dm, &refinementLimit);
189: DMPlexGetRefinementFunction(dm, &refinementFunc);
190: if (refinementLimit == 0.0 && !refinementFunc && !adaptLabel) return(0);
191: DMGetDimension(dm, &dim);
192: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
193: PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, 1024, &flg);
194: if (flg) name = genname;
196: fl = DMPlexGenerateList;
197: if (name) {
198: while (fl) {
199: PetscStrcmp(fl->name,name,&flg);
200: if (flg) {
201: refine = fl->refine;
202: goto gotit;
203: }
204: fl = fl->next;
205: }
206: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Grid refiner %g not registered",name);
207: } else {
208: while (fl) {
209: if (dim-1 == fl->dim) {
210: refine = fl->refine;
211: goto gotit;
212: }
213: fl = fl->next;
214: }
215: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No grid refiner of dimension %D registered",dim);
216: }
218: gotit: switch (dim) {
219: case 2:
220: PetscMalloc1(cEnd - cStart, &maxVolumes);
221: if (adaptLabel) {
222: DMPlexLabelToVolumeConstraint(dm, adaptLabel, cStart, cEnd, PETSC_DEFAULT, maxVolumes);
223: } else if (refinementFunc) {
224: for (c = cStart; c < cEnd; ++c) {
225: PetscReal vol, centroid[3];
226: PetscReal maxVol;
228: DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
229: (*refinementFunc)(centroid, &maxVol);
230: maxVolumes[c - cStart] = (double) maxVol;
231: }
232: } else {
233: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
234: }
235: (*refine)(dm, maxVolumes, dmRefined);
236: PetscFree(maxVolumes);
237: break;
238: case 3:
239: PetscMalloc1(cEnd - cStart, &maxVolumes);
240: if (adaptLabel) {
241: DMPlexLabelToVolumeConstraint(dm, adaptLabel, cStart, cEnd, PETSC_DEFAULT, maxVolumes);
242: } else if (refinementFunc) {
243: for (c = cStart; c < cEnd; ++c) {
244: PetscReal vol, centroid[3];
246: DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
247: (*refinementFunc)(centroid, &maxVolumes[c-cStart]);
248: }
249: } else {
250: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
251: }
252: (*refine)(dm, maxVolumes, dmRefined);
253: PetscFree(maxVolumes);
254: break;
255: default:
256: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
257: }
258: DMCopyBoundary(dm, *dmRefined);
259: if (localized) {DMLocalizeCoordinates(*dmRefined);}
260: return(0);
261: }
263: PetscErrorCode DMPlexCoarsen_Internal(DM dm, DMLabel adaptLabel, DM *dmCoarsened)
264: {
265: Vec metricVec;
266: PetscInt cStart, cEnd, vStart, vEnd;
267: DMLabel bdLabel = NULL;
268: char bdLabelName[PETSC_MAX_PATH_LEN];
269: PetscBool localized, flg;
273: DMGetCoordinatesLocalized(dm, &localized);
274: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
275: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
276: DMPlexLabelToMetricConstraint(dm, adaptLabel, cStart, cEnd, vStart, vEnd, PETSC_DEFAULT, &metricVec);
277: PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_bd_label", bdLabelName, PETSC_MAX_PATH_LEN-1, &flg);
278: if (flg) {DMGetLabel(dm, bdLabelName, &bdLabel);}
279: DMAdaptMetric_Plex(dm, metricVec, bdLabel, dmCoarsened);
280: VecDestroy(&metricVec);
281: if (localized) {DMLocalizeCoordinates(*dmCoarsened);}
282: return(0);
283: }
285: PetscErrorCode DMAdaptLabel_Plex(DM dm, DMLabel adaptLabel, DM *dmAdapted)
286: {
287: IS flagIS;
288: const PetscInt *flags;
289: PetscInt defFlag, minFlag, maxFlag, numFlags, f;
290: PetscErrorCode ierr;
293: DMLabelGetDefaultValue(adaptLabel, &defFlag);
294: minFlag = defFlag;
295: maxFlag = defFlag;
296: DMLabelGetValueIS(adaptLabel, &flagIS);
297: ISGetLocalSize(flagIS, &numFlags);
298: ISGetIndices(flagIS, &flags);
299: for (f = 0; f < numFlags; ++f) {
300: const PetscInt flag = flags[f];
302: minFlag = PetscMin(minFlag, flag);
303: maxFlag = PetscMax(maxFlag, flag);
304: }
305: ISRestoreIndices(flagIS, &flags);
306: ISDestroy(&flagIS);
307: {
308: PetscInt minMaxFlag[2], minMaxFlagGlobal[2];
310: minMaxFlag[0] = minFlag;
311: minMaxFlag[1] = -maxFlag;
312: MPI_Allreduce(minMaxFlag, minMaxFlagGlobal, 2, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
313: minFlag = minMaxFlagGlobal[0];
314: maxFlag = -minMaxFlagGlobal[1];
315: }
316: if (minFlag == maxFlag) {
317: switch (minFlag) {
318: case DM_ADAPT_DETERMINE:
319: *dmAdapted = NULL;break;
320: case DM_ADAPT_REFINE:
321: DMPlexSetRefinementUniform(dm, PETSC_TRUE);
322: DMRefine(dm, MPI_COMM_NULL, dmAdapted);break;
323: case DM_ADAPT_COARSEN:
324: DMCoarsen(dm, MPI_COMM_NULL, dmAdapted);break;
325: default:
326: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP,"DMPlex does not support refinement flag %D\n", minFlag);break;
327: }
328: } else {
329: DMPlexSetRefinementUniform(dm, PETSC_FALSE);
330: DMPlexRefine_Internal(dm, adaptLabel, dmAdapted);
331: }
332: return(0);
333: }
335: /*
336: DMAdaptMetric_Plex - Generates a new mesh conforming to a metric field.
338: Input Parameters:
339: + dm - The DM object
340: . vertexMetric - The metric to which the mesh is adapted, defined vertex-wise in a LOCAL vector
341: - bdLabel - Label for boundary tags which are preserved in dmNew, or NULL. Should not be named "_boundary_".
343: Output Parameter:
344: . dmNew - the new DM
346: Level: advanced
348: .seealso: DMCoarsen(), DMRefine()
349: */
350: PetscErrorCode DMAdaptMetric_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew)
351: {
352: #ifdef PETSC_HAVE_PRAGMATIC
353: MPI_Comm comm;
354: const char *bdName = "_boundary_";
355: #if 0
356: DM odm = dm;
357: #endif
358: DM udm, cdm;
359: DMLabel bdLabelFull;
360: const char *bdLabelName;
361: IS bdIS, globalVertexNum;
362: PetscSection coordSection;
363: Vec coordinates;
364: const PetscScalar *coords, *met;
365: const PetscInt *bdFacesFull, *gV;
366: PetscInt *bdFaces, *bdFaceIds, *l2gv;
367: PetscReal *x, *y, *z, *metric;
368: PetscInt *cells;
369: PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v;
370: PetscInt off, maxConeSize, numBdFaces, f, bdSize;
371: PetscBool flg;
372: DMLabel bdLabelNew;
373: double *coordsNew;
374: PetscInt *bdTags;
375: PetscReal *xNew[3] = {NULL, NULL, NULL};
376: PetscInt *cellsNew;
377: PetscInt d, numCellsNew, numVerticesNew;
378: PetscInt numCornersNew, fStart, fEnd;
379: PetscMPIInt numProcs;
380: PetscErrorCode ierr;
383: /* Check for FEM adjacency flags */
384: PetscObjectGetComm((PetscObject) dm, &comm);
385: MPI_Comm_size(comm, &numProcs);
386: if (bdLabel) {
387: DMLabelGetName(bdLabel, &bdLabelName);
388: PetscStrcmp(bdLabelName, bdName, &flg);
389: if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
390: }
391: /* Add overlap for Pragmatic */
392: #if 0
393: /* Check for overlap by looking for cell in the SF */
394: if (!overlapped) {
395: DMPlexDistributeOverlap(odm, 1, NULL, &dm);
396: if (!dm) {dm = odm; PetscObjectReference((PetscObject) dm);}
397: }
398: #endif
399: /* Get mesh information */
400: DMGetDimension(dm, &dim);
401: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
402: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
403: DMPlexUninterpolate(dm, &udm);
404: DMPlexGetMaxSizes(udm, &maxConeSize, NULL);
405: numCells = cEnd - cStart;
406: numVertices = vEnd - vStart;
407: PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);
408: for (c = 0, coff = 0; c < numCells; ++c) {
409: const PetscInt *cone;
410: PetscInt coneSize, cl;
412: DMPlexGetConeSize(udm, c, &coneSize);
413: DMPlexGetCone(udm, c, &cone);
414: for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
415: }
416: PetscCalloc1(numVertices, &l2gv);
417: DMPlexGetVertexNumbering(udm, &globalVertexNum);
418: ISGetIndices(globalVertexNum, &gV);
419: for (v = 0, numLocVertices = 0; v < numVertices; ++v) {
420: if (gV[v] >= 0) ++numLocVertices;
421: l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v];
422: }
423: ISRestoreIndices(globalVertexNum, &gV);
424: DMDestroy(&udm);
425: DMGetCoordinateDM(dm, &cdm);
426: DMGetDefaultSection(cdm, &coordSection);
427: DMGetCoordinatesLocal(dm, &coordinates);
428: VecGetArrayRead(coordinates, &coords);
429: for (v = vStart; v < vEnd; ++v) {
430: PetscSectionGetOffset(coordSection, v, &off);
431: x[v-vStart] = PetscRealPart(coords[off+0]);
432: if (dim > 1) y[v-vStart] = PetscRealPart(coords[off+1]);
433: if (dim > 2) z[v-vStart] = PetscRealPart(coords[off+2]);
434: }
435: VecRestoreArrayRead(coordinates, &coords);
436: /* Get boundary mesh */
437: DMLabelCreate(bdName, &bdLabelFull);
438: DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);
439: DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);
440: DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);
441: ISGetIndices(bdIS, &bdFacesFull);
442: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
443: PetscInt *closure = NULL;
444: PetscInt closureSize, cl;
446: DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
447: for (cl = 0; cl < closureSize*2; cl += 2) {
448: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
449: }
450: DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
451: }
452: PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);
453: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
454: PetscInt *closure = NULL;
455: PetscInt closureSize, cl;
457: DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
458: for (cl = 0; cl < closureSize*2; cl += 2) {
459: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart;
460: }
461: DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
462: if (bdLabel) {DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);}
463: else {bdFaceIds[f] = 1;}
464: }
465: ISDestroy(&bdIS);
466: DMLabelDestroy(&bdLabelFull);
467: /* Get metric */
468: VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");
469: VecGetArrayRead(vertexMetric, &met);
470: for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]);
471: VecRestoreArrayRead(vertexMetric, &met);
472: #if 0
473: /* Destroy overlap mesh */
474: DMDestroy(&dm);
475: #endif
476: /* Create new mesh */
477: switch (dim) {
478: case 2:
479: pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);break;
480: case 3:
481: pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);break;
482: default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim);
483: }
484: pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
485: pragmatic_set_metric(metric);
486: pragmatic_adapt(((DM_Plex *) dm->data)->remeshBd ? 1 : 0);
487: PetscFree(l2gv);
488: /* Read out mesh */
489: pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew);
490: PetscMalloc1(numVerticesNew*dim, &coordsNew);
491: switch (dim) {
492: case 2:
493: numCornersNew = 3;
494: PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);
495: pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]);
496: break;
497: case 3:
498: numCornersNew = 4;
499: PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);
500: pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]);
501: break;
502: default:
503: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim);
504: }
505: for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = (double) xNew[d][v];}
506: PetscMalloc1(numCellsNew*(dim+1), &cellsNew);
507: pragmatic_get_elements(cellsNew);
508: DMPlexCreateFromCellListParallel(comm, dim, numCellsNew, numVerticesNew, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, dmNew);
509: /* Read out boundary label */
510: pragmatic_get_boundaryTags(&bdTags);
511: DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);
512: DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);
513: DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);
514: DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);
515: DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);
516: for (c = cStart; c < cEnd; ++c) {
517: /* Only for simplicial meshes */
518: coff = (c-cStart)*(dim+1);
519: /* d is the local cell number of the vertex opposite to the face we are marking */
520: for (d = 0; d < dim+1; ++d) {
521: if (bdTags[coff+d]) {
522: 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 */
523: const PetscInt *cone;
525: /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
526: DMPlexGetCone(*dmNew, c, &cone);
527: DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);
528: }
529: }
530: }
531: /* Cleanup */
532: switch (dim) {
533: case 2: PetscFree2(xNew[0], xNew[1]);break;
534: case 3: PetscFree3(xNew[0], xNew[1], xNew[2]);break;
535: }
536: PetscFree(cellsNew);
537: PetscFree5(x, y, z, metric, cells);
538: PetscFree2(bdFaces, bdFaceIds);
539: PetscFree(coordsNew);
540: pragmatic_finalize();
541: return(0);
542: #else
544: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic.");
545: return(0);
546: #endif
547: }