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