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