Skip to content
Snippets Groups Projects
Select Git revision
  • ff42b7bfbc63b0777573559ea76ac3c1aaaf755e
  • master default protected
2 results

05-poisson-problem.cc

Blame
  • Carsten Gräser's avatar
    graeser authored
    ff42b7bf
    History
    Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    05-poisson-problem.cc 3.66 KiB
    // -*- 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>
    
    
    // included dune-fu-tutorial headers
    #include "04-gridviews.hh"
    #include "05-poisson-problem.hh"
    
    
    
    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<double>;
        using Vector = Dune::BlockVector<double>;
    
        Matrix A;
        Vector rhs;
        Vector x;
    
    
        auto rhsFunction = [](auto x) {
          return (x.two_norm() < .5)*10.0;
        };
        assemblePoissonProblemPQ1(gridView, A, rhs, rhsFunction);
    
        x.resize(rhs.size());
        x = 0;
    
        std::vector<bool> isBoundary;
        isBoundary.resize(rhs.size(), false);
        computeBoundaryVerticesPQ1(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::SeqILDL<Matrix,Vector,Vector> preconditioner(A,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-solution");
    
        return 0;
      }
      catch (Dune::Exception &e){
        std::cerr << "Dune reported error: " << e << std::endl;
      }
      catch (...){
        std::cerr << "Unknown exception thrown!" << std::endl;
      }
    }