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: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim);
167: }
168: pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
169: pragmatic_set_metric(metric);
170: pragmatic_adapt(((DM_Plex *) dm->data)->remeshBd ? 1 : 0);
171: PetscFree(l2gv);
173: /* Retrieve mesh from Pragmatic and create new plex */
174: pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew);
175: PetscMalloc1(numVerticesNew*dim, &coordsNew);
176: switch (dim) {
177: case 2:
178: numCornersNew = 3;
179: PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);
180: pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]);
181: break;
182: case 3:
183: numCornersNew = 4;
184: PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);
185: pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]);
186: break;
187: default:
188: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim);
189: }
190: for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];}
191: PetscMalloc1(numCellsNew*(dim+1), &cellsNew);
192: pragmatic_get_elements(cellsNew);
193: DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew);
195: /* Rebuild boundary label */
196: pragmatic_get_boundaryTags(&bdTags);
197: DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);
198: DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);
199: DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);
200: DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);
201: DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);
202: for (c = cStart; c < cEnd; ++c) {
204: /* Only for simplicial meshes */
205: coff = (c-cStart)*(dim+1);
207: /* d is the local cell number of the vertex opposite to the face we are marking */
208: for (d = 0; d < dim+1; ++d) {
209: if (bdTags[coff+d]) {
210: 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 */
211: const PetscInt *cone;
213: /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
214: DMPlexGetCone(*dmNew, c, &cone);
215: DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);
216: }
217: }
218: }
220: /* Clean up */
221: switch (dim) {
222: case 2: PetscFree2(xNew[0], xNew[1]);
223: break;
224: case 3: PetscFree3(xNew[0], xNew[1], xNew[2]);
225: break;
226: }
227: PetscFree(cellsNew);
228: PetscFree5(x, y, z, metric, cells);
229: PetscFree2(bdFaces, bdFaceIds);
230: PetscFree(coordsNew);
231: pragmatic_finalize();
232: return 0;
233: }