Actual source code: pragmaticadapt.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <pragmatic/cpragmatic.h>
4: PETSC_EXTERN PetscErrorCode DMAdaptMetric_Pragmatic_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DMLabel rgLabel, DM *dmNew)
5: {
6: MPI_Comm comm;
7: const char *bdName = "_boundary_";
8: #if 0
9: DM odm = dm;
10: #endif
11: DM udm, cdm;
12: DMLabel bdLabelFull;
13: const char *bdLabelName;
14: IS bdIS, globalVertexNum;
15: PetscSection coordSection;
16: Vec coordinates;
17: const PetscScalar *coords, *met;
18: const PetscInt *bdFacesFull, *gV;
19: PetscInt *bdFaces, *bdFaceIds, *l2gv;
20: PetscReal *x, *y, *z, *metric;
21: PetscInt *cells;
22: PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v;
23: PetscInt off, maxConeSize, numBdFaces, f, bdSize, i, j, Nd;
24: PetscBool flg, isotropic, uniform;
25: DMLabel bdLabelNew;
26: PetscReal *coordsNew;
27: PetscInt *bdTags;
28: PetscReal *xNew[3] = {NULL, NULL, NULL};
29: PetscInt *cellsNew;
30: PetscInt d, numCellsNew, numVerticesNew;
31: PetscInt numCornersNew, fStart, fEnd;
32: PetscMPIInt numProcs;
35: /* Check for FEM adjacency flags */
36: PetscObjectGetComm((PetscObject)dm, &comm);
37: MPI_Comm_size(comm, &numProcs);
38: if (bdLabel) {
39: PetscObjectGetName((PetscObject)bdLabel, &bdLabelName);
40: PetscStrcmp(bdLabelName, bdName, &flg);
42: }
44: #if 0
45: /* Check for overlap by looking for cell in the SF */
46: if (!overlapped) {
47: DMPlexDistributeOverlap(odm, 1, NULL, &dm);
48: if (!dm) {dm = odm; PetscObjectReference((PetscObject) dm);}
49: }
50: #endif
52: /* Get mesh information */
53: DMGetDimension(dm, &dim);
54: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
55: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
56: DMPlexUninterpolate(dm, &udm);
57: DMPlexGetMaxSizes(udm, &maxConeSize, NULL);
58: numCells = cEnd - cStart;
59: if (numCells == 0) {
60: PetscMPIInt rank;
62: MPI_Comm_rank(comm, &rank);
63: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform mesh adaptation because process %d does not own any cells.", rank);
64: }
65: numVertices = vEnd - vStart;
66: PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices * PetscSqr(dim), &metric, numCells * maxConeSize, &cells);
68: /* Get cell offsets */
69: for (c = 0, coff = 0; c < numCells; ++c) {
70: const PetscInt *cone;
71: PetscInt coneSize, cl;
73: DMPlexGetConeSize(udm, c, &coneSize);
74: DMPlexGetCone(udm, c, &cone);
75: for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
76: }
78: /* Get local-to-global vertex map */
79: PetscCalloc1(numVertices, &l2gv);
80: DMPlexGetVertexNumbering(udm, &globalVertexNum);
81: ISGetIndices(globalVertexNum, &gV);
82: for (v = 0, numLocVertices = 0; v < numVertices; ++v) {
83: if (gV[v] >= 0) ++numLocVertices;
84: l2gv[v] = gV[v] < 0 ? -(gV[v] + 1) : gV[v];
85: }
86: ISRestoreIndices(globalVertexNum, &gV);
87: DMDestroy(&udm);
89: /* Get vertex coordinate arrays */
90: DMGetCoordinateDM(dm, &cdm);
91: DMGetLocalSection(cdm, &coordSection);
92: DMGetCoordinatesLocal(dm, &coordinates);
93: VecGetArrayRead(coordinates, &coords);
94: for (v = vStart; v < vEnd; ++v) {
95: PetscSectionGetOffset(coordSection, v, &off);
96: x[v - vStart] = PetscRealPart(coords[off + 0]);
97: if (dim > 1) y[v - vStart] = PetscRealPart(coords[off + 1]);
98: if (dim > 2) z[v - vStart] = PetscRealPart(coords[off + 2]);
99: }
100: VecRestoreArrayRead(coordinates, &coords);
102: /* Get boundary mesh */
103: DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);
104: DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);
105: DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);
106: DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);
107: ISGetIndices(bdIS, &bdFacesFull);
108: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
109: PetscInt *closure = NULL;
110: PetscInt closureSize, cl;
112: DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
113: for (cl = 0; cl < closureSize * 2; cl += 2) {
114: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
115: }
116: DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
117: }
118: PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);
119: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
120: PetscInt *closure = NULL;
121: PetscInt closureSize, cl;
123: DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
124: for (cl = 0; cl < closureSize * 2; cl += 2) {
125: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart;
126: }
127: DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);
128: if (bdLabel) DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);
129: else bdFaceIds[f] = 1;
130: }
131: ISDestroy(&bdIS);
132: DMLabelDestroy(&bdLabelFull);
134: /* Get metric */
135: VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");
136: VecGetArrayRead(vertexMetric, &met);
137: DMPlexMetricIsIsotropic(dm, &isotropic);
138: DMPlexMetricIsUniform(dm, &uniform);
139: Nd = PetscSqr(dim);
140: for (v = 0; v < vEnd - vStart; ++v) {
141: for (i = 0; i < dim; ++i) {
142: for (j = 0; j < dim; ++j) {
143: if (isotropic) {
144: if (i == j) {
145: if (uniform) metric[Nd * v + dim * i + j] = PetscRealPart(met[0]);
146: else metric[Nd * v + dim * i + j] = PetscRealPart(met[v]);
147: } else metric[Nd * v + dim * i + j] = 0.0;
148: } else metric[Nd * v + dim * i + j] = PetscRealPart(met[Nd * v + dim * i + j]);
149: }
150: }
151: }
152: VecRestoreArrayRead(vertexMetric, &met);
154: #if 0
155: /* Destroy overlap mesh */
156: DMDestroy(&dm);
157: #endif
158: /* Send to Pragmatic and remesh */
159: switch (dim) {
160: case 2:
161: pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);
162: break;
163: case 3:
164: pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);
165: break;
166: default:
167: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim);
168: }
169: pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
170: pragmatic_set_metric(metric);
171: pragmatic_adapt(((DM_Plex *)dm->data)->remeshBd ? 1 : 0);
172: PetscFree(l2gv);
174: /* Retrieve mesh from Pragmatic and create new plex */
175: pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew);
176: PetscMalloc1(numVerticesNew * dim, &coordsNew);
177: switch (dim) {
178: case 2:
179: numCornersNew = 3;
180: PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);
181: pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]);
182: break;
183: case 3:
184: numCornersNew = 4;
185: PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);
186: pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]);
187: break;
188: default:
189: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim);
190: }
191: for (v = 0; v < numVerticesNew; ++v) {
192: for (d = 0; d < dim; ++d) coordsNew[v * dim + d] = xNew[d][v];
193: }
194: PetscMalloc1(numCellsNew * (dim + 1), &cellsNew);
195: pragmatic_get_elements(cellsNew);
196: DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew);
198: /* Rebuild boundary label */
199: pragmatic_get_boundaryTags(&bdTags);
200: DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);
201: DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);
202: DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);
203: DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);
204: DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);
205: for (c = cStart; c < cEnd; ++c) {
206: /* Only for simplicial meshes */
207: coff = (c - cStart) * (dim + 1);
209: /* d is the local cell number of the vertex opposite to the face we are marking */
210: for (d = 0; d < dim + 1; ++d) {
211: if (bdTags[coff + d]) {
212: const PetscInt perm[4][4] = {
213: {-1, -1, -1, -1},
214: {-1, -1, -1, -1},
215: {1, 2, 0, -1},
216: {3, 2, 1, 0 }
217: }; /* perm[d] = face opposite */
218: const PetscInt *cone;
220: /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
221: DMPlexGetCone(*dmNew, c, &cone);
222: DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff + d]);
223: }
224: }
225: }
227: /* Clean up */
228: switch (dim) {
229: case 2:
230: PetscFree2(xNew[0], xNew[1]);
231: break;
232: case 3:
233: PetscFree3(xNew[0], xNew[1], xNew[2]);
234: break;
235: }
236: PetscFree(cellsNew);
237: PetscFree5(x, y, z, metric, cells);
238: PetscFree2(bdFaces, bdFaceIds);
239: PetscFree(coordsNew);
240: pragmatic_finalize();
241: return 0;
242: }