Commit 01aefb50 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Removed tuple implementation and used wrapper instead

parent e37c96cd
......@@ -54,6 +54,8 @@ namespace AMDiS
SystemVectorType& rhs,
bool asmMatrix, bool asmVector) const;
/// Assemble operators on an element, by passing the element/intersection to
/// `elementAssembler` functor.
template <class Element, class Operators, class ElementAssembler>
void assembleElementOperators(
Element const& element,
......@@ -70,14 +72,14 @@ namespace AMDiS
bool asmMatrix, bool asmVector) const;
/// Return whether the matrix-block needs to be assembled
/// Return the element the LocalViews are bound to
template <class LocalView0, class... LovalViews>
auto const& getElement(LocalView0 const& localView, LovalViews const&...) const
{
return localView.element();
}
/// Return whether the matrix-block needs to be assembled
/// Return the gridView the localViews are bound to
template <class LocalView0, class... LovalViews>
auto const& getGridView(LocalView0 const& localView, LovalViews const&...) const
{
......
#pragma once
#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/IndexSeq.hpp>
#include <amdis/common/Loops.hpp>
namespace AMDiS
{
template <class FeSpace>
struct LocalView
{
using type = typename FeSpace::LocalView;
};
template <class FeSpace>
struct LocalIndexSet
{
using type = typename FeSpace::LocalIndexSet;
};
template <class FeSpaces>
class FiniteElementSpaces
{
template <std::size_t I>
using FeSpace = std::tuple_element_t<I, FeSpaces>;
static constexpr int nComponents = std::tuple_size<FeSpaces>::value;
static_assert( nComponents > 0, "" );
using LocalViews = MapTuple_t<LocalView, FeSpaces>;
using LocalIndexSets = MapTuple_t<LocalIndexSet, FeSpaces>;
/// The grid view the global FE basis lives on
using GridView = typename FeSpace<0>::GridView;
/// Type of the grid element we are bound to
using Element = typename GridView::template Codim<0>::Entity;
public:
explicit FiniteElementSpaces(std::shared_ptr<FeSpaces> const& feSpaces)
: feSpaces_(feSpaces)
, localViews_(mapTuple([](auto const& basis) { return basis.localView(); }, *feSpaces))
, localIndexSets_(mapTuple([](auto const& basis) { return basis.localIndexSet(); }, *feSpaces))
{}
/// Update all global bases. This will update the indexing information of the global basis.
/// NOTE: It must be called if the grid has changed.
void update(GridView const& gv)
{
forEach(range_<0, nComponents>, [&,this](auto const _i)
{
static const int I = decltype(_i)::value;
std::get<I>(*feSpaces_).update(gv);
});
}
/// Bind the \ref localViews and \ref localIndexSets to a grid element
void bind(Element const& element)
{
forEach(range_<0, nComponents>, [&,this](auto const _i)
{
static const int I = decltype(_i)::value;
auto& localView = std::get<I>(localViews_);
localView.bind(element);
auto& localIndexSet = std::get<I>(localIndexSets_);
localIndexSet.bind(localView);
});
// NOTE: maybe create element-geometry here
bound_ = true;
}
/// Unbind from the current element
void unbind()
{
forEach(range_<0, nComponents>, [&,this](auto const _i)
{
static const int I = decltype(_i)::value;
std::get<I>(localIndexSets_).unbind();
std::get<I>(localViews_).unbind();
});
bound_ = false;
}
template <std::size_t I>
auto const& feSpace(const index_t<I> _i = {}) const
{
return std::get<I>(*feSpaces_);
}
template <std::size_t I>
auto const& localView(const index_t<I> _i = {}) const
{
assert( bound_ && "localViews must be bound to an element." );
return std::get<I>(localViews_);
}
template <std::size_t I>
auto const& localIndexSet(const index_t<I> _i = {}) const
{
assert( bound_ && "localIndexSets must be bound to a localView." );
return std::get<I>(localIndexSets_);
}
auto const& element() const
{
assert( bound_ && "localViews must be bound to an element." );
return std::get<0>(localViews_).element();
}
auto const& gridView() const
{
return std::get<0>(*feSpaces_).gridView();
}
private:
/// Tuple of global functionspace bases
std::shared_ptr<FeSpaces> feSpaces_;
/// Tuple of localView objects, obtained from the tuple of global bases
LocalViews localViews_;
/// Tuple of localIndexSet objects, obtained from the tuple of global bases
LocalIndexSets localIndexSets_;
/// True, if localViews and localIndexSets are bound to an element
bool bound_ = false;
};
} // end namespace AMDiS
......@@ -16,7 +16,8 @@ namespace AMDiS
* @{
**/
/// \brief The main implementation of an operator to be used in a \ref LocalAssembler.
/// \brief The main implementation of an CRTP-base class for operators using a grid-function
/// coefficient 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
......@@ -24,6 +25,7 @@ namespace AMDiS
*
* The class is specialized, by deriving from it, in \ref GridFunctionOperator.
*
* \tparam Derived The class derived from GridFunctionOperatorBase
* \tparam LocalContext The Element or Intersection type
* \tparam GridFunction The GridFunction, a LocalFunction is created from, and
* that is evaluated at quadrature points.
......@@ -36,21 +38,26 @@ namespace AMDiS
class GridFunctionOperatorBase
: public LocalOperator<Derived, LocalContext>
{
/// The type of the localFunction associated with the GridFunction
using LocalFunction = decltype(localFunction(std::declval<GridFunction>()));
/// The Codim=0 entity of the grid, the localFunction can be bound to
using Element = typename Impl::ContextTypes<LocalContext>::Entity;
/// The geometry-type of the grid element
using Geometry = typename Element::Geometry;
/// The type of the local coordinates in the \ref Element
using LocalCoordinate = typename GridFunction::EntitySet::LocalCoordinate;
/// A factory for QuadratureRules that incooperate the order of the LocalFunction
using QuadFactory = QuadratureFactory<typename Geometry::ctype, LocalContext::mydimension, LocalFunction>;
public:
/// \brief Constructor. Stores a copy of `gridFct`.
/**
* An expression operator takes an expression, following the interface of
* \ref ExpressionBase, and stores a copy. Additionally, it gets the
* differentiation order, to calculate the quadrature degree in \ref getDegree.
* A GridFunctionOperator takes a gridFunction and the
* differentiation order of the operator, to calculate the
* quadrature degree in \ref getDegree.
**/
GridFunctionOperatorBase(GridFunction const& gridFct, int termOrder)
: gridFct_(gridFct)
......@@ -78,7 +85,8 @@ namespace AMDiS
* Since all operators need geometry information, the `Element::Geometry`
* object `geometry` is created once, during grid traversal, and passed in.
*
* By default, it binds the \ref localFct_ to the `element`.
* By default, it binds the \ref localFct_ to the `element` and the Quadrature
* factory to the localFunction.
**/
void bind_impl(Element const& element, Geometry const& geometry)
{
......@@ -93,6 +101,7 @@ namespace AMDiS
localFct_.unbind();
}
/// Create a quadrature factory from a PreQuadratureFactory, e.g. class derived from \ref QuadratureFactory
template <class PreQuadFactory>
void setQuadFactory(PreQuadFactory const& pre)
{
......@@ -107,7 +116,7 @@ namespace AMDiS
return localFct_(local);
}
/// Create a quadrature rule using the \ref QuadratureCreator by calculating the
/// Create a quadrature rule using the \ref QuadratureFactory by calculating the
/// quadrature order necessary to integrate the (bi)linear-form accurately.
template <class... Nodes>
auto const& getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
......@@ -148,23 +157,27 @@ namespace AMDiS
: transposedOp_(std::forward<Args>(args)...)
{}
/// Redirects the bind call top the transposed operator
template <class Element, class Geometry>
void bind_impl(Element const& element, Geometry const& geometry)
{
transposedOp_.bind_impl(element, geometry);
}
/// Redirects the unbind call top the transposed operator
void unbind_impl()
{
transposedOp_.unbind_impl();
}
/// Redirects the setQuadFactory call top the transposed operator
template <class PreQuadFactory>
void setQuadFactory(PreQuadFactory const& pre)
{
transposedOp_.setQuadFactory(pre);
}
/// Apply the assembling to the transposed elementMatrix with interchanged row-/colNode
template <class Context, class RowNode, class ColNode, class ElementMatrix>
void getElementMatrix(Context const& context,
RowNode const& rowNode,
......
#pragma once
#include <list>
#include <memory>
#include <type_traits>
#include <amdis/ContextGeometry.hpp>
#include <amdis/common/ConceptsBase.hpp>
// #include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/TypeDefs.hpp>
#include <amdis/utility/TreeData.hpp>
namespace AMDiS
{
......@@ -16,24 +13,28 @@ namespace AMDiS
class LocalAssemblerBase
{
public:
/// The codim=0 grid entity
using Element = typename Impl::ContextTypes<LocalContext>::Entity;
/// The geometry of the \ref Element
using Geometry = typename Element::Geometry;
static constexpr int numNodes = sizeof...(Nodes);
static_assert( numNodes == 1 || numNodes == 2,
"VectorAssembler gets 1 Node, MatrixAssembler gets 2 Nodes!");
/// Either an ElementVector or an ElementMatrix (depending on the number of nodes)
using ElementMatrixVector = std::conditional_t<
(sizeof...(Nodes)==1), Impl::ElementVector, std::conditional_t<
(sizeof...(Nodes)==2), Impl::ElementMatrix, void>>;
public:
/// Virtual destructor
virtual ~LocalAssemblerBase() {}
/// Bind the local-assembler to the grid-element with its corresponding geometry
virtual void bind(Element const& element, Geometry const& geometry) = 0;
/// Unbind from the element
virtual void unbind() = 0;
/// Assemble an element matrix or element vector on the test- (and trial-) function node(s)
......@@ -42,74 +43,4 @@ namespace AMDiS
ElementMatrixVector& elementMatrixVector) = 0;
};
template <class GridView>
class OperatorLists
{
using Element = typename GridView::template Codim<0>::Entity;
template <class OperatorType>
struct Scaled
{
std::shared_ptr<OperatorType> op;
double* factor = nullptr;
double* estFactor = nullptr;
BoundaryType b = {0};
};
template <class... Nodes>
struct Data
{
using ElementOperator = LocalAssemblerBase<Element, Nodes...>;
using IntersectionOperator = LocalAssemblerBase<typename GridView::Intersection, Nodes...>;
std::list<Scaled<ElementOperator>> element;
std::list<Scaled<IntersectionOperator>> boundary;
std::list<Scaled<IntersectionOperator>> intersection;
bool assembled = false; // if false, do reassemble
bool changing = false; // if true, or assembled false, do reassemble
bool empty() const
{
return element.empty() && boundary.empty() && intersection.empty();
}
bool doAssemble(bool flag) const
{
return flag && (!assembled || changing);
}
template <class Geo>
void bind(Element const& elem, Geo const& geo)
{
for (auto& scaled : element) scaled.op->bind(elem,geo);
for (auto& scaled : boundary) scaled.op->bind(elem,geo);
for (auto& scaled : intersection) scaled.op->bind(elem,geo);
}
void unbind()
{
for (auto& scaled : element) scaled.op->unbind();
for (auto& scaled : boundary) scaled.op->unbind();
for (auto& scaled : intersection) scaled.op->unbind();
}
};
public:
template <class RowNode, class ColNode>
using MatData = Data<RowNode, ColNode>;
template <class Node>
using VecData = Data<Node>;
};
template <class GlobalBasis>
using MatrixOperators = MatrixData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template MatData>;
template <class GlobalBasis>
using VectorOperators = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template VecData>;
} // end namespace AMDiS
#pragma once
#include <list>
#include <memory>
#include <amdis/LocalAssemblerBase.hpp>
#include <amdis/utility/TreeData.hpp>
namespace AMDiS
{
template <class GridView>
class OperatorLists
{
using Element = typename GridView::template Codim<0>::Entity;
using Intersection = typename GridView::Intersection;
template <class OperatorType>
struct DataElement
{
std::shared_ptr<OperatorType> op;
double* factor = nullptr;
double* estFactor = nullptr;
BoundaryType b = {0};
};
/// Lists of \ref DataElement on the Element, BoundaryIntersction, and InteriorIntersections
template <class... Nodes>
struct Data
{
/// The type of local operators associated with grid elements
using ElementOperator = LocalAssemblerBase<Element, Nodes...>;
/// The type of local operators associated with grid intersections
using IntersectionOperator = LocalAssemblerBase<Intersection, Nodes...>;
/// Return whether any operators are stored on the node
bool empty() const
{
return element.empty() && boundary.empty() && intersection.empty();
}
/// Test whether to assemble on the node
bool doAssemble(bool flag) const
{
return flag && (!assembled || changing);
}
/// Bind all operators to the grid element and geometry
template <class Geo>
void bind(Element const& elem, Geo const& geo)
{
for (auto& scaled : element) scaled.op->bind(elem,geo);
for (auto& scaled : boundary) scaled.op->bind(elem,geo);
for (auto& scaled : intersection) scaled.op->bind(elem,geo);
}
/// Unbind all operators from the element
void unbind()
{
for (auto& scaled : element) scaled.op->unbind();
for (auto& scaled : boundary) scaled.op->unbind();
for (auto& scaled : intersection) scaled.op->unbind();
}
/// List of operators to be assembled on grid elements
std::list<DataElement<ElementOperator>> element;
/// List of operators to be assembled on boundary intersections
std::list<DataElement<IntersectionOperator>> boundary;
/// List of operators to be assembled on interior intersections
std::list<DataElement<IntersectionOperator>> intersection;
/// if false, do reassemble of all operators
bool assembled = false;
/// if true, or \ref assembled false, do reassemble of all operators
bool changing = false;
};
public:
/// List of operators associated with matrix blocks (RowNode, ColNode)
template <class RowNode, class ColNode>
using MatData = Data<RowNode, ColNode>;
/// List of operators associated with vector blocks [Node]
template <class Node>
using VecData = Data<Node>;
};
template <class GlobalBasis>
using MatrixOperators = MatrixData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template MatData>;
template <class GlobalBasis>
using VectorOperators = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template VecData>;
} // end namespace AMDiS
......@@ -25,8 +25,11 @@ namespace AMDiS
class LocalOperator
{
public:
/// The element or intersection the operator is assembled on
using LocalContext = LocalContextType;
/// The codim=0 grid entity
using Element = typename Impl::ContextTypes<LocalContext>::Entity;
/// The geometry of the \ref Element
using Geometry = typename Element::Geometry;
/// Initialize the local operator on the current gridView
......@@ -112,19 +115,19 @@ namespace AMDiS
// default implementation called by \ref calculateElementMatrix
template <class Context, class RowNode, class ColNode, class ElementMatrix>
void getElementMatrix(Context const& context,
RowNode const& rowNode,
ColNode const& colNode,
ElementMatrix& elementMatrix)
void getElementMatrix(Context const& /*context*/,
RowNode const& /*rowNode*/,
ColNode const& /*colNode*/,
ElementMatrix& /*elementMatrix*/)
{
error_exit("Needs to be implemented by derived class!");
}
// default implementation called by \ref calculateElementVector
template <class Context, class Node, class ElementVector>
void getElementVector(Context const& context,
Node const& node,
ElementVector& elementVector)
void getElementVector(Context const& /*context*/,
Node const& /*node*/,
ElementVector& /*elementVector*/)
{
error_exit("Needs to be implemented by derived class!");
}
......@@ -186,14 +189,14 @@ namespace AMDiS
};
/// Generate an \ref GridFunctionOperator from a PreOperator.
/// Generate a \ref LocalOperator from a PreOperator.
template <class Derived, class LocalContext, class GridView>
auto makeLocalOperator(LocalOperator<Derived, LocalContext> const& localOp, GridView const& /*gridView*/)
{
return localOp.derived();
}
/// Generate a shared_ptr to \ref GridFunctionOperator from a PreOperator.
/// Generate a shared_ptr to a \ref LocalOperator from a PreOperator.
template <class Derived, class LocalContext, class GridView>
auto makeLocalOperatorPtr(LocalOperator<Derived, LocalContext> const& localOp, GridView const& /*gridView*/)
{
......
......@@ -27,7 +27,7 @@ namespace AMDiS
public:
AnalyticLocalFunction(Function const& fct)
: fct_(fct)
: fct_{fct}
{}
void bind(LocalContext const& element)
......@@ -42,7 +42,7 @@ namespace AMDiS
Range operator()(Domain const& local) const
{
assert(geometry_.has_value());
assert( bool(geometry_) );
return fct_(geometry_.value().global(local));
}
......@@ -108,8 +108,8 @@ namespace AMDiS
public:
/// \brief Constructor. Stores the function `fct` and creates an `EntitySet`.
AnalyticGridFunction(Function const& fct, GridView const& gridView)
: fct_(fct)
, entitySet_(gridView)
: fct_{fct}
, entitySet_{gridView}
{}
/// \brief Return the evaluated functor at global coordinates
......@@ -124,7 +124,6 @@ namespace AMDiS
return {fct_};
}
EntitySet const& entitySet() const
{
return entitySet_;
......
......@@ -297,6 +297,13 @@ namespace AMDiS
return *this;
}
template <class Expr>
DOFVectorMutableView& operator<<(Expr&& expr)
{
return interpolate(expr);
}
/// Return the mutable DOFVector
DOFVector<Traits>& coefficients() { return *mutableDofVector_; }
......
......@@ -45,7 +45,7 @@ namespace AMDiS
public:
/// Constructor. Stores a copy of gridFct.
explicit DerivativeGridFunction(GridFunction const& gridFct)
: gridFct_(gridFct)
: gridFct_{gridFct}
{
static_assert(isValidRange<DerivativeTraits>(), "Derivative of GridFunction not defined");
}
......@@ -58,9 +58,9 @@ namespace AMDiS
}
/// Return the derivative-localFunction of the GridFunction.
friend LocalFunction localFunction(DerivativeGridFunction const& gf)
LocalFunction localFunction() const
{
return derivative(localFunction(gf.gridFct_));
return derivative(localFunction(gridFct_));
}
/// Return the \ref EntitySet of the \ref GridFunction.
......@@ -74,7 +74,13 @@ namespace AMDiS
};
#ifndef DOXYGEN
template <class GF>
auto localFunction(DerivativeGridFunction<GF> const& gf)
{
return gf.localFunction();
}
template <class GridFct,