Commit 90c3d73b authored by Praetorius, Simon's avatar Praetorius, Simon

improved error messages, added quadrature degree parameter

parent e4baaf0f
......@@ -3,7 +3,8 @@
#include <cassert>
#include <type_traits>
#include <dune/amdis/gridfunctions/GridFunctionConcepts.hpp>
#include <dune/amdis/GridFunctions.hpp>
#include <dune/amdis/common/Utility.hpp>
#include <dune/amdis/utility/GetDegree.hpp>
namespace AMDiS
......@@ -24,10 +25,11 @@ namespace AMDiS
* \ref ExpressionBase, and stores a copy. Additionally, it gets the
* differentiation order, to calculate the quadrature degree in \ref getDegree.
**/
GridFunctionOperatorBase(GridFunction const& gridFct, int order)
GridFunctionOperatorBase(GridFunction const& gridFct, int order, int degree = -1)
: gridFct_(gridFct)
, localFct_(localFunction(gridFct_))
, order_(order)
, degree_(degree)
{}
/// \brief Binds operator to `element` and `geometry`.
......@@ -75,7 +77,7 @@ namespace AMDiS
assert( bound_ );
int psiDegree = getPolynomialDegree(node);
int coeffDegree = order(localFct_);
int coeffDegree = getGridFctDegree(localFct_);
int degree = psiDegree + coeffDegree;
if (isSimplex_)
......@@ -100,7 +102,7 @@ namespace AMDiS
int psiDegree = getPolynomialDegree(rowNode);
int phiDegree = getPolynomialDegree(colNode);
int coeffDegree = order(localFct_);
int coeffDegree = getGridFctDegree(localFct_);
int degree = psiDegree + phiDegree + coeffDegree;
if (isSimplex_)
......@@ -111,11 +113,27 @@ namespace AMDiS
return degree;
}
protected:
template <class LocalFct,
REQUIRES(Concepts::HasOrder<LocalFct>)>
int getGridFctDegree(LocalFct const& localFct) const
{
return degree_ >= 0 ? degree_ : order(localFct_);
}
template <class LocalFct,
REQUIRES(not Concepts::HasOrder<LocalFct>)>
int getGridFctDegree(LocalFct const& localFct) const
{
return degree_;
}
private:
GridFunction gridFct_;
LocalFunction localFct_;
int order_; //< the derivative order of this operator
int degree_; //< the polynomial order of the gridFunction
bool isSimplex_ = false; //< the bound element is a simplex
bool isAffine_ = false; //< the bound geometry is affine
......@@ -166,50 +184,84 @@ namespace AMDiS
}
};
namespace Concepts
{
namespace Definition
{
struct HasGridFunctionOrder
{
template <class F, class GridView>
auto requires_(F&& f, GridView const& gridView)
-> decltype( order(localFunction(makeGridFunction(f, gridView))) );
};
}
template <class F, class GridView = Dune::YaspGrid<2>::LeafGridView>
constexpr bool HasGridFunctionOrder = models<Definition::HasGridFunctionOrder(F, GridView)>;
} // end namespace Concepts
template <class Tag, class Expr>
struct ExpressionPreOperator
{
Tag tag;
Expr expr;
int order = -1;
};
/// Store tag and expression in struct
template <class Tag, class Expr>
auto makeOperator(Tag t, Expr const& expr)
{
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 in `makeOperator()`.");
return ExpressionPreOperator<Tag, Expr>{t, expr};
}
template <class Tag, class Expr>
auto makeOperator(Tag t, Expr const& expr, int order)
{
return ExpressionPreOperator<Tag, Expr>{t, expr, order};
}
/// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
/// @{
template <class Tag, class Expr, class GridView>
auto makeGridOperator(ExpressionPreOperator<Tag, Expr> const& op, GridView const& gridView)
{
auto gridFct = makeGridFunction(op.expr, gridView, Dune::PriorityTag<42>{});
return GridFunctionOperator<Tag, decltype(gridFct)>{op.tag, gridFct};
auto gridFct = makeGridFunction(op.expr, gridView);
return GridFunctionOperator<Tag, decltype(gridFct)>{op.tag, gridFct, op.order};
}
template <class Tag, class Expr, class GridView>
auto makeGridOperator(std::reference_wrapper<ExpressionPreOperator<Tag, Expr>> op, GridView const& gridView)
{
ExpressionPreOperator<Tag, Expr> const& op_ref = op;
auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView, Dune::PriorityTag<42>{});
return GridFunctionOperator<Tag, decltype(gridFct)>{op_ref.tag, gridFct};
auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
return GridFunctionOperator<Tag, decltype(gridFct)>{op_ref.tag, gridFct, op_ref.order};
}
/// @}
/// Generate a shared_ptr to \ref GridFunctionOperator from a PreOperator (tag, expr).
/// @{
template <class Tag, class Expr, class GridView>
auto makeGridOperatorPtr(ExpressionPreOperator<Tag, Expr> const& op, GridView const& gridView)
{
auto gridFct = makeGridFunction(op.expr, gridView, Dune::PriorityTag<42>{});
return std::make_shared<GridFunctionOperator<Tag, decltype(gridFct)>>(op.tag, gridFct);
auto gridFct = makeGridFunction(op.expr, gridView);
return std::make_shared<GridFunctionOperator<Tag, decltype(gridFct)>>(op.tag, gridFct, op.order);
}
template <class Tag, class Expr, class GridView>
auto makeGridOperatorPtr(std::reference_wrapper<ExpressionPreOperator<Tag, Expr>> op, GridView const& gridView)
{
ExpressionPreOperator<Tag, Expr> const& op_ref = op;
auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView, Dune::PriorityTag<42>{});
return std::make_shared<GridFunctionOperator<Tag, decltype(gridFct)>>(op_ref.tag, gridFct);
auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
return std::make_shared<GridFunctionOperator<Tag, decltype(gridFct)>>(op_ref.tag, gridFct, op_ref.order);
}
/// @}
} // end namespace AMDiS
......@@ -6,3 +6,14 @@
#include <dune/amdis/gridfunctions/DerivativeGridFunction.hpp>
#include <dune/amdis/gridfunctions/FunctorGridFunction.hpp>
#include <dune/amdis/gridfunctions/OperationsGridFunction.hpp>
namespace AMDiS
{
// Generator for Gridfunctions from Pre-Gridfunctions
template <class PreGridFct, class GridView>
decltype(auto) makeGridFunction(PreGridFct&& preGridFct, GridView const& gridView)
{
return Impl::makeGridFunctionImpl(std::forward<PreGridFct>(preGridFct), gridView, Dune::PriorityTag<10>{});
}
} // end namespace AMDiS
......@@ -31,6 +31,7 @@
#include <dune/amdis/common/TypeDefs.hpp>
#include <dune/amdis/common/Utility.hpp>
#include <dune/amdis/GridFunctions.hpp>
#include <dune/amdis/gridfunctions/DOFVectorView.hpp>
#include <dune/amdis/io/FileWriterInterface.hpp>
......
......@@ -216,7 +216,7 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val
auto i = child(globalBasis->localView().tree(), makeTreePath(row));
auto j = child(globalBasis->localView().tree(), makeTreePath(col));
auto valueGridFct = makeGridFunction(values, globalBasis->gridView(), Dune::PriorityTag<42>{});
auto valueGridFct = makeGridFunction(values, globalBasis->gridView());
using Range = typename decltype(valueGridFct)::Range;
using BcType = DirichletBC<WorldVector,Range>;
......
......@@ -20,8 +20,8 @@ namespace AMDiS
using Transposed = GridFunctionOperator<tag::test_divtrialvec, GridFct>;
public:
GridFunctionOperator(tag::divtestvec_trial, GridFct const& expr)
: Transposed(tag::test_divtrialvec{}, expr)
GridFunctionOperator(tag::divtestvec_trial, GridFct const& expr, int degree)
: Transposed(tag::test_divtrialvec{}, expr, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -20,8 +20,8 @@ namespace AMDiS
using Transposed = GridFunctionOperator<tag::test_gradtrial, GridFct>;
public:
GridFunctionOperator(tag::gradtest_trial, GridFct const& expr)
: Transposed(tag::test_gradtrial{}, expr)
GridFunctionOperator(tag::gradtest_trial, GridFct const& expr, int degree)
: Transposed(tag::test_gradtrial{}, expr, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -20,8 +20,8 @@ namespace AMDiS
using Transposed = GridFunctionOperator<tag::testvec_gradtrial, GridFct>;
public:
GridFunctionOperator(tag::gradtest_trialvec, GridFct const& expr)
: Transposed(tag::testvec_gradtrial{}, expr)
GridFunctionOperator(tag::gradtest_trialvec, GridFct const& expr, int degree)
: Transposed(tag::testvec_gradtrial{}, expr, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -23,8 +23,8 @@ namespace AMDiS
using Transposed = GridFunctionOperator<tag::test_partialtrial, GridFct>;
public:
GridFunctionOperator(tag::partialtest_trial tag, GridFct const& expr)
: Transposed(tag::test_partialtrial{tag.comp}, expr)
GridFunctionOperator(tag::partialtest_trial tag, GridFct const& expr, int degree)
: Transposed(tag::test_partialtrial{tag.comp}, expr, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -24,8 +24,8 @@ namespace AMDiS
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
GridFunctionOperator(tag::test_divtrialvec, GridFct const& expr)
: Super(expr, 1)
GridFunctionOperator(tag::test_divtrialvec, GridFct const& expr, int degree)
: Super(expr, 1, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -24,8 +24,8 @@ namespace AMDiS
static_assert( Category::Vector<typename GridFct::Range>, "Expression must be of vector type." );
public:
GridFunctionOperator(tag::test_gradtrial, GridFct const& expr)
: Super(expr, 1)
GridFunctionOperator(tag::test_gradtrial, GridFct const& expr, int degree)
: Super(expr, 1, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -27,8 +27,8 @@ namespace AMDiS
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
GridFunctionOperator(tag::test_partialtrial tag, GridFct const& expr)
: Super(expr, 1)
GridFunctionOperator(tag::test_partialtrial tag, GridFct const& expr, int degree)
: Super(expr, 1, degree)
, comp_(tag.comp)
{}
......
......@@ -24,8 +24,8 @@ namespace AMDiS
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
GridFunctionOperator(tag::testvec_gradtrial, GridFct const& expr)
: Super(expr, 1)
GridFunctionOperator(tag::testvec_gradtrial, GridFct const& expr, int degree)
: Super(expr, 1, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -24,8 +24,8 @@ namespace AMDiS
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
GridFunctionOperator(tag::divtestvec_divtrialvec, GridFct const& expr)
: Super(expr, 2)
GridFunctionOperator(tag::divtestvec_divtrialvec, GridFct const& expr, int degree)
: Super(expr, 2, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -26,8 +26,8 @@ namespace AMDiS
"Expression must be of scalar or matrix type." );
public:
GridFunctionOperator(tag::gradtest_gradtrial, GridFct const& expr)
: Super(expr, 2)
GridFunctionOperator(tag::gradtest_gradtrial, GridFct const& expr, int degree)
: Super(expr, 2, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -28,8 +28,8 @@ namespace AMDiS
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
GridFunctionOperator(tag::partialtest_partialtrial tag, GridFct const& expr)
: Super(expr, 2)
GridFunctionOperator(tag::partialtest_partialtrial tag, GridFct const& expr, int degree)
: Super(expr, 2, degree)
, compTest_(tag.comp_test)
, compTrial_(tag.comp_trial)
{}
......
......@@ -23,8 +23,8 @@ namespace AMDiS
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
GridFunctionOperator(tag::test, GridFct const& expr)
: Super(expr, 0)
GridFunctionOperator(tag::test, GridFct const& expr, int degree)
: Super(expr, 0, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -23,8 +23,8 @@ namespace AMDiS
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
GridFunctionOperator(tag::test_trial, GridFct const& expr)
: Super(expr, 0)
GridFunctionOperator(tag::test_trial, GridFct const& expr, int degree)
: Super(expr, 0, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -23,8 +23,8 @@ namespace AMDiS
static_assert( Category::Vector<typename GridFct::Range>, "Expression must be of vector type." );
public:
GridFunctionOperator(tag::test_trialvec, GridFct const& expr)
: Super(expr, 0)
GridFunctionOperator(tag::test_trialvec, GridFct const& expr, int degree)
: Super(expr, 0, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -23,8 +23,8 @@ namespace AMDiS
static_assert( Category::Vector<typename GridFct::Range>, "Expression must be of vector type." );
public:
GridFunctionOperator(tag::testvec, GridFct const& expr)
: Super(expr, 0)
GridFunctionOperator(tag::testvec, GridFct const& expr, int degree)
: Super(expr, 0, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -20,8 +20,8 @@ namespace AMDiS
using Transposed = GridFunctionOperator<tag::test_trialvec, GridFct>;
public:
GridFunctionOperator(tag::testvec_trial, GridFct const& expr)
: Transposed(tag::test_trialvec{}, expr)
GridFunctionOperator(tag::testvec_trial, GridFct const& expr, int degree)
: Transposed(tag::test_trialvec{}, expr, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -25,8 +25,8 @@ namespace AMDiS
"Expression must be of scalar or matrix type." );
public:
GridFunctionOperator(tag::testvec_trialvec, GridFct const& expr)
: Super(expr, 0)
GridFunctionOperator(tag::testvec_trialvec, GridFct const& expr, int degree)
: Super(expr, 0, degree)
{}
template <class Context, class QuadratureRule,
......
......@@ -76,25 +76,33 @@ namespace AMDiS
/// Return the polynomial order of the function f, if a free function
/// order(f) exists, otherwise return a default value.
#ifdef AMDIS_HAS_CXX_CONSTEXPR_IF
template <class Sig, class LocalContext, class F>
template <class Sig, class LocalContext, class F,
REQUIRES(Concepts::HasOrder<F>)>
int order(AnalyticLocalFunction<Sig,LocalContext,F> const& lf)
{
if constexpr (Concepts::HasOrder<F>)
return order(lf.fct(),1);
else
return AMDIS_CALLABLE_DEFAULT_ORDER;
return order(lf.fct(),1);
}
#else
template <class Sig, class LC, class F>
int order(AnalyticLocalFunction<Sig,LC,F> const& lf)
{
return Dune::Hybrid::ifElse( bool_<Concepts::HasOrder<F>>,
[&lf](auto id) { return id( order(lf.fct(),1) ); },
[] (auto id) { return AMDIS_CALLABLE_DEFAULT_ORDER; }
);
}
#endif
// #ifdef AMDIS_HAS_CXX_CONSTEXPR_IF
// template <class Sig, class LocalContext, class F>
// int order(AnalyticLocalFunction<Sig,LocalContext,F> const& lf)
// {
// if constexpr (Concepts::HasOrder<F>)
// return order(lf.fct(),1);
// else
// return AMDIS_CALLABLE_DEFAULT_ORDER;
// }
// #else
// template <class Sig, class LC, class F>
// int order(AnalyticLocalFunction<Sig,LC,F> const& lf)
// {
// return Dune::Hybrid::ifElse( bool_<Concepts::HasOrder<F>>,
// [&lf](auto id) { return id( order(lf.fct(),1) ); },
// [] (auto id) { return AMDIS_CALLABLE_DEFAULT_ORDER; }
// );
// }
// #endif
/// A Gridfunction that evaluates a function with global coordinates
......@@ -180,22 +188,24 @@ namespace AMDiS
return AnalyticPreGridFunction<Function>{f};
}
/** @} **/
/// Generator function for \ref AnalyticGridFunction expressions from callable functions.
template <class Function, class GridView,
REQUIRES(Concepts::CallableDomain<Function>)>
auto makeGridFunction(Function const& f, GridView const& gridView, Dune::PriorityTag<3>)
{
return AnalyticGridFunction<Function, GridView>{f, gridView};
}
/// Generator function for \ref AnalyticGridFunction expressions from \ref AnalyticPreGridFunction.
template <class Function, class GridView>
auto makeGridFunction(AnalyticPreGridFunction<Function> const& pre, GridView const& gridView, Dune::PriorityTag<3>)
namespace Impl
{
return AnalyticGridFunction<Function, GridView>{pre.fct_, gridView};
}
/// Generator function for \ref AnalyticGridFunction expressions from callable functions.
template <class Function, class GridView,
REQUIRES(Concepts::CallableDomain<Function>)>
auto makeGridFunctionImpl(Function const& f, GridView const& gridView, Dune::PriorityTag<3>)
{
return AnalyticGridFunction<Function, GridView>{f, gridView};
}
/** @} **/
/// Generator function for \ref AnalyticGridFunction expressions from \ref AnalyticPreGridFunction.
template <class Function, class GridView>
auto makeGridFunctionImpl(AnalyticPreGridFunction<Function> const& pre, GridView const& gridView, Dune::PriorityTag<3>)
{
return AnalyticGridFunction<Function, GridView>{pre.fct_, gridView};
}
} // end namespace Impl
} // end namespace AMDiS
......@@ -57,21 +57,19 @@ namespace AMDiS
return 0;
}
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};
}
private:
T value_;
};
template <class R, class D, class LC, class T>
auto derivative(ConstantLocalFunction<R(D),LC,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),LC,DerivativeRange>{diff};
}
/// A Gridfunction that evaluates a function with global coordinates
template <class T, class GridView>
class ConstantGridFunction
......@@ -157,20 +155,15 @@ namespace AMDiS
} // end namespace Concepts
/**
* \addtogroup GridFunctions
* @{
**/
/// Generator for a GridFunction representing a constant value
template <class T, class GridView,
REQUIRES(Concepts::ConstantToGridFunction<T>)>
auto makeGridFunction(T const& value, GridView const& gridView, Dune::PriorityTag<2>)
namespace Impl
{
return ConstantGridFunction<T,GridView>{value, gridView};
}
/** @} **/
/// Generator for a GridFunction representing a constant value
template <class T, class GridView,
REQUIRES(Concepts::ConstantToGridFunction<T>)>
auto makeGridFunctionImpl(T const& value, GridView const& gridView, Dune::PriorityTag<2>)
{
return ConstantGridFunction<T,GridView>{value, gridView};
}
} // end namespace Impl
} // end namespace AMDiS
......@@ -16,6 +16,17 @@ namespace AMDiS
/// A functor that evaluates to the global coordinates
struct CoordsFunction
{
template <class T, int N>
Dune::FieldVector<T, N> const& operator()(Dune::FieldVector<T, N> const& x) const
{
return x;
}
friend int order(CoordsFunction const& /*f*/, int /*d*/)
{
return 1;
}
struct Derivative
{
template <class T, int N>
......@@ -24,28 +35,32 @@ namespace AMDiS
return Dune::DiagonalMatrix<T,N>{T(1)};
}
};
template <class T, int N>
Dune::FieldVector<T, N> const& operator()(Dune::FieldVector<T, N> const& x) const
{
return x;
}
friend Derivative partial(CoordsFunction const& /*f*/, index_t<0>)
{
return Derivative{};
}
friend int order(CoordsFunction const& /*f*/, int /*d*/)
{
return 1;
}
};
/// A functor that evaluates to a component of the global coordinates
struct CoordsCompFunction
{
/// Constructor. Stores the component `comp` of the coordinates.
explicit CoordsCompFunction(int comp)
: comp_(comp)
{}
template <class T, int N>
T const& operator()(Dune::FieldVector<T, N> const& x) const
{
return x[comp_];
}
friend int order(CoordsCompFunction const& /*f*/, int /*d*/)
{
return 1;
}
struct Derivative
{
explicit Derivative(int comp)
......@@ -65,27 +80,11 @@ namespace AMDiS
int comp_;
};
/// Constructor. Stores the component `comp` of the coordinates.
explicit CoordsCompFunction(int comp)
: comp_(comp)
{}
template <class T, int N>
T const& operator()(Dune::FieldVector<T, N> const& x) const
{
return x[comp_];
}
friend Derivative partial(CoordsCompFunction const& f, index_t<0>)
{
return Derivative{f.comp_};
}
friend int order(CoordsCompFunction const& /*f*/, int /*d*/)
{
return 1;
}
private:
int comp_;
};
......@@ -117,25 +116,21 @@ namespace AMDiS
: std::true_type {};
}
/**
* \addtogroup GridFunctions
* @{
**/
/// Generator for GridFunction from the pre-GridFunction \ref CoordsFunction
template <class GridView>
auto makeGridFunction(Operation::CoordsFunction const& f, GridView const& gridView, Dune::PriorityTag<1>)
{
return makeAnalyticGridFunction(f, gridView);
}
/// Generator for GridFunction from the pre-GridFunction \ref CoordsCompFunction
template <class GridView>
auto makeGridFunction(Operation::CoordsCompFunction const& f, GridView const& gridView, Dune::PriorityTag<1>)
namespace Impl
{
return makeAnalyticGridFunction(f, gridView);
}
/// Generator for GridFunction from the pre-GridFunction \ref CoordsFunction
template <class GridView>
auto makeGridFunctionImpl(Operation::CoordsFunction const& f, GridView const& gridView, Dune::PriorityTag<1>)
{
return makeAnalyticGridFunction(f, gridView);
}
/** @} **/
/// Generator for GridFunction from the pre-GridFunction \ref CoordsCompFunction
template <class GridView>
auto makeGridFunctionImpl(Operation::CoordsCompFunction const& f, GridView const& gridView, Dune::PriorityTag<1>)
{
return makeAnalyticGridFunction(f, gridView);
}
} // end namespace Impl
} // end namespace AMDiS
......@@ -201,7 +201,6 @@ namespace AMDiS
public:
/// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
DOFVectorView(DOFVector<GlobalBasis> const& dofVector, TreePath const& treePath)
: dofVector_(&dofVector)
......@@ -228,9 +227,7 @@ namespace AMDiS
return entitySet_;
}
public: