diff --git a/src/05-poisson-problem.cc b/src/05-poisson-problem.cc new file mode 100644 index 0000000000000000000000000000000000000000..23baf461bc04a1657844c90375df02b0ff362e22 --- /dev/null +++ b/src/05-poisson-problem.cc @@ -0,0 +1,338 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + + + +// included standard library headers +#include <iostream> +#include <array> + +// included dune-common headers +#include <dune/common/parallel/mpihelper.hh> +#include <dune/common/exceptions.hh> +#include <dune/common/fvector.hh> +#include <dune/common/stringutility.hh> + +// included dune-geometry headers +#include <dune/geometry/quadraturerules.hh> + +// included dune-grid headers +#include <dune/grid/io/file/vtk/vtkwriter.hh> +#include <dune/grid/utility/structuredgridfactory.hh> +#include <dune/grid/yaspgrid.hh> +#include <dune/grid/uggrid.hh> + +// included dune-istl headers +#include <dune/istl/matrix.hh> +#include <dune/istl/bcrsmatrix.hh> +#include <dune/istl/bvector.hh> +#include <dune/istl/preconditioners.hh> +#include <dune/istl/solvers.hh> + + +// included dune-localfunctions headers +#include <dune/localfunctions/lagrange/pqkfactory.hh> + +// included dune-fu-tutorial headers +#include <dune/fu-tutorial/referenceelementutility.hh> + + + +template<class Grid> +auto createCubeGrid() +{ + static const int dim = Grid::dimension; + + auto lowerLeft = Dune::FieldVector<double,dim>(0); + auto upperRight = Dune::FieldVector<double,dim>(1); + auto elementsPerDirection = std::array<unsigned int,dim>(); + for(auto& e : elementsPerDirection) + e = 2; + + return Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, elementsPerDirection); +} + +template<class Grid> +auto createSimplexGrid() +{ + static const int dim = Grid::dimension; + + auto lowerLeft = Dune::FieldVector<double,dim>(0); + auto upperRight = Dune::FieldVector<double,dim>(1); + auto elementsPerDirection = std::array<unsigned int,dim>(); + for(auto& e : elementsPerDirection) + e = 2; + + return Dune::StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, elementsPerDirection); +} + + +template<class GridView> +void writeGridView(const GridView& gridView, std::string postFix) +{ + Dune::VTKWriter<GridView> vtkWriter(gridView); + vtkWriter.write(std::string("04-gridviews-")+postFix); +} + + + +// Compute local per element matrix +template <class Element, class LocalFiniteElement, class MatrixType> +void assembleLocalStiffnessMatrix(const Element& element, const LocalFiniteElement& localFiniteElement, MatrixType& elementMatrix) +{ + auto geometryType = element.type(); + const auto& geometry = element.geometry(); + const int dim = Element::dimension; + std::size_t localSize = localFiniteElement.localBasis().size(); + + // Resize and initialize matrix + elementMatrix.setSize(localSize, localSize); + elementMatrix = 0; // fills the entire matrix with zeroes + + // Get a quadrature rule + int order = 2*(dim*localFiniteElement.localBasis().order()-1); + const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order); + + // Loop over all quadrature points + for (size_t k=0; k < quad.size(); ++k) + { + const auto& quadPos = quad[k].position(); + const auto& weight = quad[k].weight(); + + // The transposed inverse Jacobian of the map from the reference element to the element + const auto& jacobian = geometry.jacobianInverseTransposed(quadPos); + + // The multiplicative factor in the integral transformation formula + const auto& integrationElement = geometry.integrationElement(quadPos); + + // The gradients of the shape functions on the reference element + std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients; + localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients); + + // Compute the shape function gradients on the real basis functions + std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size()); + for (size_t i=0; i<gradients.size(); i++) + jacobian.mv(referenceGradients[i][0], gradients[i]); + + // Compute the actual matrix entries + for (size_t i=0; i<elementMatrix.N(); i++) + for (size_t j=0; j<elementMatrix.M(); j++ ) + elementMatrix[i][j] += ( gradients[i] * gradients[j] ) * weight * integrationElement; + } +} + + + +// Compute local per element right hand side +template <class Element, class LocalFiniteElement, class LocalRHSVector, class RHSFunction> +void assembleLocalRHS( + const Element& element, + const LocalFiniteElement& localFiniteElement, + LocalRHSVector& localRhs, + const RHSFunction& rhsFunction) +{ + auto geometryType = element.type(); + const auto& geometry = element.geometry(); + const int dim = Element::dimension; + std::size_t localSize = localFiniteElement.localBasis().size(); + + // Resize and initialize vector + localRhs.resize(localSize); + localRhs = 0; + + // Get a quadrature rule + int order = 2*(dim*localFiniteElement.localBasis().order()); + const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order); + + // Loop over all quadrature points + for (size_t k=0; k < quad.size(); ++k) + { + const auto& quadPos = quad[k].position(); + const auto& weight = quad[k].weight(); + + // Quadrature position in global world coordinates + const auto& globalQuadPos = geometry.global(quadPos); + + // The multiplicative factor in the integral transformation formula + const auto& integrationElement = geometry.integrationElement(quadPos); + + double functionValue = rhsFunction(globalQuadPos); + + // Evaluate all shape function values at this point + std::vector<Dune::FieldVector<double,1> > shapeFunctionValues; + localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues); + + // Actually compute the vector entries + for (size_t i=0; i<localSize; i++) + localRhs[i] += shapeFunctionValues[i] * functionValue * weight * integrationElement; + } +} + + +template<class GridView, class RHSFunction> +void assemblePoissonProblem( + const GridView& gridView, + Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>& matrix, + Dune::BlockVector<Dune::FieldVector<double,1>>& rhs, + const RHSFunction& rhsFunction) +{ + static const int dim = GridView::dimension; + const auto& indexSet = gridView.indexSet(); + std::size_t size = indexSet.size(dim); + + using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>; + using ElementMatrix = Dune::Matrix<Dune::FieldMatrix<double,1,1>>; + using ElementRhs = Dune::BlockVector<Dune::FieldVector<double,1>>; + + rhs.resize(size); + rhs = 0; + Dune::ImplicitMatrixBuilder<Matrix> matrixBuilder(matrix, size, size, 10, 0.1); + + Dune::PQkLocalFiniteElementCache<double, double, dim, 1> feCache; + + ElementMatrix elementMatrix; + ElementRhs elementRhs; + + for(const auto& e: Dune::elements(gridView)) + { + const auto& localFiniteElement = feCache.get(e.type()); + std::size_t localSize = localFiniteElement.localBasis().size(); + + assembleLocalStiffnessMatrix(e, localFiniteElement, elementMatrix); + assembleLocalRHS(e, localFiniteElement, elementRhs, rhsFunction); + + for(std::size_t i=0; i<localSize; ++i) + { + std::size_t globalI = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(i).subEntity(), dim); + for(std::size_t j=0; j<localSize; ++j) + { + std::size_t globalJ = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(j).subEntity(), dim); + matrixBuilder[globalI][globalJ] += elementMatrix[i][j]; + } + + rhs[globalI] = elementRhs[i]; + } + } + matrix.compress(); +} + + + +template<class GridView, class BitVector> +void computeBoundaryVertices(const GridView& gridView, BitVector& isBoundary) +{ + using namespace Dune; + using namespace Dune::FuTutorial; + + const auto& indexSet = gridView.indexSet(); + for(const auto& element: elements(gridView)) + for(const auto& intersection: intersections(gridView, element)) + if (intersection.boundary()) + for(const auto& vertex: subVertices(referenceElement(element), insideFacet(intersection))) + isBoundary[subIndex(indexSet, element, vertex)] = true; +} + + + +int main(int argc, char** argv) +{ + try{ + // Maybe initialize MPI + Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv); + + // Print process rank + if(Dune::MPIHelper::isFake) + std::cout<< "This is a sequential program." << std::endl; + else + std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size() + <<" processes!"<<std::endl; + +// using Grid = Dune::YaspGrid<2>; +// auto gridPointer = createCubeGrid<Grid>(); + +// using Grid = Dune::YaspGrid<3>; +// auto gridPointer = createCubeGrid<Grid>(); + +// using Grid = Dune::UGGrid<2>; +// auto gridPointer = createCubeGrid<Grid>(); + + using Grid = Dune::UGGrid<2>; + auto gridPointer = createSimplexGrid<Grid>(); + + Grid& grid = *gridPointer; + + grid.globalRefine(5); + for(int i=0; i<2; ++i) + { + auto gridView = grid.levelGridView(grid.maxLevel()); + for(const auto& e : Dune::elements(gridView)) + if (e.geometry().center().two_norm() < .5) + grid.mark(1, e); + grid.adapt(); + } + auto gridView = grid.leafGridView(); + using GridView = decltype(gridView); + + + using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>; + using Vector = Dune::BlockVector<Dune::FieldVector<double,1>>; + + Matrix A; + Vector rhs; + Vector x; + + + auto rhsFunction = [](auto x) { + return x.two_norm() < .5; + }; + assemblePoissonProblem(gridView, A, rhs, rhsFunction); + + x.resize(rhs.size(), 0); + + std::vector<bool> isBoundary; + isBoundary.resize(rhs.size(), false); + computeBoundaryVertices(gridView, isBoundary); + + for(std::size_t i=0; i<rhs.size(); ++i) + { + if (isBoundary[i]) + { + for(auto& entry: A[i]) + entry = 0; + A[i][i] = 1; + rhs[i] = x[i]; + } + } + + Dune::MatrixAdapter<Matrix,Vector,Vector> op(A); + Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A,1.0); + Dune::SeqSSOR<Matrix,Vector,Vector> preconditioner(A,1,1.0); + Dune::CGSolver<Vector> cg(op, + preconditioner, // preconditione + 1e-4, // desired residual reduction factor + 200, // maximum number of iterations + 2); // verbosity of the solver + + // Object storing some statistics about the solving process + Dune::InverseOperatorResult statistics; + + // Solve! + cg.apply(x, rhs, statistics); + + Dune::VTKWriter<GridView> vtkWriter(gridView); + vtkWriter.addVertexData(x, "solution"); + vtkWriter.write("05-poisson-solution"); + + return 0; + } + catch (Dune::Exception &e){ + std::cerr << "Dune reported error: " << e << std::endl; + } + catch (...){ + std::cerr << "Unknown exception thrown!" << std::endl; + } +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0b3e826c9da805c3a0f636cdf7b7eee59b5fc40d..d50032e6f78b2324df4517d536d7ca4fa9d21b1a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,3 +12,6 @@ target_link_dune_default_libraries("03-function-integration") add_executable("04-gridviews" 04-gridviews.cc) target_link_dune_default_libraries("04-gridviews") + +add_executable("05-poisson-problem" 05-poisson-problem.cc) +target_link_dune_default_libraries("05-poisson-problem")