Select Git revision
build_kernel.sh
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;
}
}