Actual source code: plexadapt.c
1: #include <petsc/private/dmpleximpl.h>
2: #if defined(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);
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: DMGetLocalSection(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 PlexGeneratorFunctionList DMPlexGenerateList;
165: PetscErrorCode DMPlexRefine_Internal(DM dm, DMLabel adaptLabel, DM *dmRefined)
166: {
167: PlexGeneratorFunctionList fl;
168: PetscErrorCode (*refine)(DM,PetscReal*,DM*);
169: PetscErrorCode (*adapt)(DM,DMLabel,DM*);
170: PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *);
171: char genname[PETSC_MAX_PATH_LEN], *name = NULL;
172: PetscReal refinementLimit;
173: PetscReal *maxVolumes;
174: PetscInt dim, cStart, cEnd, c;
175: PetscBool flg, flg2, localized;
176: PetscErrorCode ierr;
179: DMGetCoordinatesLocalized(dm, &localized);
180: DMPlexGetRefinementLimit(dm, &refinementLimit);
181: DMPlexGetRefinementFunction(dm, &refinementFunc);
182: if (refinementLimit == 0.0 && !refinementFunc && !adaptLabel) return(0);
183: DMGetDimension(dm, &dim);
184: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
185: PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_adaptor", genname, sizeof(genname), &flg);
186: if (flg) name = genname;
187: else {
188: PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, sizeof(genname), &flg2);
189: if (flg2) name = genname;
190: }
192: fl = DMPlexGenerateList;
193: if (name) {
194: while (fl) {
195: PetscStrcmp(fl->name,name,&flg);
196: if (flg) {
197: refine = fl->refine;
198: adapt = fl->adaptlabel;
199: goto gotit;
200: }
201: fl = fl->next;
202: }
203: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Grid refiner %s not registered",name);
204: } else {
205: while (fl) {
206: if (dim-1 == fl->dim) {
207: refine = fl->refine;
208: adapt = fl->adaptlabel;
209: goto gotit;
210: }
211: fl = fl->next;
212: }
213: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No grid refiner of dimension %D registered",dim);
214: }
216: gotit:
217: switch (dim) {
218: case 2:
219: case 3:
220: if (adapt) {
221: (*adapt)(dm, adaptLabel, dmRefined);
222: } else {
223: PetscMalloc1(cEnd - cStart, &maxVolumes);
224: if (adaptLabel) {
225: DMPlexLabelToVolumeConstraint(dm, adaptLabel, cStart, cEnd, PETSC_DEFAULT, maxVolumes);
226: } else if (refinementFunc) {
227: for (c = cStart; c < cEnd; ++c) {
228: PetscReal vol, centroid[3];
230: DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);
231: (*refinementFunc)(centroid, &maxVolumes[c-cStart]);
232: }
233: } else {
234: for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit;
235: }
236: (*refine)(dm, maxVolumes, dmRefined);
237: PetscFree(maxVolumes);
238: }
239: break;
240: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %D is not supported.", dim);
241: }
242: DMCopyDisc(dm, *dmRefined);
243: if (localized) {DMLocalizeCoordinates(*dmRefined);}
244: return(0);
245: }
247: PetscErrorCode DMPlexCoarsen_Internal(DM dm, DMLabel adaptLabel, DM *dmCoarsened)
248: {
249: Vec metricVec;
250: PetscInt cStart, cEnd, vStart, vEnd;
251: DMLabel bdLabel = NULL;
252: char bdLabelName[PETSC_MAX_PATH_LEN];
253: PetscBool localized, flg;
257: DMGetCoordinatesLocalized(dm, &localized);
258: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
259: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
260: DMPlexLabelToMetricConstraint(dm, adaptLabel, cStart, cEnd, vStart, vEnd, PETSC_DEFAULT, &metricVec);
261: PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_bd_label", bdLabelName, sizeof(bdLabelName), &flg);
262: if (flg) {DMGetLabel(dm, bdLabelName, &bdLabel);}
263: DMAdaptMetric_Plex(dm, metricVec, bdLabel, dmCoarsened);
264: VecDestroy(&metricVec);
265: if (localized) {DMLocalizeCoordinates(*dmCoarsened);}
266: return(0);
267: }
269: PetscErrorCode DMAdaptLabel_Plex(DM dm, DMLabel adaptLabel, DM *dmAdapted)
270: {
271: IS flagIS;
272: const PetscInt *flags;
273: PetscInt defFlag, minFlag, maxFlag, numFlags, f;
274: PetscErrorCode ierr;
277: DMLabelGetDefaultValue(adaptLabel, &defFlag);
278: minFlag = defFlag;
279: maxFlag = defFlag;
280: DMLabelGetValueIS(adaptLabel, &flagIS);
281: ISGetLocalSize(flagIS, &numFlags);
282: ISGetIndices(flagIS, &flags);
283: for (f = 0; f < numFlags; ++f) {
284: const PetscInt flag = flags[f];
286: minFlag = PetscMin(minFlag, flag);
287: maxFlag = PetscMax(maxFlag, flag);
288: }
289: ISRestoreIndices(flagIS, &flags);
290: ISDestroy(&flagIS);
291: {
292: PetscInt minMaxFlag[2], minMaxFlagGlobal[2];
294: minMaxFlag[0] = minFlag;
295: minMaxFlag[1] = -maxFlag;
296: MPI_Allreduce(minMaxFlag, minMaxFlagGlobal, 2, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
297: minFlag = minMaxFlagGlobal[0];
298: maxFlag = -minMaxFlagGlobal[1];
299: }
300: if (minFlag == maxFlag) {
301: switch (minFlag) {
302: case DM_ADAPT_DETERMINE:
303: *dmAdapted = NULL;break;
304: case DM_ADAPT_REFINE:
305: DMPlexSetRefinementUniform(dm, PETSC_TRUE);
306: DMRefine(dm, MPI_COMM_NULL, dmAdapted);break;
307: case DM_ADAPT_COARSEN:
308: DMCoarsen(dm, MPI_COMM_NULL, dmAdapted);break;
309: default:
310: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP,"DMPlex does not support refinement flag %D\n", minFlag);
311: }
312: } else {
313: DMPlexSetRefinementUniform(dm, PETSC_FALSE);
314: DMPlexRefine_Internal(dm, adaptLabel, dmAdapted);
315: }
316: return(0);
317: }
319: /*
320: DMAdaptMetric_Plex - Generates a new mesh conforming to a metric field.
322: Input Parameters:
323: + dm - The DM object
324: . vertexMetric - The metric to which the mesh is adapted, defined vertex-wise in a LOCAL vector
325: - bdLabel - Label for boundary tags which are preserved in dmNew, or NULL. Should not be named "_boundary_".
327: Output Parameter:
328: . dmNew - the new DM
330: Level: advanced
332: .seealso: DMCoarsen(), DMRefine()
333: */
334: PetscErrorCode DMAdaptMetric_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew)
335: {
336: #if defined(PETSC_HAVE_PRAGMATIC)
337: MPI_Comm comm;
338: const char *bdName = "_boundary_";
339: #if 0
340: DM odm = dm;
341: #endif
342: DM udm, cdm;
343: DMLabel bdLabelFull;
344: const char *bdLabelName;
345: IS bdIS, globalVertexNum;
346: PetscSection coordSection;
347: Vec coordinates;
348: const PetscScalar *coords, *met;
349: const PetscInt *bdFacesFull, *gV;
350: PetscInt *bdFaces, *bdFaceIds, *l2gv;
351: PetscReal *x, *y, *z, *metric;
352: PetscInt *cells;
353: PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v;
354: PetscInt off, maxConeSize, numBdFaces, f, bdSize;
355: PetscBool flg;
356: DMLabel bdLabelNew;
357: PetscReal *coordsNew;
358: PetscInt *bdTags;
359: PetscReal *xNew[3] = {NULL, NULL, NULL};
360: PetscInt *cellsNew;
361: PetscInt d, numCellsNew, numVerticesNew;
362: PetscInt numCornersNew, fStart, fEnd;
363: PetscMPIInt numProcs;
364: PetscErrorCode ierr;
367: /* Check for FEM adjacency flags */
368: PetscObjectGetComm((PetscObject) dm, &comm);
369: MPI_Comm_size(comm, &numProcs);
370: if (bdLabel) {
371: PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);
372: PetscStrcmp(bdLabelName, bdName, &flg);
373: if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
374: }
375: /* Add overlap for Pragmatic */
376: #if 0
377: /* Check for overlap by looking for cell in the SF */
378: if (!overlapped) {
379: DMPlexDistributeOverlap(odm, 1, NULL, &dm);
380: if (!dm) {dm = odm; PetscObjectReference((PetscObject) dm);}
381: }
382: #endif
383: /* Get mesh information */
384: DMGetDimension(dm, &dim);
385: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
386: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
387: DMPlexUninterpolate(dm, &udm);
388: DMPlexGetMaxSizes(udm, &maxConeSize, NULL);
389: numCells = cEnd - cStart;
390: numVertices = vEnd - vStart;
391: PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);
392: for (c = 0, coff = 0; c < numCells; ++c) {
393: const PetscInt *cone;
394: PetscInt coneSize, cl;
396: DMPlexGetConeSize(udm, c, &coneSize);
397: DMPlexGetCone(udm, c, &cone);
398: for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
399: }
400: PetscCalloc1(numVertices, &l2gv);
401: DMPlexGetVertexNumbering(udm, &globalVertexNum);
402: ISGetIndices(globalVertexNum, &gV);
403: for (v = 0, numLocVertices = 0; v < numVertices; ++v) {
404: if (gV[v] >= 0) ++numLocVertices;
405: l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v];
406: }
407: ISRestoreIndices(globalVertexNum, &gV);
408: DMDestroy(&udm);
409: DMGetCoordinateDM(dm, &cdm);
410: DMGetLocalSection(cdm, &coordSection);
411: DMGetCoordinatesLocal(dm, &coordinates);
412: VecGetArrayRead(coordinates, &coords);
413: for (v = vStart; v < vEnd; ++v) {
414: PetscSectionGetOffset(coordSection, v, &off);
415: x[v-vStart] = PetscRealPart(coords[off+0]);
416: if (dim > 1) y[v-vStart] = PetscRealPart(coords[off+1]);
417: if (dim > 2) z[v-vStart] = PetscRealPart(coords[off+2]);
418: }
419: VecRestoreArrayRead(coordinates, &coords);
420: /* Get boundary mesh */
421: DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);
422: DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);
423: DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);
424: DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);
425: ISGetIndices(bdIS, &bdFacesFull);
426: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
427: PetscInt *closure = NULL;
428: PetscInt closureSize, cl;
430: DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
431: for (cl = 0; cl < closureSize*2; cl += 2) {
432: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
433: }
434: DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
435: }
436: PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);
437: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
438: PetscInt *closure = NULL;
439: PetscInt closureSize, cl;
441: DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
442: for (cl = 0; cl < closureSize*2; cl += 2) {
443: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart;
444: }
445: DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
446: if (bdLabel) {DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);}
447: else {bdFaceIds[f] = 1;}
448: }
449: ISDestroy(&bdIS);
450: DMLabelDestroy(&bdLabelFull);
451: /* Get metric */
452: VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");
453: VecGetArrayRead(vertexMetric, &met);
454: for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]);
455: VecRestoreArrayRead(vertexMetric, &met);
456: #if 0
457: /* Destroy overlap mesh */
458: DMDestroy(&dm);
459: #endif
460: /* Create new mesh */
461: switch (dim) {
462: case 2:
463: pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);break;
464: case 3:
465: pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);break;
466: default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim);
467: }
468: pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
469: pragmatic_set_metric(metric);
470: pragmatic_adapt(((DM_Plex *) dm->data)->remeshBd ? 1 : 0);
471: PetscFree(l2gv);
472: /* Read out mesh */
473: pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew);
474: PetscMalloc1(numVerticesNew*dim, &coordsNew);
475: switch (dim) {
476: case 2:
477: numCornersNew = 3;
478: PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);
479: pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]);
480: break;
481: case 3:
482: numCornersNew = 4;
483: PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);
484: pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]);
485: break;
486: default:
487: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim);
488: }
489: for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];}
490: PetscMalloc1(numCellsNew*(dim+1), &cellsNew);
491: pragmatic_get_elements(cellsNew);
492: DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, dmNew);
493: /* Read out boundary label */
494: pragmatic_get_boundaryTags(&bdTags);
495: DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);
496: DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);
497: DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);
498: DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);
499: DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);
500: for (c = cStart; c < cEnd; ++c) {
501: /* Only for simplicial meshes */
502: coff = (c-cStart)*(dim+1);
503: /* d is the local cell number of the vertex opposite to the face we are marking */
504: for (d = 0; d < dim+1; ++d) {
505: if (bdTags[coff+d]) {
506: 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 */
507: const PetscInt *cone;
509: /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
510: DMPlexGetCone(*dmNew, c, &cone);
511: DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);
512: }
513: }
514: }
515: /* Cleanup */
516: switch (dim) {
517: case 2: PetscFree2(xNew[0], xNew[1]);break;
518: case 3: PetscFree3(xNew[0], xNew[1], xNew[2]);break;
519: }
520: PetscFree(cellsNew);
521: PetscFree5(x, y, z, metric, cells);
522: PetscFree2(bdFaces, bdFaceIds);
523: PetscFree(coordsNew);
524: pragmatic_finalize();
525: return(0);
526: #else
528: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic.");
529: #endif
530: }