Actual source code: plexcoarsen.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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: }