Commit 00b09ad2 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

allow to pass quadrature rule to makeOperator, cleanup of documentation

parent 90c3d73b
......@@ -10,25 +10,30 @@ INTERNAL_DOCS = NO
MARKDOWN_SUPPORT = YES
EXCLUDE_SYMBOLS = AMDiS::Impl \
AMDiS::traits::Impl \
AMDiS::Math::Impl_ \
AMDiS::Concepts::Impl_ \
AMDiS::detail \
itl::details
PREDEFINED += HAVE_UMFPACK \
HAVE_ALBERTA \
HAVE_UG \
AMDIS_BACKEND_MTL
# 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/amdis \
@top_srcdir@/dune/amdis/assembler \
@top_srcdir@/dune/amdis/common \
@top_srcdir@/dune/amdis/utility \
@top_srcdir@/dune/amdis/gridfunctions \
@top_srcdir@/dune/amdis/io \
@top_srcdir@/dune/amdis/linear_algebra \
@top_srcdir@/dune/amdis/linear_algebra/mtl
@top_srcdir@/dune/amdis/linear_algebra/mtl \
@top_srcdir@/dune/amdis/operations \
@top_srcdir@/dune/amdis/utility
# see e.g. dune-grid for the examples of mainpage and modules
#INPUT += @srcdir@/mainpage \
# @srcdir@/modules
......@@ -50,4 +55,4 @@ EXAMPLE_PATH += @top_srcdir@/src
# directories that contain image that are included in the documentation (see
# the \image command).
# IMAGE_PATH += @top_srcdir@/dune/amdis/pics
# IMAGE_PATH += @top_srcdir@/doc/pics
......@@ -30,7 +30,7 @@ namespace AMDiS
* Must be implemented by sub classes of CreatorInterface.
* Creates a new instance of the sub class of BaseClass.
*/
virtual shared_ptr<BaseClass> create() = 0;
virtual std::shared_ptr<BaseClass> create() = 0;
};
/**
......
......@@ -9,7 +9,31 @@
namespace AMDiS
{
template <class GridFunction>
/**
* \addtogroup operators
* @{
**/
/// \brief The main implementation of an operator to be used in a \ref LocalAssembler.
/**
* An Operator that takes a GridFunction as coefficient.
* Provides quadrature rules and handles the evaluation of the GridFunction at
* local coordinates.
*
* The class is specialized, by deriving from it, in \ref GridFunctionOperator.
*
* \tparam GridFunction The GridFunction, a LocalFunction is created from, and
* that is evaluated at quadrature points.
* \tparam QuadratureCreator A functor that provides a \ref Dune::QuadratureRule.
*
* **Requirements:**
* - `GridFunction` models the \ref Concepts::GridFunction
* - `QuadratureCreator` models \ref Concepts::Callable<Dune::GeometryType, LocalFunction, F>
* where `F` is a functor of the signature `int(int)` that calculates the
* degree of the (bi)linear-form. The argument passed to `F` is the polynomial
* order of the GridFunction.
**/
template <class GridFunction, class QuadratureCreator>
class GridFunctionOperatorBase
{
using LocalFunction = decltype(localFunction(std::declval<GridFunction>()));
......@@ -25,11 +49,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, int degree = -1)
GridFunctionOperatorBase(GridFunction const& gridFct, QuadratureCreator creator, int order)
: gridFct_(gridFct)
, localFct_(localFunction(gridFct_))
, quadCreator_(creator)
, order_(order)
, degree_(degree)
{}
/// \brief Binds operator to `element` and `geometry`.
......@@ -57,7 +81,7 @@ namespace AMDiS
localFct_.unbind();
}
/// Return expressions iq'th value. Must be initialized in \ref init before.
/// Return expression value at LocalCoordinates
auto coefficient(LocalCoordinate const& local) const
{
assert( bound_ );
......@@ -65,6 +89,16 @@ namespace AMDiS
}
/// Create a quadrature rule using the \ref QuadratureCreator by calculating the
/// quadrature order necessary to integrate the (bi)linear-form accurately.
template <class... Nodes>
decltype(auto) getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
{
auto degreeFunctor = [&,this](int order) -> int { return this->getDegree(order, nodes...); };
return quadCreator_(type, localFct_, degreeFunctor);
}
protected:
/// \brief Return the quadrature degree for a vector operator.
/**
* The quadrature degree that is necessary, to integrate the expression
......@@ -72,12 +106,11 @@ namespace AMDiS
* the order of derivatives, this operator implements.
**/
template <class Node>
int getDegree(Node const& node) const
int getDegree(int coeffDegree, Node const& node) const
{
assert( bound_ );
int psiDegree = getPolynomialDegree(node);
int coeffDegree = getGridFctDegree(localFct_);
int degree = psiDegree + coeffDegree;
if (isSimplex_)
......@@ -96,13 +129,12 @@ namespace AMDiS
* the order of derivatives, this operator implements.
**/
template <class RowNode, class ColNode>
int getDegree(RowNode const& rowNode, ColNode const& colNode) const
int getDegree(int coeffDegree, RowNode const& rowNode, ColNode const& colNode) const
{
assert( bound_ );
int psiDegree = getPolynomialDegree(rowNode);
int phiDegree = getPolynomialDegree(colNode);
int coeffDegree = getGridFctDegree(localFct_);
int degree = psiDegree + phiDegree + coeffDegree;
if (isSimplex_)
......@@ -113,27 +145,12 @@ 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
QuadratureCreator quadCreator_; //< a creator to provide a quadrature rule
int order_; //< the derivative order of this operator
bool isSimplex_ = false; //< the bound element is a simplex
bool isAffine_ = false; //< the bound geometry is affine
......@@ -147,13 +164,13 @@ namespace AMDiS
* \ref calculateElementMatrix, if it is a vector or matrix operator,
* respectively.
**/
template <class Tag, class GridFct>
template <class Tag, class GridFct, class QuadratureCreator>
class GridFunctionOperator
: public GridFunctionOperatorBase<GridFct>
: public GridFunctionOperatorBase<GridFct, QuadratureCreator>
{
public:
GridFunctionOperator(Tag, GridFct const& gridFct)
: GridFunctionOperatorBase<GridFct>(gridFct, 0)
GridFunctionOperator(Tag, GridFct const& gridFct, QuadratureCreator const& quadCreator)
: GridFunctionOperatorBase<GridFct, QuadratureCreator>(gridFct, 0, quadCreator)
{}
/// Assemble a local element vector on the element that is bound.
......@@ -184,6 +201,8 @@ namespace AMDiS
}
};
#ifndef DOXYGEN
namespace Concepts
{
namespace Definition
......@@ -202,66 +221,157 @@ namespace AMDiS
} // end namespace Concepts
namespace tag
{
struct deduce {};
struct order {};
struct rule {};
} // end namespace tag
template <class Tag, class Expr, class QuadTag, class Rule = void>
struct ExpressionPreOperator;
template <class Tag, class Expr>
struct ExpressionPreOperator
struct ExpressionPreOperator<Tag, Expr, tag::deduce>
{
Tag tag;
Expr expr;
};
template <class Tag, class Expr>
struct ExpressionPreOperator<Tag, Expr, tag::order>
{
Tag tag;
Expr expr;
int order;
Dune::QuadratureType::Enum qt;
};
template <class Tag, class Expr, class Rule>
struct ExpressionPreOperator<Tag, Expr, tag::rule, Rule>
{
Tag tag;
Expr expr;
int order = -1;
Rule const& rule;
};
#endif
/// Store tag and expression in struct
/// Store tag and expression to create a \ref GridFunctionOperator
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()`.");
"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>{t, expr};
return ExpressionPreOperator<Tag, Expr, tag::deduce>{t, expr};
}
/// Store tag and expression and polynomial order of expression to create a \ref GridFunctionOperator
template <class Tag, class Expr>
auto makeOperator(Tag t, Expr const& expr, int order)
auto makeOperator(Tag t, Expr const& expr, int order,
Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
{
return ExpressionPreOperator<Tag, Expr, tag::order>{t, expr, order, qt};
}
/// Store tag and expression and a quadrature rule to create a \ref GridFunctionOperator
template <class Tag, class Expr, class ctype, int dim>
auto makeOperator(Tag t, Expr const& expr, Dune::QuadratureRule<ctype,dim> const& rule)
{
return ExpressionPreOperator<Tag, Expr>{t, expr, order};
return ExpressionPreOperator<Tag, Expr, tag::rule, Dune::QuadratureRule<ctype,dim>>{t, expr, rule};
}
/** @} **/
#ifndef DOXYGEN
namespace Impl
{
// Deduce polynomial order of expression automatically. Create standard quadrature rule with degree
// build up of operator derivative order, and polynomial degrees of trial/test functions.
template <class Tag, class Expr, class GridView>
auto makeQuadCreator(ExpressionPreOperator<Tag, Expr, tag::deduce> const& /*op*/, GridView const& /*gridView*/)
{
using QuadratureRules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
return [](auto type, auto&& localFct, auto getDegree) -> auto const&
{
return QuadratureRules::rule(type, getDegree(order(localFct)));
};
}
// Provide polynomial order of expression explicitly
template <class Tag, class Expr, class GridView>
auto makeQuadCreator(ExpressionPreOperator<Tag, Expr, tag::order> const& op, GridView const& /*gridView*/)
{
using QuadratureRules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
return [order=op.order,qt=op.qt](auto type, auto&& /*localFct*/, auto getDegree) -> auto const&
{
return QuadratureRules::rule(type, getDegree(order), qt);
};
}
// Provide quadrature rule explicitly
template <class Tag, class Expr, class Rule, class GridView>
auto makeQuadCreator(ExpressionPreOperator<Tag, Expr, tag::rule, Rule> const& op, GridView const& /*gridView*/)
{
return [&rule=op.rule](auto&&...) -> auto const&
{
return rule;
};
}
} // end namespace Impl
/// 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)
template <class Tag, class... Args, class GridView>
auto makeGridOperator(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
{
auto gridFct = makeGridFunction(op.expr, gridView);
return GridFunctionOperator<Tag, decltype(gridFct)>{op.tag, gridFct, op.order};
auto quadCreator = Impl::makeQuadCreator(op, gridView);
using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
return GridFctOp{op.tag, gridFct, quadCreator};
}
template <class Tag, class Expr, class GridView>
auto makeGridOperator(std::reference_wrapper<ExpressionPreOperator<Tag, Expr>> op, GridView const& gridView)
template <class Tag, class... Args, class GridView>
auto makeGridOperator(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
{
ExpressionPreOperator<Tag, Expr> const& op_ref = op;
ExpressionPreOperator<Tag, Args...> const& op_ref = op;
auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
return GridFunctionOperator<Tag, decltype(gridFct)>{op_ref.tag, gridFct, op_ref.order};
auto quadCreator = Impl::makeQuadCreator(op_ref, gridView);
using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
return GridFctOp{op_ref.tag, gridFct, quadCreator};
}
/// @}
/// 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)
template <class Tag, class... Args, class GridView>
auto makeGridOperatorPtr(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
{
auto gridFct = makeGridFunction(op.expr, gridView);
return std::make_shared<GridFunctionOperator<Tag, decltype(gridFct)>>(op.tag, gridFct, op.order);
auto quadCreator = Impl::makeQuadCreator(op, gridView);
using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
return std::make_shared<GridFctOp>(op.tag, gridFct, quadCreator);
}
template <class Tag, class Expr, class GridView>
auto makeGridOperatorPtr(std::reference_wrapper<ExpressionPreOperator<Tag, Expr>> op, GridView const& gridView)
template <class Tag, class... Args, class GridView>
auto makeGridOperatorPtr(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
{
ExpressionPreOperator<Tag, Expr> const& op_ref = op;
ExpressionPreOperator<Tag, Args...> const& op_ref = op;
auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
return std::make_shared<GridFunctionOperator<Tag, decltype(gridFct)>>(op_ref.tag, gridFct, op_ref.order);
auto quadCreator = Impl::makeQuadCreator(op_ref, gridView);
using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
return std::make_shared<GridFctOp>(op_ref.tag, gridFct, quadCreator);
}
/// @}
#endif // DOXYGEN
} // end namespace AMDiS
#pragma once
/**
* \defgroup GridFunctions GridFunction module
* \brief Defines GridFunctions to be used in operators, boundary-conditions,
* interpolation and integration.
*
* GridFunctions are expressions build up of some elementary terms and can be
* used to construct a \ref GridFunctionOperator, can be interpolated to a
* \ref DOFVector, and can be integrated over a GridView.
*
* Thus, GridFunctions are an important incredient to formulate the bilinear
* and linear forms und to postprocess the solutions.
*
* **Examples:**
* 1. Usage of GridFunctions to build Operators:
* ```
* ProblemStat<Traits> prob("name");
* prob.initialize(INIT_ALL);
*
* auto opB = makeOperator(BiLinearForm, Expression);
* prob.addMatrixOperator(opB, Row, Col);
*
* auto opL = makeOperator(LinearForm, Expression);
* prob.addVectorOperator(opL, Row);
* ```
*
* 2. Usage of GridFunctions in BoundaryConditions:
* ```
* prob.addDirichletBC(Predicate, Row, Col, Expression);
* ```
*
* 3. Interpolate a GridFunction to a DOFVector:
* ```
* prob.getSolution(_0).interpol(Expression);
* ```
*
* 4. Integrate a GridFunction on a GridView:
* ```
* auto value = integrate(Expression, prob.leafGridView());
* ```
*
* **Remarks:**
* - An `Expression` is anything, a GridFunction can be created from, sometimes
* also called PreGridFunction. It includes constants, functors callable with
* GlobalCoordinates, and any combination of GridFunctions.
* - Anything that needs a quadrature formula, e.g., makeOperator() and
* integrate(), needs to determine the (approximative) polynomial degree of
* the GridFunctions. If the Gridfunction builds a polynomial expression, it
* can be deduced automatically, i.e. if it includes constants, DOFVectors,
* and arithmetic operator operator+, operator-, or operator*.
*
* If the polynomial order can not be deduced, the compiler gives an error.
* Then, these functions, accept an additional argument, to provide either the
* polynomial degree of the expression, or a quadrature rule explicitly.
*
* *Examples:*
* + `auto op1 = makeOperator(B, 1.0 + pow<2>(prob.getSolution(_0)));`
* + `auto op2 = makeOperator(B, sin(X(0)), 4);`
* + `auto op3 = makeOperator(B, sin(X(0)), Dune::QuadratureRules(Dune::GeometryType::simplex, 4));`
* + `auto value1 = integrate(sin(X(0)), 4);`
**/
#include <dune/amdis/gridfunctions/AnalyticGridFunction.hpp>
#include <dune/amdis/gridfunctions/ConstantGridFunction.hpp>
#include <dune/amdis/gridfunctions/CoordsGridFunction.hpp>
......@@ -9,7 +70,40 @@
namespace AMDiS
{
// Generator for Gridfunctions from Pre-Gridfunctions
/// \brief Generator for Gridfunctions from Expressions (PreGridfunctions)
/**
* \ingroup GridFunctions
* Create an evaluable GridFunction from an expression that itself can not be
* evaluated. Therefore, it binds the GridFunction to a GridView.
*
* **Example:**
* ```
* ProblemStat<Traits> prob("name");
* prob.initialize(INIT_ALL);
*
* auto gridFct = makeGridFunction(Expression, prob.leafGridView());
*
* // eval GridFunction at GlobalCoordinates
* auto value = gridFct(Dune::FieldVector<double,2>{1.0, 2.0});
*
* auto localFct = localFunction(gridFct);
* for (auto const& element : elements(prob.leafGridView())) {
* localFct.bind(element);
* // eval LocalFunction at local coordinates
* auto x = localFct(element.geometry().center());
* localFct.unbind();
* }
* ```
*
* In contrast to Expressions, GridFunctions can be evaluated, and
* - have the free-function \ref localFunction() to obtain a LocalFunction
* - its LocalFunctions have the free-function \ref order() to obtain the
* polynomial order of the Expression (if available)
* - its LocalFunctions have the free-function \ref derivative() to
* differentiate the Expression with respect to global Coordinates.
* A derivative Expression can be created, using \ref gradientAtQP() that
* can be converted to a GridFunction afterwards.
**/
template <class PreGridFct, class GridView>
decltype(auto) makeGridFunction(PreGridFct&& preGridFct, GridView const& gridView)
{
......
......@@ -83,14 +83,11 @@ namespace AMDiS
ElementMatrixVector& elementMatrixVector,
Nodes const&... nodes) final
{
decltype(auto) localGeometry = getLocalGeometry(localContext);
using QuadratureRules = Dune::QuadratureRules<typename Geometry::ctype, LocalContext::mydimension>;
int degree = op_.getDegree(nodes...);
auto const& quad = QuadratureRules::rule(localGeometry.type(), degree);
auto&& localGeometry = getLocalGeometry(localContext);
auto&& quad = op_.getQuadratureRule(localGeometry.type(), nodes...);
Context data{localContext, getGeometry(), localGeometry};
assembleImpl(data, quad, elementMatrixVector, nodes...);
assembleImpl(data, std::forward<decltype(quad)>(quad), elementMatrixVector, nodes...);
return true;
}
......@@ -199,7 +196,7 @@ namespace AMDiS
/// Generate a \ref LocalAssembler on a given `LocalContext` (element or intersection)
template <class LocalContext, class Operator, class... Nodes,
std::enable_if_t<not traits::is_reference_wrapper<Operator>::value, int> = 0>
std::enable_if_t<not Traits::IsReferenceWrapper<Operator>::value, int> = 0>
auto makeLocalAssembler(Operator const& op, Nodes const&...)
{
return LocalAssembler<LocalContext, Operator, Nodes...>{op};
......@@ -213,7 +210,7 @@ namespace AMDiS
/// Generate a shared_ptr to \ref LocalAssembler on a given `LocalContext` (element or intersection)
template <class LocalContext, class Operator, class... Nodes,
std::enable_if_t<not traits::is_reference_wrapper<Operator>::value, int> = 0>
std::enable_if_t<not Traits::IsReferenceWrapper<Operator>::value, int> = 0>
auto makeLocalAssemblerPtr(Operator const& op, Nodes const&...)
{
return std::make_shared<LocalAssembler<LocalContext, Operator, Nodes...>>(op);
......
......@@ -70,7 +70,7 @@ namespace AMDiS
template <class Grid>
class MeshCreator
{
static unique_ptr<Grid> create(std::string meshName)
static std::unique_ptr<Grid> create(std::string meshName)
{
error_exit("Creator not yet implemented for this mesh type.");
}
......@@ -82,7 +82,7 @@ namespace AMDiS
{
using Grid = Dune::AlbertaGrid<dim, dimworld>;
static unique_ptr<Grid> create(std::string meshName)
static std::unique_ptr<Grid> create(std::string meshName)
{
std::string macro_filename = "";
Parameters::get(meshName + "->macro file name", macro_filename);
......@@ -90,7 +90,7 @@ namespace AMDiS
// TODO: if filename_extension is ".2d" or ".3d" read it directly from file
// otherwise use a factory method
return make_unique<Grid>(macro_filename);
return std::make_unique<Grid>(macro_filename);
}
};
#endif
......@@ -102,7 +102,7 @@ namespace AMDiS
{
using Grid = Dune::UGGrid<dim>;
static unique_ptr<Grid> create(std::string meshName)
static std::unique_ptr<Grid> create(std::string meshName)
{
std::string filename = "";
......@@ -117,12 +117,12 @@ namespace AMDiS
Dune::GridFactory<Grid> factory;
Dune::AlbertaReader<Grid> reader;
reader.readGrid(filename, factory);
return unique_ptr<Grid>{factory.createGrid()};
return std::unique_ptr<Grid>{factory.createGrid()};
}
#endif
if (ext == "msh") {
Dune::GmshReader<Grid> reader;
return unique_ptr<Grid>{reader.read(filename)};
return std::unique_ptr<Grid>{reader.read(filename)};
}
} else {
error_exit("Construction of UGGrid without filename not yet implemented!");
......@@ -139,7 +139,7 @@ namespace AMDiS
{
using Grid = Dune::YaspGrid<dim, Dune::EquidistantCoordinates<T,dim>>;
static unique_ptr<Grid> create(std::string meshName)
static std::unique_ptr<Grid> create(std::string meshName)
{
Dune::FieldVector<double, dim> L; L = 1.0; // extension of the domain
Parameters::get(meshName + "->dimension", L);
......@@ -151,7 +151,7 @@ namespace AMDiS
// TODO: add more parameters for yasp-grid (see constructor)
return make_unique<Grid>(L, s);
return std::make_unique<Grid>(L, s);
}
};
......@@ -161,7 +161,7 @@ namespace AMDiS
{
using Grid = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<T, dim>>;
static unique_ptr<Grid> create(std::string meshName)
static std::unique_ptr<Grid> create(std::string meshName)
{
Dune::FieldVector<double, dim> lowerleft; lowerleft = 0.0; // Lower left corner of the domain