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 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, numCoarseBdFaces, f, maxConeSize, size, bdSize, coff; 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 #ifdef PETSC_HAVE_PRAGMATIC 30 if (!mesh->coarseMesh) { 31 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 32 ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr); 33 ierr = DMGetDefaultSection(coordDM, &coordSection);CHKERRQ(ierr); 34 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 35 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 36 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 37 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 38 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); 39 numCells = cEnd - cStart; 40 numVertices = vEnd - vStart; 41 ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr); 42 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 43 for (v = vStart; v < vEnd; ++v) { 44 PetscInt off; 45 46 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 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 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 52 for (c = 0, coff = 0; c < numCells; ++c) { 53 const PetscInt *cone; 54 PetscInt coneSize, cl; 55 56 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 57 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 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 ierr = DMLabelCreate("boundary", &bd);CHKERRQ(ierr); 72 ierr = DMPlexMarkBoundaryFaces(dm, bd);CHKERRQ(ierr); 73 ierr = DMLabelGetStratumIS(bd, 1, &bdIS);CHKERRQ(ierr); 74 ierr = DMLabelGetStratumSize(bd, 1, &numBdFaces);CHKERRQ(ierr); 75 ierr = ISGetIndices(bdIS, &faces);CHKERRQ(ierr); 76 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 77 PetscInt *closure = NULL; 78 PetscInt closureSize, cl; 79 80 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 81 for (cl = 0; cl < closureSize*2; cl += 2) { 82 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 83 } 84 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 85 } 86 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 87 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 88 PetscInt *closure = NULL; 89 PetscInt closureSize, cl; 90 91 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 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 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 98 } 99 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 100 ierr = DMLabelDestroy(&bd);CHKERRQ(ierr); 101 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 102 /* Create metric */ 103 size = (dim*(dim+1))/2; 104 ierr = PetscMalloc1(PetscSqr(size), &eqns);CHKERRQ(ierr); 105 ierr = MatCreateSeqDense(PETSC_COMM_SELF, size, size, eqns, &A);CHKERRQ(ierr); 106 ierr = MatCreateVecs(A, &mx, &mb);CHKERRQ(ierr); 107 ierr = VecSet(mb, 1.0);CHKERRQ(ierr); 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; 114 115 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr); 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 ierr = MatSetUnfactored(A);CHKERRQ(ierr); 136 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr); 137 ierr = MatLUFactor(A, NULL, NULL, NULL);CHKERRQ(ierr); 138 ierr = MatSolve(A, mb, mx);CHKERRQ(ierr); 139 ierr = VecGetArrayRead(mx, &sol);CHKERRQ(ierr); 140 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr); 141 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 142 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 143 for (cl = 0; cl < coneSize; ++cl) { 144 const PetscInt v = cone[cl] - vStart; 145 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 ierr = VecRestoreArrayRead(mx, &sol);CHKERRQ(ierr); 164 } 165 for (v = 0; v < numVertices; ++v) { 166 const PetscInt *support; 167 PetscInt supportSize, s; 168 PetscReal vol, totVol = 0.0; 169 170 ierr = DMPlexGetSupport(udm, v+vStart, &support);CHKERRQ(ierr); 171 ierr = DMPlexGetSupportSize(udm, v+vStart, &supportSize);CHKERRQ(ierr); 172 for (s = 0; s < supportSize; ++s) {ierr = DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL);CHKERRQ(ierr); totVol += vol;} 173 for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol; 174 } 175 ierr = VecDestroy(&mx);CHKERRQ(ierr); 176 ierr = VecDestroy(&mb);CHKERRQ(ierr); 177 ierr = MatDestroy(&A);CHKERRQ(ierr); 178 ierr = DMDestroy(&udm);CHKERRQ(ierr); 179 ierr = PetscFree(eqns);CHKERRQ(ierr); 180 pragmatic_set_metric(metric); 181 pragmatic_adapt(); 182 /* Read out mesh */ 183 pragmatic_get_info(&numCoarseVertices, &numCoarseCells); 184 ierr = PetscMalloc1(numCoarseVertices*dim, &coarseCoords);CHKERRQ(ierr); 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 ierr = DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCoarseCells, numCoarseVertices, numCorners, PETSC_TRUE, cells, dim, coarseCoords, &mesh->coarseMesh);CHKERRQ(ierr); 202 pragmatic_finalize(); 203 ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr); 204 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 205 ierr = PetscFree(coarseCoords);CHKERRQ(ierr); 206 } 207 #endif 208 ierr = PetscObjectReference((PetscObject) mesh->coarseMesh);CHKERRQ(ierr); 209 *dmCoarsened = mesh->coarseMesh; 210 PetscFunctionReturn(0); 211 } 212