Commit c9ad5da6 authored by Praetorius, Simon's avatar Praetorius, Simon

merged origin/develop

parents 5d3b30c0 076a7e14
......@@ -7,7 +7,6 @@ if(NOT (dune-common_DIR OR dune-common_ROOT OR
${PROJECT_BINARY_DIR})
endif()
#find dune-common and set the module path
find_package(dune-common REQUIRED)
list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"
......
......@@ -8,4 +8,4 @@ Version: 0.1
Maintainer: simon.praetorius@tu-dresden.de
#depending on
Depends: dune-common dune-geometry dune-localfunctions dune-typetree dune-grid dune-functions
Suggests: dune-uggrid dune-alugrid dune-foamgrid
#Suggests: dune-uggrid dune-alugrid dune-foamgrid
......@@ -60,14 +60,22 @@ namespace AMDiS
bool asmMatrix, bool asmVector) const;
template <class ElementContainer, class Container, class Operators, class Geometry, class... Bases>
template <class ElementContainer, class Container, class Operators, class Geometry, class Basis>
void assembleElementOperators(
ElementContainer& elementContainer,
Container& container,
Operators& operators,
Geometry const& geometry,
Bases const&... subBases) const;
Basis const& subBasis) const;
template <class ElementContainer, class Container, class Operators, class Geometry, class RowBasis, class ColBasis>
void assembleElementOperators(
ElementContainer& elementContainer,
Container& container,
Operators& operators,
Geometry const& geometry,
RowBasis const& rowBasis,
ColBasis const& colBasis) const;
/// Finish insertion into the matrix and assembles boundary conditions
/// Return the number of nonzeros assembled into the matrix
......
......@@ -94,23 +94,63 @@ void Assembler<Traits>::assemble(
template <class Traits>
template <class ElementContainer, class Container, class Operators, class Geometry, class... LocalViews>
template <class ElementContainer, class Container, class Operators, class Geometry, class LocalView>
void Assembler<Traits>::assembleElementOperators(
ElementContainer& elementContainer,
Container& container,
Operators& operators,
Geometry const& geometry,
LocalViews const&... localViews) const
LocalView const& localView) const
{
auto const& element = getElement(localViews...);
auto const& gridView = getGridView(localViews...);
auto const& element = getElement(localView);
auto const& gridView = getGridView(localView);
bool add = false;
auto assemble_operators = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list) {
scaled.op->bind(element, geometry);
bool add_op = scaled.op->assemble(context, elementContainer, localViews.tree()...);
bool add_op = scaled.op->assemble(context, elementContainer, localView.tree());
scaled.op->unbind();
add = add || add_op;
}
};
// assemble element operators
assemble_operators(element, operators.element);
// assemble intersection operators
if (!operators.intersection.empty()
|| (!operators.boundary.empty() && element.hasBoundaryIntersections()))
{
for (auto const& intersection : intersections(gridView, element)) {
if (intersection.boundary())
assemble_operators(intersection, operators.boundary);
else
assemble_operators(intersection, operators.intersection);
}
}
}
template <class Traits>
template <class ElementContainer, class Container, class Operators, class Geometry, class RowLocalView, class ColLocalView>
void Assembler<Traits>::assembleElementOperators(
ElementContainer& elementContainer,
Container& container,
Operators& operators,
Geometry const& geometry,
RowLocalView const& rowLocalView, ColLocalView const& colLocalView) const
{
auto const& element = getElement(rowLocalView, colLocalView);
auto const& gridView = getGridView(rowLocalView, colLocalView);
bool add = false;
auto assemble_operators = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list) {
scaled.op->bind(element, geometry);
bool add_op = scaled.op->assemble(context, elementContainer, rowLocalView.tree(), colLocalView.tree());
scaled.op->unbind();
add = add || add_op;
}
......
......@@ -35,12 +35,13 @@ install(FILES
LocalAssembler.hpp
LocalAssemblerBase.hpp
Mesh.hpp
Operations.hpp
Operators.hpp
Output.hpp
ProblemInstat.hpp
ProblemInstat.inc.hpp
ProblemInstatBase.hpp
ProblemInterationInterface.hpp
ProblemIterationInterface.hpp
ProblemStat.hpp
ProblemStat.inc.hpp
ProblemStatBase.hpp
......@@ -55,4 +56,4 @@ add_subdirectory("gridfunctions")
add_subdirectory("io")
add_subdirectory("linear_algebra")
add_subdirectory("operations")
add_subdirectory("utility")
\ No newline at end of file
add_subdirectory("utility")
......@@ -278,7 +278,7 @@ namespace AMDiS
{
using RawExpr = Underlying_t<Expr>;
static_assert(Concepts::HasGridFunctionOrder<RawExpr>,
"Polynomial degree of expression can not be deduced. You need to provide an explicit value for polynomial order or a quadrature rule in `makeOperator()`.");
"Polynomial degree of expression can not be deduced. You need to provide an explicit value for polynomial order or a quadrature rule in 'makeOperator()'.");
return ExpressionPreOperator<Tag, Expr, tag::deduce>{tag, expr};
}
......
......@@ -137,6 +137,11 @@ void ProblemStat<Traits>::addMatrixOperator(
RowTreePath row, ColTreePath col,
double* factor, double* estFactor)
{
static_assert( Concepts::PreTreePath<RowTreePath>,
"row must be a valid treepath, or an integer/index-constant");
static_assert( Concepts::PreTreePath<ColTreePath>,
"col must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(row));
auto j = child(globalBasis->localView().tree(), makeTreePath(col));
......@@ -156,6 +161,11 @@ void ProblemStat<Traits>::addMatrixOperator(
RowTreePath row, ColTreePath col,
double* factor, double* estFactor)
{
static_assert( Concepts::PreTreePath<RowTreePath>,
"row must be a valid treepath, or an integer/index-constant");
static_assert( Concepts::PreTreePath<ColTreePath>,
"col must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(row));
auto j = child(globalBasis->localView().tree(), makeTreePath(col));
using Intersection = typename GridView::Intersection;
......@@ -175,6 +185,9 @@ void ProblemStat<Traits>::addVectorOperator(
TreePath path,
double* factor, double* estFactor)
{
static_assert( Concepts::PreTreePath<TreePath>,
"path must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(path));
auto op = makeGridOperator(preOp, globalBasis->gridView());
......@@ -193,6 +206,9 @@ void ProblemStat<Traits>::addVectorOperator(
TreePath path,
double* factor, double* estFactor)
{
static_assert( Concepts::PreTreePath<TreePath>,
"path must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(path));
using Intersection = typename GridView::Intersection;
......
......@@ -18,4 +18,4 @@ install(FILES
ZeroOrderTestvec.hpp
ZeroOrderTestvecTrial.hpp
ZeroOrderTestvecTrialvec.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/operators)
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/assembler)
#install headers
install(FILES
Assigner.hpp
ClonablePtr.hpp
Concepts.hpp
ConceptsBase.hpp
......@@ -10,7 +11,6 @@ install(FILES
Loops.hpp
Math.hpp
Mpl.hpp
Operations.hpp
OptionalDelete.hpp
ScalarTypes.hpp
Timer.hpp
......
......@@ -341,4 +341,15 @@ namespace AMDiS
return C;
}
template <class T, int M, int N>
Dune::FieldMatrix<T,M,N>& multiplies_ABt(Dune::FieldMatrix<T, M, N> const& A, Dune::DiagonalMatrix<T, N> const& B, Dune::FieldMatrix<T,M,N>& C)
{
for (int m = 0; m < M; ++m) {
for (int n = 0; n < N; ++n) {
C[m][n] = A[m][n] * B.diagonal(n);
}
}
return C;
}
} // end namespace AMDiS
......@@ -15,14 +15,16 @@ namespace AMDiS
return unsigned(c) - unsigned('0');
}
template <std::size_t N>
constexpr std::size_t string2num(const char (&arr)[N])
template <char... digits>
constexpr std::size_t string2num()
{
const char arr[] = {digits...};
assert(arr[0] != '-' && "Negative integral constant");
std::size_t result = 0;
std::size_t power = 1;
const int N = sizeof...(digits);
for (std::size_t i = 0; i < N; ++i) {
char c = arr[N - 1 - i];
result+= char2digit(c) * power;
......@@ -38,7 +40,7 @@ namespace AMDiS
template <char... digits>
constexpr auto operator"" _c()
{
return index_<Impl::string2num<sizeof...(digits)>({digits...})>;
return index_<Impl::string2num<digits...>()>;
}
} // end namespace AMDiS
......@@ -45,18 +45,6 @@ namespace AMDiS
return fct_(geometry_.value().global(local));
}
friend auto derivative(AnalyticLocalFunction const& lf)
{
static_assert(Concepts::HasPartial<Function>,
"No partial(_0,...) defined for Functor F of AnalyticLocalFunction.");
auto df = partial(lf.fct(), index_<0>);
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
using DerivativeSignature = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range(D);
return AnalyticLocalFunction<DerivativeSignature,LocalContext,decltype(df)>{df};
}
Function const& fct() const
{
return fct_;
......@@ -77,6 +65,20 @@ namespace AMDiS
}
template <class R, class D, class LocalContext, class F>
auto derivative(AnalyticLocalFunction<R(D),LocalContext,F> const& lf)
{
static_assert(Concepts::HasPartial<F>,
"No partial(_0,...) defined for Functor F of AnalyticLocalFunction.");
auto df = partial(lf.fct(), index_<0>);
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
using DerivativeSignature = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range(D);
return AnalyticLocalFunction<DerivativeSignature,LocalContext,decltype(df)>{df};
}
/// \class AnalyticGridFunction
/// \brief A Gridfunction that evaluates a function with global coordinates.
/**
......@@ -116,21 +118,11 @@ namespace AMDiS
}
/// \brief Return the LocalFunction of the AnalyticGridFunction.
friend auto localFunction(AnalyticGridFunction const& gf)
friend LocalFunction localFunction(AnalyticGridFunction const& gf)
{
return LocalFunction{gf.fct_};
}
/// \brief Return a GridFunction representing the derivative of a functor.
// [expects: Functor f has free function derivative(f)]
friend auto derivative(AnalyticGridFunction const& gf)
{
static_assert(Concepts::HasPartial<Function>,
"No partial(_0,...) defined for Functor of AnalyticLocalFunction.");
auto df = partial(gf.fct(), index_<0>);
return AnalyticGridFunction<decltype(df), GridView>{df, gf.entitySet().gridView()};
}
EntitySet const& entitySet() const
{
......@@ -146,6 +138,19 @@ namespace AMDiS
};
/// \brief Return a GridFunction representing the derivative of a functor.
// [expects: Functor f has free function derivative(f)]
template <class F, class GV>
auto derivative(AnalyticGridFunction<F,GV> const& gf)
{
static_assert(Concepts::HasPartial<F>,
"No partial(_0,...) defined for Functor of AnalyticLocalFunction.");
auto df = partial(gf.fct(), index_<0>);
return AnalyticGridFunction<decltype(df), GV>{df, gf.entitySet().gridView()};
}
#ifndef DOXYGEN
// A pre-GridFunction that just stores the function
template <class Function>
......
......@@ -47,19 +47,23 @@ namespace AMDiS
return 0;
}
/// \relates ConstantLocalFunction
friend auto derivative(ConstantLocalFunction const& lf)
{
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range;
DerivativeRange diff(0);
return ConstantLocalFunction<DerivativeRange(D),LocalContext,DerivativeRange>{diff};
}
template <class R_, class D_, class L_, class T_>
friend auto derivative(ConstantLocalFunction<R_(D_), L_, T_> const& lf);
private:
T value_;
};
/// \relates ConstantLocalFunction
template <class R, class D, class LocalContext, class T>
auto derivative(ConstantLocalFunction<R(D), LocalContext, T> const& lf)
{
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range;
DerivativeRange diff(0);
return ConstantLocalFunction<DerivativeRange(D),LocalContext,DerivativeRange>{diff};
}
/**
* \addtogroup GridFunctions
......@@ -101,7 +105,7 @@ namespace AMDiS
}
/// Return the LocalFunction of the \ref ConstantGridFunction. \relates ConstantGridFunction
friend auto localFunction(ConstantGridFunction const& gf)
friend LocalFunction localFunction(ConstantGridFunction const& gf)
{
return LocalFunction{gf.value_};
}
......
......@@ -179,9 +179,10 @@ namespace AMDiS
using RawRange = std::decay_t<Range>;
using LocalDomain = typename EntitySet::LocalCoordinate;
using LocalFunction = FunctorLocalFunction<RawRange(LocalDomain), Functor, LocalFct<GridFunctions>...>;
public:
using LocalFunction = FunctorLocalFunction<RawRange(LocalDomain), Functor, LocalFct<GridFunctions>...>;
/// Constructor. Stores copies of the functor and gridfunctions.
template <class... GridFcts>
explicit FunctorGridFunction(Functor const& fct, GridFcts&&... gridFcts)
......@@ -196,22 +197,15 @@ namespace AMDiS
MakeSeq_t<sizeof...(GridFunctions)>{});
}
/// \brief Creates a LocalFunction from the LocalFunctions of the GridFunctions. \relates FunctorLocalFunction
friend LocalFunction localFunction(FunctorGridFunction const& gf)
{
return Dune::Std::apply([&gf](auto const&... gridFcts)
{
return LocalFunction{gf.fct_, localFunction(gridFcts)...};
},
gf.gridFcts_);
}
/// Return the stored \ref EntitySet of the first GridFunction
EntitySet const& entitySet() const
{
return std::get<0>(gridFcts_).entitySet();
}
auto const& fct() const { return fct_; }
auto const& gridFcts() const { return gridFcts_; }
private:
template <class Outer, class Inner, std::size_t... I>
auto eval(Outer outer, Inner inner, Indices<I...>) const
......@@ -225,6 +219,19 @@ namespace AMDiS
};
/// \brief Creates a LocalFunction from the LocalFunctions of the GridFunctions. \relates FunctorLocalFunction
template <class F, class... GFs>
auto localFunction(FunctorGridFunction<F,GFs...> const& gf)
{
return Dune::Std::apply([&gf](auto const&... gridFcts)
{
using LocalFunction = typename FunctorGridFunction<F,GFs...>::LocalFunction;
return LocalFunction{gf.fct(), localFunction(gridFcts)...};
},
gf.gridFcts());
}
#ifndef DOXYGEN
// Generator function for FunctorGridFunction expressions
template <class Functor, class... GridFcts>
......
......@@ -2,28 +2,11 @@
install(FILES
HierarchicWrapper.hpp
LinearAlgebraBse.hpp
LinearAlgebraBase.hpp
LinearSolverInterface.hpp
PreconditionerInterface.hpp
RunnerInterface.hpp
SolverInfo.hpp
mtl/BITL_Preconditioner.hpp
mtl/BITL_Solver.hpp
mtl/BlockMTLMatrix.hpp
mtl/BlockMTLVector.hpp
mtl/Copy.hpp
mtl/DOFMatrix.hpp
mtl/DOFVector.hpp
mtl/ITL_Preconditioner.hpp
mtl/ITL_Solver.hpp
mtl/KrylovRunner.hpp
mtl/LinearSolver.hpp
mtl/Mapper.hpp
mtl/MTLDenseVector.hpp
mtl/Preconditioner.hpp
mtl/SystemMatrix.hpp
mtl/SystemVector.hpp
mtl/UmfpackRunner.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/linear_algebra)
add_subdirectory("mtl")
\ No newline at end of file
add_subdirectory("mtl")
......@@ -22,3 +22,5 @@ install(FILES
SystemVector.hpp
UmfpackRunner.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/linear_algebra/mtl)
add_subdirectory("itl")
......@@ -8,7 +8,7 @@
#include <boost/numeric/itl/pc/ic_0.hpp>
#include <boost/numeric/mtl/vector/assigner.hpp>
#include "itl/masslumping.hpp"
#include <dune/amdis/linear_algebra/mtl/itl/masslumping.hpp>
#include <dune/amdis/CreatorMap.hpp>
#include <dune/amdis/linear_algebra/mtl/Preconditioner.hpp>
......
......@@ -18,13 +18,13 @@
// #include <boost/numeric/itl/krylov/tfqmr.hpp>
// more solvers defined in AMDiS
#include "itl/minres.hpp"
#include "itl/gcr.hpp"
#include "itl/fgmres.hpp"
#include "itl/fgmres_householder.hpp"
#include "itl/gmres2.hpp"
#include "itl/gmres_householder.hpp"
#include "itl/preonly.hpp"
#include <dune/amdis/linear_algebra/mtl/itl/minres.hpp>
#include <dune/amdis/linear_algebra/mtl/itl/gcr.hpp>
#include <dune/amdis/linear_algebra/mtl/itl/fgmres.hpp>
#include <dune/amdis/linear_algebra/mtl/itl/fgmres_householder.hpp>
#include <dune/amdis/linear_algebra/mtl/itl/gmres2.hpp>
#include <dune/amdis/linear_algebra/mtl/itl/gmres_householder.hpp>
#include <dune/amdis/linear_algebra/mtl/itl/preonly.hpp>
#include <dune/amdis/CreatorMap.hpp>
#include <dune/amdis/Initfile.hpp>
......
install(FILES
block_diagonal.hpp
details.hpp
fgmres_householder.hpp
fgmres.hpp
gcr.hpp
gmres2.hpp
gmres_householder.hpp
hypre.hpp
masslumping.hpp
minres.hpp
preonly.hpp
umfpack2_solve.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/linear_algebra/mtl/itl)
......@@ -10,6 +10,68 @@
namespace AMDiS
{
namespace Concepts
{
/** \addtogroup Concepts
* @{
**/
namespace Definition
{
template <class TP>
struct IsPreTreePath
: any_of_t<std::is_same<TP, int>::value, std::is_same<TP,std::size_t>::value>
{};
template <int I>
struct IsPreTreePath<int_t<I>>
: std::true_type
{};
template <std::size_t I>
struct IsPreTreePath<index_t<I>>
: std::true_type
{};
template <std::size_t... I>
struct IsPreTreePath<Indices<I...>>
: std::true_type
{};
} // end namespace Definition
template <class TP>
constexpr bool PreTreePath = Dune::TypeTree::IsTreePath<TP>::value || Definition::IsPreTreePath<TP>::value;
/** @} **/
} // end namespace Concepts
#ifdef DOXYGEN
/// \brief Converts a (pre)TreePath into a HybridTreePath
/**
* Converts an integer, an integralconstant or a Dune TreePath into an
* \ref Dune::TypeTree::HybridTreePath that is used in GlobalBasis traversal.
*
* **Requirements:**
* - `PreTreePath` one of
* + integer type (`int, std::size_t`),
* + integral constant (`std::integral_constant<[int|std::size_t], i>`)
* + any Dune TreePath (`TreePath<std::size_t...>, DynamicTreePath, HybridTreePath<class... T>`)
*
* **Example:**
* ```
* makeTreePath(1),
* makeTreePath(int_<2>),
* makeTreePath(treepath(1, int_<2>))
* ```
**/
template <class PreTreePath>
auto makeTreePath(PreTreePath tp);
#else // DOXYGEN
auto makeTreePath(int i) { return Dune::TypeTree::hybridTreePath(std::size_t(i)); }
auto makeTreePath(std::size_t i) { return Dune::TypeTree::hybridTreePath(i); }
......@@ -30,11 +92,12 @@ namespace AMDiS
template <class TP>
auto makeTreePath(TP const&) {
static_assert( Dune::TypeTree::IsTreePath<TP>::value,
static_assert( Concepts::PreTreePath<TP>,
"Argument must be a valid treepath, or an integer/index-constant");
return Dune::TypeTree::hybridTreePath();
}
<<<<<<< HEAD
template <class... T>
auto makeDynamicTreePath(Dune::TypeTree::HybridTreePath<T...> const& tp)
{
......@@ -44,57 +107,43 @@ namespace AMDiS
});
return tmp.view();
}
=======
#endif // DOXYGEN
>>>>>>> 076a7e1479a879871f8599eade9f911c41ca09a8
namespace Impl
{
template <std::size_t I, std::size_t N>
struct PrintTreePath
{
template <class... T>
static void apply(std::ostream& os, Dune::TypeTree::HybridTreePath<T...> const& tp)
{
os << "," << treePathIndex(tp, index_<I>);
PrintTreePath<I+1,N>::apply(os, tp);
}
};
template <std::size_t N>
struct PrintTreePath<N,N>
{
template <class... T>
static void apply(std::ostream&, Dune::TypeTree::HybridTreePath<T...> const&) {}
};