Actual source code: plexcoarsen.c
petsc-3.7.7 2017-09-25
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #ifdef PETSC_HAVE_PRAGMATIC
3: #include <pragmatic/cpragmatic.h>
4: #endif
9: PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened)
10: {
11: #ifdef PETSC_HAVE_PRAGMATIC
12: DM udm, coordDM;
13: DMLabel bd;
14: Mat A;
15: Vec coordinates, mb, mx;
16: PetscSection coordSection;
17: const PetscScalar *coords;
18: double *coarseCoords;
19: IS bdIS;
20: PetscReal *x, *y, *z, *eqns, *metric;
21: PetscReal coarseRatio = PetscSqr(0.5);
22: const PetscInt *faces;
23: PetscInt *cells, *bdFaces, *bdFaceIds;
24: PetscInt dim, numCorners, cStart, cEnd, numCells, numCoarseCells, c, vStart, vEnd, numVertices, numCoarseVertices, v, numBdFaces, f, maxConeSize, size, bdSize, coff;
25: #endif
29: #ifdef PETSC_HAVE_PRAGMATIC
30: if (!dm->coarseMesh) {
31: DMGetDimension(dm, &dim);
32: DMGetCoordinateDM(dm, &coordDM);
33: DMGetDefaultSection(coordDM, &coordSection);
34: DMGetCoordinatesLocal(dm, &coordinates);
35: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
36: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
37: DMPlexUninterpolate(dm, &udm);
38: DMPlexGetMaxSizes(udm, &maxConeSize, NULL);
39: numCells = cEnd - cStart;
40: numVertices = vEnd - vStart;
41: PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);
42: VecGetArrayRead(coordinates, &coords);
43: for (v = vStart; v < vEnd; ++v) {
44: PetscInt off;
46: PetscSectionGetOffset(coordSection, v, &off);
47: x[v-vStart] = coords[off+0];
48: y[v-vStart] = coords[off+1];
49: if (dim > 2) z[v-vStart] = coords[off+2];
50: }
51: VecRestoreArrayRead(coordinates, &coords);
52: for (c = 0, coff = 0; c < numCells; ++c) {
53: const PetscInt *cone;
54: PetscInt coneSize, cl;
56: DMPlexGetConeSize(udm, c, &coneSize);
57: DMPlexGetCone(udm, c, &cone);
58: for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
59: }
60: switch (dim) {
61: case 2:
62: pragmatic_2d_init(&numVertices, &numCells, cells, x, y);
63: break;
64: case 3:
65: pragmatic_3d_init(&numVertices, &numCells, cells, x, y, z);
66: break;
67: default:
68: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic coarsening defined for dimension %d", dim);
69: }
70: /* Create boundary mesh */
71: DMLabelCreate("boundary", &bd);
72: DMPlexMarkBoundaryFaces(dm, bd);
73: DMLabelGetStratumIS(bd, 1, &bdIS);
74: DMLabelGetStratumSize(bd, 1, &numBdFaces);
75: ISGetIndices(bdIS, &faces);
76: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
77: PetscInt *closure = NULL;
78: PetscInt closureSize, cl;
80: DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);
81: for (cl = 0; cl < closureSize*2; cl += 2) {
82: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
83: }
84: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
85: }
86: PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);
87: for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
88: PetscInt *closure = NULL;
89: PetscInt closureSize, cl;
91: DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);
92: for (cl = 0; cl < closureSize*2; cl += 2) {
93: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart;
94: }
95: /* TODO Fix */
96: bdFaceIds[f] = 1;
97: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
98: }
99: ISDestroy(&bdIS);
100: DMLabelDestroy(&bd);
101: pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
102: /* Create metric */
103: size = (dim*(dim+1))/2;
104: PetscMalloc1(PetscSqr(size), &eqns);
105: MatCreateSeqDense(PETSC_COMM_SELF, size, size, eqns, &A);
106: MatCreateVecs(A, &mx, &mb);
107: VecSet(mb, 1.0);
108: for (c = 0; c < numCells; ++c) {
109: const PetscScalar *sol;
110: PetscScalar *cellCoords = NULL;
111: PetscReal e[3], vol;
112: const PetscInt *cone;
113: PetscInt coneSize, cl, i, j, d, r;
115: DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);
116: /* Only works for simplices */
117: for (i = 0, r = 0; i < dim+1; ++i) {
118: for (j = 0; j < i; ++j, ++r) {
119: for (d = 0; d < dim; ++d) e[d] = cellCoords[i*dim+d] - cellCoords[j*dim+d];
120: /* FORTRAN ORDERING */
121: if (dim == 2) {
122: eqns[0*size+r] = PetscSqr(e[0]);
123: eqns[1*size+r] = 2.0*e[0]*e[1];
124: eqns[2*size+r] = PetscSqr(e[1]);
125: } else {
126: eqns[0*size+r] = PetscSqr(e[0]);
127: eqns[1*size+r] = 2.0*e[0]*e[1];
128: eqns[2*size+r] = 2.0*e[0]*e[2];
129: eqns[3*size+r] = PetscSqr(e[1]);
130: eqns[4*size+r] = 2.0*e[1]*e[2];
131: eqns[5*size+r] = PetscSqr(e[2]);
132: }
133: }
134: }
135: MatSetUnfactored(A);
136: DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);
137: MatLUFactor(A, NULL, NULL, NULL);
138: MatSolve(A, mb, mx);
139: VecGetArrayRead(mx, &sol);
140: DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);
141: DMPlexGetCone(udm, c, &cone);
142: DMPlexGetConeSize(udm, c, &coneSize);
143: for (cl = 0; cl < coneSize; ++cl) {
144: const PetscInt v = cone[cl] - vStart;
146: if (dim == 2) {
147: metric[v*4+0] += vol*coarseRatio*sol[0];
148: metric[v*4+1] += vol*coarseRatio*sol[1];
149: metric[v*4+2] += vol*coarseRatio*sol[1];
150: metric[v*4+3] += vol*coarseRatio*sol[2];
151: } else {
152: metric[v*9+0] += vol*coarseRatio*sol[0];
153: metric[v*9+1] += vol*coarseRatio*sol[1];
154: metric[v*9+3] += vol*coarseRatio*sol[1];
155: metric[v*9+2] += vol*coarseRatio*sol[2];
156: metric[v*9+6] += vol*coarseRatio*sol[2];
157: metric[v*9+4] += vol*coarseRatio*sol[3];
158: metric[v*9+5] += vol*coarseRatio*sol[4];
159: metric[v*9+7] += vol*coarseRatio*sol[4];
160: metric[v*9+8] += vol*coarseRatio*sol[5];
161: }
162: }
163: VecRestoreArrayRead(mx, &sol);
164: }
165: for (v = 0; v < numVertices; ++v) {
166: const PetscInt *support;
167: PetscInt supportSize, s;
168: PetscReal vol, totVol = 0.0;
170: DMPlexGetSupport(udm, v+vStart, &support);
171: DMPlexGetSupportSize(udm, v+vStart, &supportSize);
172: for (s = 0; s < supportSize; ++s) {DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL); totVol += vol;}
173: for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol;
174: }
175: VecDestroy(&mx);
176: VecDestroy(&mb);
177: MatDestroy(&A);
178: DMDestroy(&udm);
179: PetscFree(eqns);
180: pragmatic_set_metric(metric);
181: pragmatic_adapt();
182: /* Read out mesh */
183: pragmatic_get_info(&numCoarseVertices, &numCoarseCells);
184: PetscMalloc1(numCoarseVertices*dim, &coarseCoords);
185: switch (dim) {
186: case 2:
187: pragmatic_get_coords_2d(x, y);
188: numCorners = 3;
189: for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*2+0] = x[v]; coarseCoords[v*2+1] = y[v];}
190: break;
191: case 3:
192: pragmatic_get_coords_3d(x, y, z);
193: numCorners = 4;
194: for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*3+0] = x[v]; coarseCoords[v*3+1] = y[v]; coarseCoords[v*3+2] = z[v];}
195: break;
196: default:
197: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic coarsening defined for dimension %d", dim);
198: }
199: pragmatic_get_elements(cells);
200: /* TODO Read out markers for boundary */
201: DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCoarseCells, numCoarseVertices, numCorners, PETSC_TRUE, cells, dim, coarseCoords, &dm->coarseMesh);
202: pragmatic_finalize();
203: PetscFree5(x, y, z, metric, cells);
204: PetscFree2(bdFaces, bdFaceIds);
205: PetscFree(coarseCoords);
206: }
207: #endif
208: PetscObjectReference((PetscObject) dm->coarseMesh);
209: *dmCoarsened = dm->coarseMesh;
210: return(0);
211: }
215: PetscErrorCode DMCoarsenHierarchy_Plex(DM dm, PetscInt nlevels, DM dmCoarsened[])
216: {
217: DM rdm = dm;
218: PetscInt c;
222: for (c = nlevels-1; c >= 0; --c) {
223: DMCoarsen(rdm, PetscObjectComm((PetscObject) dm), &dmCoarsened[c]);
224: DMCopyBoundary(rdm, dmCoarsened[c]);
225: DMSetCoarseDM(rdm, dmCoarsened[c]);
226: rdm = dmCoarsened[c];
227: }
228: return(0);
229: }