diff --git a/src/07-poisson-problem-functions.cc b/src/07-poisson-problem-functions.cc new file mode 100644 index 0000000000000000000000000000000000000000..8039bea5ecbfc77112df1bf626b3b03d953b4c6f --- /dev/null +++ b/src/07-poisson-problem-functions.cc @@ -0,0 +1,258 @@ +// -*- 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/io/file/vtk/subsamplingvtkwriter.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-functions headers +#include <dune/functions/functionspacebases/pqknodalbasis.hh> +#include <dune/functions/functionspacebases/interpolate.hh> +#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh> + +// included dune-fu-tutorial headers +#include <dune/fu-tutorial/referenceelementutility.hh> + + +// included dune-fu-tutorial headers +#include "04-gridviews.hh" +#include "05-poisson-problem.hh" + + +template<class GridView, class RHSFunction, class Basis> +void assemblePoissonProblem( + const GridView& gridView, + Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>& matrix, + Dune::BlockVector<Dune::FieldVector<double,1>>& rhs, + const RHSFunction& rhsFunction, + const Basis& basis) +{ + std::size_t size = basis.size(); + + 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, 100, 0.1); + + ElementMatrix elementMatrix; + ElementRhs elementRhs; + + auto localView = basis.localView(); + auto localIndexSet = basis.localIndexSet(); + + for(const auto& e: Dune::elements(gridView)) + { + localView.bind(e); + localIndexSet.bind(localView); + + const auto& localFiniteElement = localView.tree().finiteElement(); + + 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 = localIndexSet.index(i); + for(std::size_t j=0; j<localSize; ++j) + { + std::size_t globalJ = localIndexSet.index(j); + matrixBuilder[globalI][globalJ] += elementMatrix[i][j]; + } + + rhs[globalI] += elementRhs[i]; + } + } + matrix.compress(); +} + +template<class GridView, class BitVector, class Basis> +void computeBoundaryVertices(const GridView& gridView, BitVector& isBoundary, const Basis& basis) +{ + using namespace Dune; + using namespace Dune::FuTutorial; + + auto localView = basis.localView(); + auto localIndexSet = basis.localIndexSet(); + + for(const auto& element: elements(gridView)) + { + localView.bind(element); + localIndexSet.bind(localView); + + const auto& localFiniteElement = localView.tree().finiteElement(); + std::size_t localSize = localFiniteElement.localBasis().size(); + + for(const auto& intersection: intersections(gridView, element)) + { + if (intersection.boundary()) + { + for(std::size_t i=0; i<localSize; ++i) + { + auto localKey = localFiniteElement.localCoefficients().localKey(i); + if (localKey.codim() > 0) + for(const auto& subEntity: subEntities(referenceElement(element), insideFacet(intersection), codim(localKey.codim()))) + if (subEntity == localKey.subEntity()) + isBoundary[localIndexSet.index(i)] = 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(1); + 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)*10.0; + }; + + + using Basis = typename Dune::Functions::PQkNodalBasis<GridView,4>; + Basis basis(gridView); + + assemblePoissonProblem(gridView, A, rhs, rhsFunction, basis); + + x.resize(basis.size(), 0); + + auto dirichletFunction = [](auto x) { + return sin(x[0] * 2.0*M_PI)*.1; + }; +// Dune::Functions::interpolate(basis, Dune::TypeTree::hybridTreePath(), dirichletNodes, dirichletFunction); + Dune::Functions::interpolate(basis, x, dirichletFunction); + + std::vector<bool> isBoundary; + isBoundary.resize(basis.size(), false); + computeBoundaryVertices(gridView, isBoundary, basis); + + 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); + + auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(basis, x); + + { + Dune::VTKWriter<GridView> vtkWriter(gridView); + vtkWriter.addVertexData(xFunction, Dune::VTK::FieldInfo("solution", Dune::VTK::FieldInfo::Type::scalar, 1)); + vtkWriter.write("07-poisson-solution-functions"); + } { + Dune::SubsamplingVTKWriter<GridView> vtkWriter(gridView,4); + vtkWriter.addVertexData(xFunction, Dune::VTK::FieldInfo("solution", Dune::VTK::FieldInfo::Type::scalar, 1)); + vtkWriter.write("07-poisson-solution-functions-subsamples"); + } + + + 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 067c003e6a4aec0120a46dd8f0145d27de950c21..51a4e0cbce7dbf96765f369ffe9f38e1b5d79083 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,3 +18,6 @@ target_link_dune_default_libraries("05-poisson-problem") add_executable("06-interpolation" 06-interpolation.cc) target_link_dune_default_libraries("06-interpolation") + +add_executable("07-poisson-problem-functions" 07-poisson-problem-functions.cc) +target_link_dune_default_libraries("07-poisson-problem-functions")