Commit 87233e85 authored by Praetorius, Simon's avatar Praetorius, Simon

Initial commit

Pipeline #2483 failed with stages
in 26 seconds
"files.associations": {
"string_view": "cpp",
"*.tcc": "cpp",
"cctype": "cpp",
"clocale": "cpp",
"cmath": "cpp",
"cstdarg": "cpp",
"cstdint": "cpp",
"cstdio": "cpp",
"cstdlib": "cpp",
"cwchar": "cpp",
"cwctype": "cpp",
"exception": "cpp",
"initializer_list": "cpp",
"iosfwd": "cpp",
"iostream": "cpp",
"istream": "cpp",
"limits": "cpp",
"new": "cpp",
"ostream": "cpp",
"stdexcept": "cpp",
"streambuf": "cpp",
"system_error": "cpp",
"type_traits": "cpp",
"typeinfo": "cpp"
\ No newline at end of file
# We require version CMake version 3.1 to prevent issues
# with dune_enable_all_packages and older CMake versions.
cmake_minimum_required(VERSION 3.1)
project(amg_error CXX)
if(NOT (dune-common_DIR OR dune-common_ROOT OR
"${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
string(REPLACE ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR
#find dune-common and set the module path
find_package(dune-common REQUIRED)
#include the dune macros
# start a dune project with information from dune.module
# finalize the dune project, e.g. generating config.h etc.
Preparing the Sources
Additional to the software mentioned in README you'll need the
following programs installed on your system:
cmake >= 2.8.12
Getting started
If these preliminaries are met, you should run
dunecontrol all
which will find all installed dune modules as well as all dune modules
(not installed) which sources reside in a subdirectory of the current
directory. Note that if dune is not installed properly you will either
have to add the directory where the dunecontrol script resides (probably
./dune-common/bin) to your path or specify the relative path of the script.
Most probably you'll have to provide additional information to dunecontrol
(e. g. compilers, configure options) and/or make options.
The most convenient way is to use options files in this case. The files
define four variables:
CMAKE_FLAGS flags passed to cmake (during configure)
An example options file might look like this:
#use this options to configure and make if no other options are given
-DCMAKE_CXX_FLAGS='-Wall -pedantic' \
-DCMAKE_INSTALL_PREFIX=/install/path" #Force g++-5 and set compiler flags
If you save this information into example.opts you can pass the opts file to
dunecontrol via the --opts option, e. g.
dunecontrol --opts=example.opts all
More info
dunecontrol --help
for further options.
The full build system is described in the dune-common/doc/buildsystem (Git version) or under share/doc/dune-common/buildsystem if you installed DUNE!
Version: @VERSION@
Description: amg_error module
Requires: dune-common dune-geometry dune-localfunctions dune-typetree dune-grid dune-istl dune-functions
Libs: -L${libdir}
Cflags: -I${includedir}
# File for module specific CMake tests.
set(modules "Amg_errorMacros.cmake")
/* begin amg_error
put the definitions for config.h specific to
your project here. Everything above will be
/* begin private */
/* Name of package */
/* Define to the address where bug reports for this package should be sent. */
/* Define to the full name of this package. */
/* Define to the full name and version of this package. */
/* Define to the one symbol short name of this package. */
/* Define to the home page for this package. */
/* Define to the version of this package. */
/* end private */
/* Define to the version of amg_error */
/* Define to the major version of amg_error */
/* Define to the minor version of amg_error */
/* Define to the revision of amg_error */
/* end amg_error
Everything below here will be overwritten
# shortcut for creating the and Doxyfile
# This file contains local changes to the doxygen configuration
# please use '+=' to add files/directories to the lists
# The INPUT tag can be used to specify the files and/or directories that contain
# documented source files. You may enter file names like "myfile.cpp" or
# directories like "/usr/src/myproject". Separate the files or directories
# with spaces.
INPUT += @top_srcdir@/dune/
# see e.g. dune-grid for the examples of mainpage and modules
# INPUT += @srcdir@/mainpage \
# @srcdir@/modules
# The EXCLUDE tag can be used to specify files and/or directories that should
# excluded from the INPUT source files. This way you can easily exclude a
# subdirectory from a directory tree whose root is specified with the INPUT tag.
# EXCLUDE += @top_srcdir@/dune/amg_error/test
# The EXAMPLE_PATH tag can be used to specify one or more files or
# directories that contain example code fragments that are included (see
# the \include command).
# EXAMPLE_PATH += @top_srcdir@/src
# The IMAGE_PATH tag can be used to specify one or more files or
# directories that contain image that are included in the documentation (see
# the \image command).
# IMAGE_PATH += @top_srcdir@/dune/amg_error/pics
# Dune module information file #
#Name of the module
Module: amg_error
Version: 0.1
#depending on
Depends: dune-common dune-geometry dune-localfunctions dune-typetree dune-grid dune-istl dune-functions
#pragma once
#include <tuple>
#include <type_traits>
#include <utility>
namespace AMDiS
namespace Tools
namespace Impl_
template <class Functor, class Tuple, std::size_t... I>
constexpr decltype(auto) apply_impl(Functor&& f, Tuple&& t, std::index_sequence<I...>)
using std::get;
return f(get<I>(std::forward<Tuple>(t))...);
template <class Functor, std::size_t I0, std::size_t... I>
constexpr decltype(auto) apply_indices_impl(Functor&& f, std::integral_constant<std::size_t, I0>, std::index_sequence<I...>)
return f(std::integral_constant<std::size_t, I0+I>{}...);
} // namespace Impl_
template <class Functor, class Tuple>
constexpr decltype(auto) apply(Functor&& f, Tuple&& t)
return Impl_::apply_impl(std::forward<Functor>(f), std::forward<Tuple>(t),
template <class Functor, class... Args>
constexpr decltype(auto) apply_variadic(Functor&& f, Args&&... args)
return apply(std::forward<Functor>(f), std::forward_as_tuple(args...));
template <std::size_t N, class Functor>
constexpr decltype(auto) apply_indices(Functor&& f)
return Impl_::apply_indices_impl(std::forward<Functor>(f), std::integral_constant<std::size_t, 0>{},
template <class Functor, std::size_t N>
constexpr decltype(auto) apply_indices(Functor&& f, std::integral_constant<std::size_t, N>)
return Impl_::apply_indices_impl(std::forward<Functor>(f), std::integral_constant<std::size_t, 0>{},
template <class Functor, std::size_t I0, std::size_t I1>
constexpr decltype(auto) apply_indices(Functor&& f, std::integral_constant<std::size_t, I0>, std::integral_constant<std::size_t, I1>)
return Impl_::apply_indices_impl(std::forward<Functor>(f), std::integral_constant<std::size_t, I0>{},
} // end namespace AMDiS
#install headers
#pragma once
#include <initializer_list>
#include <tuple>
#include <type_traits>
#include <utility>
namespace AMDiS
namespace Tools
namespace Impl_
template <class T>
void ignored_evaluation(std::initializer_list<T>&&) { /* do nothing */ }
template <class Functor, class... Args>
constexpr void for_variadic(Functor&& f, Args&&... args)
Impl_::ignored_evaluation<int>({0, (f(std::forward<Args>(args)), 0)...});
template <std::size_t... I, class Tuple, class Functor>
constexpr void for_each(std::index_sequence<I...>, Tuple&& tuple, Functor&& f)
using std::get;
Impl_::ignored_evaluation<int>({0, (f(get<I>(tuple)), 0)...});
template <class Tuple, class Functor>
constexpr void for_each(Tuple&& tuple, Functor&& f)
Tools::for_each(std::make_index_sequence<std::tuple_size<std::remove_reference_t<Tuple>>::value>{}, std::forward<Tuple>(tuple), std::forward<Functor>(f));
template <std::size_t I0 = 0, std::size_t... I, class Functor>
constexpr void for_range(std::index_sequence<I...>, Functor&& f)
Impl_::ignored_evaluation<int>({0, (f(std::integral_constant<std::size_t, I0+I>{}), 0)...});
template <std::size_t I0, std::size_t I1, class Functor>
constexpr void for_range(std::integral_constant<std::size_t, I0> i0, std::integral_constant<std::size_t, I1> i1, Functor&& f)
Tools::for_range<I0>(std::make_index_sequence<std::size_t(I1-I0)>{}, std::forward<Functor>(f));
template <std::size_t N, class Functor>
constexpr void for_range(std::integral_constant<std::size_t, N>, Functor&& f)
Tools::for_range(std::make_index_sequence<N>{}, std::forward<Functor>(f));
template <std::size_t I0, std::size_t I1, class Functor>
constexpr void for_range(Functor&& f)
Tools::for_range<I0>(std::make_index_sequence<std::size_t(I1-I0)>{}, std::forward<Functor>(f));
} // end namespace Tools
} // end namespace AMDiS
This diff is collapsed.
#pragma once
#include <dune/common/version.hh>
namespace AMDiS
// dune version independent extraction of node type from preBasis
template <class PB, class TP>
using Node_t =
typename PB::template Node<TP>;
typename PB::Node;
// dune version independent extraction of node indexset type from preBasis
template <class PB, class TP>
using NodeIndexSet_t =
typename PB::template IndexSet<TP>;
typename PB::IndexSet;
// dune version independent creation of node from preBasis
template <class PB, class TP>
auto makeNode(PB const& preBasis, TP const& treePath)
return preBasis.node(treePath);
return preBasis.makeNode();
} // end namespace AMDiS
#pragma once
#include <tuple>
#include <type_traits>
#include <utility>
#include <dune/common/hash.hh>
namespace AMDiS
namespace Impl
// Recursive template code derived from Matthieu M.
template <class Tuple, std::size_t I = std::tuple_size<Tuple>::value - 1>
struct HashTupleImpl
static void apply(size_t& seed, Tuple const& tuple)
HashTupleImpl<Tuple, I-1>::apply(seed, tuple);
Dune::hash_combine(seed, std::get<I>(tuple));
template <class Tuple>
struct HashTupleImpl<Tuple, 0>
static void apply(std::size_t& seed, Tuple const& tuple)
Dune::hash_combine(seed, std::get<0>(tuple));
} // end namespace Impl
template <class Tuple, template <class> class Map>
struct MapTuple;
template <class Tuple, template <class> class Map>
using MapTuple_t = typename MapTuple<Tuple,Map>::type;
template <class... T, template <class> class Map>
struct MapTuple<std::tuple<T...>, Map>
using type = std::tuple<Map<T>...>;
template <class Indices, template <std::size_t> class Map>
struct IndexMapTuple;
template <class Indices, template <std::size_t> class Map>
using IndexMapTuple_t = typename IndexMapTuple<Indices,Map>::type;
template <std::size_t... I, template <std::size_t> class Map>
struct IndexMapTuple<std::index_sequence<I...>, Map>
using type = std::tuple<Map<I>...>;
} // end namespace AMDiS
#pragma once
#include <cassert>
#include <set>
#include <dune/common/unused.hh>
#include <dune/grid/common/datahandleif.hh>
namespace AMDiS
/// \brief Determine for each border entity which processor owns it
* All entities must be uniquely owned by exactly one processor, but they can
* exist on multiple processors. For interior, overlap, and ghost entities the
* assignment is trivial: interior: owner, otherwise not owner.
* For border entities (codim != 0) the ownership is not known a priori and must
* be communicated. Here we assign the entity to the processor with the lowest rank.
template <class Grid>
class UniqueBorderPartitionDataHandle
: public Dune::CommDataHandleIF<UniqueBorderPartitionDataHandle<Grid>, int>
using IdSet = typename Grid::GlobalIdSet;
using IdType = typename IdSet::IdType;
using EntitySet = std::set<IdType>;
/// \brief Construct a UniqueBorderPartition DataHandle to be used in a GridView
/// communicator.
* \param rank The own processor rank
* \param borderEntities A set of entity ids filled with all border entities
* owned by this processor after communication.
* \param idSet The id set of entity ids to store in borderEntities,
* typically the grid globalIdSet.
* NOTE: Since idSet is stored by reference it must not go out of scope
* until all calls to \ref gather and \ref scatter are finished.
UniqueBorderPartitionDataHandle(int rank, EntitySet& borderEntities, IdSet const& idSet)
: myrank_(rank)
, borderEntities_(&borderEntities)
, idSet_(idSet)
// Communicate all entities except for grid elements
bool contains(int /*dim*/, int codim) const { return codim != 0; }
// communicate exactly one integer, the rank
bool fixedSize(int /*dim*/, int /*codim*/) const { return true; }
// Always contains one int, the rank
template <class Entity>
std::size_t size(Entity const& e) const { return 1; }
template <class MessageBuffer, class Entity>
void gather(MessageBuffer& buff, Entity const& e) const
// insert all border entities
template <class MessageBuffer, class Entity>
void scatter(MessageBuffer& buff, Entity const& e, std::size_t n)
assert(n == 1);
int rank = 0;;
// remove all border entities with rank < myrank_
if (rank < myrank_)
int myrank_;
EntitySet* borderEntities_;
IdSet const& idSet_;
} // end namespace AMDiS
#ifndef AMG_ERROR_HH
#define AMG_ERROR_HH
// add your classes here
#endif // AMG_ERROR_HH
# include "config.h"
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/amg_error/GlobalBasisIdSet.hpp>
#include "assemble.hh"
#include "communication.hh"
// load function
double f(Dune::FieldVector<double,2> const& x)
double r2 =;
double ux = std::exp(-10.0 * r2);
return -(400.0 * r2 - 20.0 * x.size()) * ux;
// boundary function
double g(Dune::FieldVector<double,2> const& x)
return std::exp(-10.0 *;
int main(int argc, char** argv)
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
using Grid = Dune::YaspGrid<2>;
auto grid = Grid({1.0, 1.0}, {8,8}, 0, 1); // 8x8 domain with overlap=1
using namespace Dune::Functions::BasisFactory;
auto basis = makeBasis(grid.leafGridView(), lagrange<2>());
using Basis = decltype(basis);
// assemble the linear system
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double,1>>;
Matrix mat;
buildMatrixPattern(basis, mat);
Vector rhs(basis.dimension());
assemble(basis, mat, rhs, f);
Vector x(rhs);
fillBoundaryCondition(basis, mat, x, rhs, g);
// setup the communication object
using DofIdType = typename AMDiS::GlobalBasisIdSet<Basis>::IdType;
using Comm = Dune::OwnerOverlapCopyCommunication<DofIdType, typename Basis::size_type>;
Comm comm(Dune::MPIHelper::getCommunicator(), Dune::SolverCategory::Category::overlapping);
buildCommunication(basis, comm);
// setup the linear solver and preconditioner
using LinOp = Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,Comm>;
LinOp linOp(mat, comm);
using ScalarProduct = Dune::OverlappingSchwarzScalarProduct<Vector,Comm>;
ScalarProduct sp(comm);
using Smoother = Dune::SeqSSOR<Matrix, Vector, Vector>;
using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
SmootherArgs smootherArgs;
using Norm = Dune::Amg::FirstDiagonal;
using SymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Matrix,Norm>>;
using UnSymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Matrix,Norm>>;
Dune::Amg::Parameters parameter;
SymCriterion criterion(parameter);
using Precon = Dune::Amg::AMG<LinOp, Vector, Smoother, Comm>;
Precon precon(linOp, criterion, smootherArgs, comm);
using Solver = Dune::CGSolver<Vector>;
Solver solver(linOp, sp, precon, 1.e-8, 1000, 2, false);
// solver the linear system
Dune::InverseOperatorResult statistics;
solver.apply(x, rhs, statistics);
#pragma once
#include <type_traits>
#include <vector>
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/istl/matrixindexset.hh>
// fill element matrix and element vector
template <class LocalView, class ElementMatrix, class ElementVector, class F>
void assembleLocal(LocalView const& localView, ElementMatrix& elementMatrix, ElementVector& elementVector, F const& f)
auto const& element = localView.element();
auto const& node = localView.tree();
auto const& localFE = node.finiteElement();
std::size_t size = localFE.size();
// quadrature rule on the element
using QuadratureRules = Dune::QuadratureRules<double, 2>;
auto quad = QuadratureRules::rule(element.type(), 4);
auto geo = element.geometry();
auto const& localBasis = localFE.localBasis();
using LocalBasisType = std::decay_t<decltype(localBasis)>;
using JacobianType = typename LocalBasisType::Traits::JacobianType;
using RangeType = typename LocalBasisType::Traits::RangeType;
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
const auto& local = quad[iq].position();
const double dx = geo.integrationElement(local) * quad[iq].weight();
const double coeff = f(;
thread_local std::vector<RangeType> shapeValues;
localBasis.evaluateFunction(local, shapeValues);
thread_local std::vector<JacobianType> shapeGradients;
localBasis.evaluateJacobian(local, shapeGradients);
const auto jacobian = geo.jacobianInverseTransposed(local);
thread_local std::vector<Dune::FieldVector<double,2>> gradients;
for (std::size_t i = 0; i < gradients.size(); ++i)[i][0], gradients[i]);
for (std::size_t i = 0; i < size; ++i) {
const auto local_i = node.localIndex(i);
elementMatrix[local_i][local_i] += dx * gradients[i].two_norm2();
elementVector[local_i] += dx * coeff * shapeValues[i];
for (std::size_t j = i+1; j < size; ++j) {
const auto local_j = node.localIndex(j);
const auto value = dx * gradients[i].dot(gradients[j]);