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;
34: PetscFunctionBegin;
36: /* Check for FEM adjacency flags */
37: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
38: PetscCallMPI(MPI_Comm_size(comm, &numProcs));
39: if (bdLabel) {
40: PetscCall(PetscObjectGetName((PetscObject)bdLabel, &bdLabelName));
41: PetscCall(PetscStrcmp(bdLabelName, bdName, &flg));
42: PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
43: }
44: PetscCheck(!rgLabel, comm, PETSC_ERR_ARG_WRONG, "Cannot currently preserve cell tags with Pragmatic");
45: #if 0
46: /* Check for overlap by looking for cell in the SF */
47: if (!overlapped) {
48: PetscCall(DMPlexDistributeOverlap(odm, 1, NULL, &dm));
49: if (!dm) {dm = odm; PetscCall(PetscObjectReference((PetscObject) dm));}
50: }
51: #endif
53: /* Get mesh information */
54: PetscCall(DMGetDimension(dm, &dim));
55: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
56: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
57: PetscCall(DMPlexUninterpolate(dm, &udm));
58: PetscCall(DMPlexGetMaxSizes(udm, &maxConeSize, NULL));
59: numCells = cEnd - cStart;
60: if (numCells == 0) {
61: PetscMPIInt rank;
63: PetscCallMPI(MPI_Comm_rank(comm, &rank));
64: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform mesh adaptation because process %d does not own any cells.", rank);
65: }
66: numVertices = vEnd - vStart;
67: PetscCall(PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices * PetscSqr(dim), &metric, numCells * maxConeSize, &cells));
69: /* Get cell offsets */
70: for (c = 0, coff = 0; c < numCells; ++c) {
71: const PetscInt *cone;
72: PetscInt coneSize, cl;
74: PetscCall(DMPlexGetConeSize(udm, c, &coneSize));
75: PetscCall(DMPlexGetCone(udm, c, &cone));
76: for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
77: }
79: /* Get local-to-global vertex map */
80: PetscCall(PetscCalloc1(numVertices, &l2gv));
81: PetscCall(DMPlexGetVertexNumbering(udm, &globalVertexNum));
82: PetscCall(ISGetIndices(globalVertexNum, &gV));
83: for (v = 0, numLocVertices = 0; v < numVertices; ++v) {
84: if (gV[v] >= 0) ++numLocVertices;
85: l2gv[v] = gV[v] < 0 ? -(gV[v] + 1) : gV[v];
86: }
87: PetscCall(ISRestoreIndices(globalVertexNum, &gV));
88: PetscCall(DMDestroy(&udm));
90: /* Get vertex coordinate arrays */
91: PetscCall(DMGetCoordinateDM(dm, &cdm));
92: PetscCall(DMGetLocalSection(cdm, &coordSection));
93: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
94: PetscCall(VecGetArrayRead(coordinates, &coords));
95: for (v = vStart; v < vEnd; ++v) {
96: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
97: x[v - vStart] = PetscRealPart(coords[off + 0]);
98: if (dim > 1) y[v - vStart] = PetscRealPart(coords[off + 1]);
99: if (dim > 2) z[v - vStart] = PetscRealPart(coords[off + 2]);
100: }
101: PetscCall(VecRestoreArrayRead(coordinates, &coords));
103: /* Get boundary mesh */
104: PetscCall(DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull));
105: PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull));
106: PetscCall(DMLabelGetStratumIS(bdLabelFull, 1, &bdIS));
107: PetscCall(DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces));
108: PetscCall(ISGetIndices(bdIS, &bdFacesFull));
109: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
110: PetscInt *closure = NULL;
111: PetscInt closureSize, cl;
113: PetscCall(DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure));
114: for (cl = 0; cl < closureSize * 2; cl += 2) {
115: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
116: }
117: PetscCall(DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure));
118: }
119: PetscCall(PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds));
120: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
121: PetscInt *closure = NULL;
122: PetscInt closureSize, cl;
124: PetscCall(DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure));
125: for (cl = 0; cl < closureSize * 2; cl += 2) {
126: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart;
127: }
128: PetscCall(DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure));
129: if (bdLabel) PetscCall(DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]));
130: else bdFaceIds[f] = 1;
131: }
132: PetscCall(ISDestroy(&bdIS));
133: PetscCall(DMLabelDestroy(&bdLabelFull));
135: /* Get metric */
136: PetscCall(VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view"));
137: PetscCall(VecGetArrayRead(vertexMetric, &met));
138: PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
139: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
140: Nd = PetscSqr(dim);
141: for (v = 0; v < vEnd - vStart; ++v) {
142: for (i = 0; i < dim; ++i) {
143: for (j = 0; j < dim; ++j) {
144: if (isotropic) {
145: if (i == j) {
146: if (uniform) metric[Nd * v + dim * i + j] = PetscRealPart(met[0]);
147: else metric[Nd * v + dim * i + j] = PetscRealPart(met[v]);
148: } else metric[Nd * v + dim * i + j] = 0.0;
149: } else metric[Nd * v + dim * i + j] = PetscRealPart(met[Nd * v + dim * i + j]);
150: }
151: }
152: }
153: PetscCall(VecRestoreArrayRead(vertexMetric, &met));
155: #if 0
156: /* Destroy overlap mesh */
157: PetscCall(DMDestroy(&dm));
158: #endif
159: /* Send to Pragmatic and remesh */
160: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
161: switch (dim) {
162: case 2:
163: PetscStackCallExternalVoid("pragmatic_2d_mpi_init", pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm));
164: break;
165: case 3:
166: PetscStackCallExternalVoid("pragmatic_3d_mpi_init", pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm));
167: break;
168: default:
169: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim);
170: }
171: PetscStackCallExternalVoid("pragmatic_set_boundary", pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds));
172: PetscStackCallExternalVoid("pragmatic_set_metric", pragmatic_set_metric(metric));
173: PetscStackCallExternalVoid("pragmatic_adapt", pragmatic_adapt(((DM_Plex *)dm->data)->remeshBd ? 1 : 0));
174: PetscCall(PetscFree(l2gv));
175: PetscCall(PetscFPTrapPop());
177: /* Retrieve mesh from Pragmatic and create new plex */
178: PetscStackCallExternalVoid("pragmatic_get_info_mpi", pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew));
179: PetscCall(PetscMalloc1(numVerticesNew * dim, &coordsNew));
180: switch (dim) {
181: case 2:
182: numCornersNew = 3;
183: PetscCall(PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]));
184: PetscStackCallExternalVoid("pragmatic_get_coords_2d_mpi", pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]));
185: break;
186: case 3:
187: numCornersNew = 4;
188: PetscCall(PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]));
189: PetscStackCallExternalVoid("pragmatic_get_coords_3d_mpi", pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]));
190: break;
191: default:
192: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim);
193: }
194: for (v = 0; v < numVerticesNew; ++v) {
195: for (d = 0; d < dim; ++d) coordsNew[v * dim + d] = xNew[d][v];
196: }
197: PetscCall(PetscMalloc1(numCellsNew * (dim + 1), &cellsNew));
198: PetscStackCallExternalVoid("pragmatic_get_elements", pragmatic_get_elements(cellsNew));
199: PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew));
201: /* Rebuild boundary label */
202: PetscStackCallExternalVoid("pragmatic_get_boundaryTags", pragmatic_get_boundaryTags(&bdTags));
203: PetscCall(DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName));
204: PetscCall(DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew));
205: PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd));
206: PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd));
207: PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd));
208: for (c = cStart; c < cEnd; ++c) {
209: /* Only for simplicial meshes */
210: coff = (c - cStart) * (dim + 1);
212: /* d is the local cell number of the vertex opposite to the face we are marking */
213: for (d = 0; d < dim + 1; ++d) {
214: if (bdTags[coff + d]) {
215: const PetscInt perm[4][4] = {
216: {-1, -1, -1, -1},
217: {-1, -1, -1, -1},
218: {1, 2, 0, -1},
219: {3, 2, 1, 0 }
220: }; /* perm[d] = face opposite */
221: const PetscInt *cone;
223: /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
224: PetscCall(DMPlexGetCone(*dmNew, c, &cone));
225: PetscCall(DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff + d]));
226: }
227: }
228: }
230: /* Clean up */
231: switch (dim) {
232: case 2:
233: PetscCall(PetscFree2(xNew[0], xNew[1]));
234: break;
235: case 3:
236: PetscCall(PetscFree3(xNew[0], xNew[1], xNew[2]));
237: break;
238: }
239: PetscCall(PetscFree(cellsNew));
240: PetscCall(PetscFree5(x, y, z, metric, cells));
241: PetscCall(PetscFree2(bdFaces, bdFaceIds));
242: PetscCall(PetscFree(coordsNew));
243: PetscStackCallExternalVoid("pragmatic_finalize", pragmatic_finalize());
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }