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