Commit 7e44b707 authored by Praetorius, Simon's avatar Praetorius, Simon

new operator style implemented, using tags

parent 0057b587
......@@ -14,9 +14,6 @@
namespace AMDiS
{
const Flag EQUAL_BASES = 0x01L;
const Flag EQUAL_LOCALBASES = 0x02L;
template <class Traits>
class Assembler
{
......
......@@ -110,7 +110,7 @@ void Assembler<Traits>::assembleElementOperators(
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()..., scaled.factor);
bool add_op = scaled.op->assemble(context, elementContainer, localViews.tree()...);
scaled.op->unbind();
add = add || add_op;
}
......@@ -148,18 +148,16 @@ void Assembler<Traits>::initMatrixVector(
auto localView = globalBasis_.localView();
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{
#ifdef HAVE_EXTENDED_DUNE_FUNCTIONS
if (rowNode.isLeaf)
msg(globalBasis_.dimension(rowTreePath), " DOFs for Basis[", to_string(rowTreePath), "]");
#endif
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
//if (rhsOperators_[rowNode].assemble(asmVector))
// rhsOperators_[rowNode].init(rowBasis);
forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
//if (matrixOperators_[rowNode][colNode].assemble(asmMatrix))
// matrixOperators_[rowNode][colNode].init(rowBasis, colBasis);
for (auto bc : constraints_[rowNode][colNode].scalar)
bc->init(matrix, solution, rhs, rowBasis, colBasis);
......
#pragma once
namespace AMDiS
{
template <class LocalContext, class Geometry, class LocalGeometry>
struct ContextGeometry
{
enum {
dim = Geometry::mydimension, //< the dimension of the grid element
dow = Geometry::coorddimension //< the dimension of the world
};
LocalContext const& localContext;
Geometry const& geometry;
LocalGeometry const& localGeometry;
/// Coordinate `p` given in `localGeometry`, transformed to coordinate in `geometry`.
template <class Coordinate>
decltype(auto) position(Coordinate const& p) const
{
return position(p, std::is_same<Geometry, LocalGeometry>{});
}
/// The integration element from the `localGeometry`, the quadrature points are
/// defined in.
template <class Coordinate>
auto integrationElement(Coordinate const& p) const
{
return localGeometry.integrationElement(p);
}
/// Transformation of coordinate `p` given in `localGeometry` to world space coordinates.
template <class Coordinate>
decltype(auto) global(Coordinate const& p) const
{
return geometry.global(p);
}
private: // implementation detail
// position for elements
template <class Coordinate>
Coordinate const& position(Coordinate const& p, std::true_type) const
{
return p;
}
// position for intersection
template <class Coordinate>
auto position(Coordinate const& p, std::false_type) const
{
return localGeometry.global(p);
}
};
} // end namespace AMDiS
......@@ -91,7 +91,7 @@ namespace AMDiS
{
using WorldVector = typename GlobalBasis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
template <class Node>
template <class RowNode, class ColNode>
struct type
{
std::list<std::shared_ptr<DirichletBC<WorldVector, double>>> scalar;
......
This diff is collapsed.
This diff is collapsed.
#pragma once
#include <list>
#include <memory>
#include <type_traits>
#include <dune/amdis/common/TypeDefs.hpp>
#include <dune/amdis/utility/GetEntity.hpp>
#include <dune/amdis/utility/TreeData.hpp>
namespace AMDiS
{
/// Abstract base-class of a \ref LocalAssembler
template <class LocalContext, class... Nodes>
class LocalAssemblerBase
{
public:
using Element = typename Impl::Get<LocalContext>::Entity;
using Geometry = typename Element::Geometry;
using LocalGeometry = typename Impl::Get<LocalContext>::Geometry;
static constexpr int numNodes = sizeof...(Nodes);
static_assert( numNodes == 1 || numNodes == 2,
"VectorAssembler gets 1 Node, MatrixAssembler gets 2 Nodes!");
using ElementMatrixVector = std::conditional_t<
(sizeof...(Nodes)==1), Impl::ElementVector, std::conditional_t<
(sizeof...(Nodes)==2), Impl::ElementMatrix, void>>;
public:
virtual void bind(Element const& element, Geometry const& geometry) = 0;
virtual void unbind() = 0;
/// Assemble an element matrix or element vector on the test- (and trial-) function node(s)
virtual bool assemble(LocalContext const& localContext,
ElementMatrixVector& elementMatrixVector,
Nodes const&... nodes) = 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 assemble(bool flag) const
{
return flag && (!assembled || changing);
}
};
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
This diff is collapsed.
#pragma once
#include <algorithm>
#include <list>
namespace AMDiS {
template <class GridView, class LocalContext>
template <class RowNode, class ColNode>
void Operator<GridView, LocalContext>::init(
RowNode const& rowNode,
ColNode const& colNode)
{
psiDegree_ = getPolynomialDegree(rowNode);
phiDegree_ = getPolynomialDegree(colNode);
}
template <class GridView, class LocalContext>
template <class RowNode>
void Operator<GridView, LocalContext>::init(RowNode const& node)
{
psiDegree_ = getPolynomialDegree(node);
}
template <class GridView, class LocalContext>
int Operator<GridView, LocalContext>::
getQuadratureDegree(Dune::GeometryType t, LocalGeometry const& geometry, int order, FirstOrderType type)
{
std::list<OperatorTermType*>* terms = NULL;
switch(order)
{
case 0:
terms = &zeroOrder;
break;
case 1:
if (type == GRD_PHI)
terms = &firstOrderGrdPhi;
else
terms = &firstOrderGrdPsi;
break;
case 2:
terms = &secondOrder;
break;
}
int maxTermDegree = 0;
for (OperatorTermType* term : *terms)
maxTermDegree = std::max(maxTermDegree, term->getDegree());
int degree = psiDegree_ + phiDegree_ + maxTermDegree;
if (t.isSimplex())
degree -= order;
if (!geometry.affine())
degree += 1; // oder += (order+1)
return degree;
}
template <class GridView, class LocalContext>
template <class RowNode, class ColNode>
bool Operator<GridView, LocalContext>::assemble(
LocalContext const& context,
ElementMatrix& elementMatrix,
RowNode const& rowNode,
ColNode const& colNode,
double* factor)
{
double fac = factor ? *factor : 1.0;
if (fac == 0.0)
return false;
init(rowNode, colNode);
auto&& geometry = getGeometry();
auto&& quad_geometry = getLocalGeometry(context);
if (!zeroOrder.empty()) {
zeroOrderAssembler.bind(fac, *this, context, geometry, quad_geometry);
zeroOrderAssembler.calculateElementMatrix(*this, elementMatrix, rowNode, colNode);
zeroOrderAssembler.unbind();
}
if (!firstOrderGrdPhi.empty()) {
firstOrderGrdPhiAssembler.bind(fac, *this, context, geometry, quad_geometry);
firstOrderGrdPhiAssembler.calculateElementMatrix(*this, elementMatrix, rowNode, colNode);
firstOrderGrdPhiAssembler.unbind();
}
if (!firstOrderGrdPsi.empty()) {
firstOrderGrdPsiAssembler.bind(fac, *this, context, geometry, quad_geometry);
firstOrderGrdPsiAssembler.calculateElementMatrix(*this, elementMatrix, rowNode, colNode);
firstOrderGrdPsiAssembler.unbind();
}
if (!secondOrder.empty()) {
secondOrderAssembler.bind(fac, *this, context, geometry, quad_geometry);
secondOrderAssembler.calculateElementMatrix(*this, elementMatrix, rowNode, colNode);
secondOrderAssembler.unbind();
}
return true;
}
template <class GridView, class LocalContext>
template <class RowNode>
bool Operator<GridView, LocalContext>::assemble(
LocalContext const& context,
ElementVector& elementVector,
RowNode const& node,
double* factor)
{
double fac = factor ? *factor : 1.0;
if (fac == 0.0)
return false;
test_exit(firstOrderGrdPhi.empty() && secondOrder.empty(),
"Derivatives on ansatz-functions not allowed on the vector-side!");
init(node);
auto&& geometry = getGeometry();
auto&& quad_geometry = getLocalGeometry(context);
if (!zeroOrder.empty()) {
zeroOrderAssembler.bind(fac, *this, context, geometry, quad_geometry);
zeroOrderAssembler.calculateElementVector(*this, elementVector, node);
zeroOrderAssembler.unbind();
}
if (!firstOrderGrdPsi.empty()) {
firstOrderGrdPsiAssembler.bind(fac, *this, context, geometry, quad_geometry);
firstOrderGrdPsiAssembler.calculateElementVector(*this, elementVector, node);
firstOrderGrdPsiAssembler.unbind();
}
return true;
}
template <class GridView, class LocalContext>
template <class Term>
Operator<GridView, LocalContext>& Operator<GridView, LocalContext>::
addZOTImpl(Term const& term)
{
zeroOrder.push_back(new GenericOperatorTerm<LocalContext, Term>(term));
return *this;
}
template <class GridView, class LocalContext>
template <class Term>
Operator<GridView, LocalContext>& Operator<GridView, LocalContext>::
addFOTImpl(Term const& term, FirstOrderType type)
{
using OpTerm = GenericOperatorTerm<LocalContext, Term>;
if (type == GRD_PHI)
firstOrderGrdPhi.push_back(new OpTerm(term));
else
firstOrderGrdPsi.push_back(new OpTerm(term));
return *this;
}
template <class GridView, class LocalContext>
template <class Term>
Operator<GridView, LocalContext>& Operator<GridView, LocalContext>::
addFOTImpl(Term const& term, std::size_t i, FirstOrderType type)
{
using OpTerm = GenericOperatorTerm<LocalContext, Term, VectorComponent>;
if (type == GRD_PHI)
firstOrderGrdPhi.push_back(new OpTerm(term, {i}));
else
firstOrderGrdPsi.push_back(new OpTerm(term, {i}));
return *this;
}
template <class GridView, class LocalContext>
template <class Term>
Operator<GridView, LocalContext>& Operator<GridView, LocalContext>::
addSOTImpl(Term const& term)
{
secondOrder.push_back(new GenericOperatorTerm<LocalContext, Term>(term));
return *this;
}
template <class GridView, class LocalContext>
template <class Term>
Operator<GridView, LocalContext>& Operator<GridView, LocalContext>::
addSOTImpl(Term const& term, std::size_t i, std::size_t j)
{
using OpTerm = GenericOperatorTerm<LocalContext, Term, MatrixComponent>;
secondOrder.push_back(new OpTerm(term, {i,j}));
return *this;
}
} // end namespace AMDiS
#pragma once
#include <cassert>
#include <dune/amdis/terms/TermGenerator.hpp>
#include <dune/amdis/utility/GetDegree.hpp>
namespace AMDiS
{
template <class Expr>
class ExpressionOperatorBase
{
public:
/// \brief Constructor. Stores a copy of `expr`.
/**
* 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.
**/
ExpressionOperatorBase(Expr const& expr, int order)
: expr_(expr)
, order_(order)
{}
/// \brief Binds operator to `element` and `geometry`.
/**
* Binding an operator to the currently visited element in grid traversal.
* 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 expr_ to the `element`.
**/
template <class Element, class Geometry>
void bind(Element const& element, Geometry const& geometry)
{
if (!bound_) {
expr_.bind(element);
isSimplex_ = geometry.type().isSimplex();
isAffine_ = geometry.affine();
bound_ = true;
}
}
/// Unbinds operator from element.
void unbind()
{
bound_ = false;
initSize_ = 0;
expr_.unbind();
}
/// \brief Evaluate expression at quadrature points.
/**
* Initialization at all points `quad[iq].position()` where
* `iq` in `[0, quad.size())`.
*
* \param context Container for localContext, localGeometry and geometry, see \ref LocalAssembler
* \param quad Liste of evaluations points. Follows interface of \ref Dune::QuadratureRule
**/
template <class Context, class Quad>
void init(Context const& context, Quad const& quad)
{
if (initSize_ != quad.size()) {
expr_.init(context, quad);
initSize_ = quad.size();
}
}
/// Return expressions iq'th value. Must be initialized in \ref init before.
auto eval(std::size_t iq) const
{
assert( initSize_ > 0u );
return expr_[iq];
}
/// \brief Return the quadrature degree for a vector operator.
/**
* The quadrature degree that is necessary, to integrate the expression
* multiplied with (possibly derivatives of) shape functions. Thus it needs
* the order of derivatives, this operator implements.
**/
template <class Node>
int getDegree(Node const& node) const
{
assert( bound_ );
int psiDegree = getPolynomialDegree(node);
int termDegree = expr_.getDegree();
int degree = psiDegree + termDegree;
if (isSimplex_)
degree -= order_;
if (isAffine_)
degree += 1; // oder += (order+1)
return degree;
}
/// \brief Return the quadrature degree for a matrix operator.
/**
* The quadrature degree that is necessary, to integrate the expression
* multiplied with (possibly derivatives of) shape functions. Thus it needs
* the order of derivatives, this operator implements.
**/
template <class RowNode, class ColNode>
int getDegree(RowNode const& rowNode, ColNode const& colNode) const
{
assert( bound_ );
int psiDegree = getPolynomialDegree(rowNode);
int phiDegree = getPolynomialDegree(colNode);
int termDegree = expr_.getDegree();
int degree = psiDegree + phiDegree + termDegree;
if (isSimplex_)
degree -= order_;
if (isAffine_)
degree += 1; // oder += (order+1)
return degree;
}
protected:
Expr& expr()
{
return expr_;
}
Expr const& expr() const
{
return expr_;
}
private:
Expr expr_; //< the stores expression, see \ref ExpressionBase
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
bool bound_ = false; //< an element is bound to the operator
std::size_t initSize_ = 0; //< number of quadrature points, the expression was evaluated at.
};
/// A default implementation of an ExpressionOperator if no specialization is available.
/**
* An operator must implement either \ref calculateElementVector, or
* \ref calculateElementMatrix, if it is a vector or matrix operator,
* respectively.
**/
template <class Tag, class Expr>
class ExpressionOperator
: public ExpressionOperatorBase<Expr>
{
public:
ExpressionOperator(Tag, Expr const& expr)
: ExpressionOperatorBase<Expr>(expr, 0)
{}
/// Assemble a local element vector on the element that is bound.
template <class Context, class QuadratureRule,
class ElementVector, class Node>
void calculateElementVector(Context const&, //< container for context related data
QuadratureRule const&, //< a quadrature formula
ElementVector&, //< the output element vector
Node const& //< the node of the test-basis
)
{
error_exit("Needs to be implemented!");
}
/// Assemble a local element matrix on the element that is bound.
template <class Context, class QuadratureRule,
class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
void calculateElementMatrix(Context const&, //< container for context related data
QuadratureRule const&, //< a quadrature formula
ElementMatrix&, //< the output element matrix
RowNode const&, //< the node of the test-basis
ColNode const&, //< the node of the trial-basis
std::integral_constant<bool, sameFE>, //< finiteelement is the same for test and trial
std::integral_constant<bool, sameNode> //< nodes are the same in the basis-tree
)
{
error_exit("Needs to be implemented!");
}
};
/// Generate an \ref ExpressionOperator with a given tag `t` and pre-expression `epxr`.
template <class Tag, class Expr_>
auto makeOperator(Tag t, Expr_ const& expr)
{
using Expr = ToTerm_t<Expr_>;
return ExpressionOperator<Tag, Expr>{t, toTerm(expr)};
}
/// Generate a shared_ptr to \ref ExpressionOperator with a given tag `t` and pre-expression `epxr`.
template <class Tag, class Expr_>
auto makeOperatorPtr(Tag t, Expr_ const& expr)
{
using Expr = ToTerm_t<Expr_>;
return std::make_shared<ExpressionOperator<Tag, Expr>>(t, toTerm(expr));
}
} // end namespace AMDiS
#pragma once
#include <cassert>
#include <dune/amdis/common/FieldMatVec.hpp>
#include <dune/amdis/common/ValueCategory.hpp>
namespace AMDiS {
struct Evaluate
{
using Scalar = Dune::FieldVector<double, 1>;
template <class... Args>
static auto zot(Args...) { assert( false ); return 0; }
template <class ValueType>
static auto zot(tag::scalar, tag::none, ValueType const& v,
Scalar const& test, Scalar const& trial)
{
return v * test * trial;
}
template <class ValueType>
static auto zot(tag::scalar, tag::none, ValueType const& v,
Scalar const& test)
{
return v * test;
}
template <class... Args>
static auto fot(Args...) { assert( false ); return 0; }
template <class ValueType, class Gradient>
static auto fot(tag::scalar, VectorComponent comp, ValueType const& v,
Gradient const& grad_test, Scalar const& trial)
{
return v * grad_test[comp.index] * trial;
}
template <class ValueType, class Gradient>
static auto fot(tag::vector, tag::none, ValueType const& v,
Gradient const& grad_test, Scalar const& trial)
{
return (v * grad_test) * trial;
}
template <class ValueType, class Gradient>
static auto fot(tag::scalar, tag::none, ValueType const& v,
Gradient const& grad_test, Scalar const& trial)
{
return v * sum(grad_test) * trial;
}
template <class ValueType, class Gradient>
static auto fot(tag::scalar, VectorComponent comp, ValueType const& v,
Gradient const& grad_test)
{
return v * grad_test[comp.index] ;
}
template <class ValueType, class Gradient>
static auto fot(tag::vector, tag::none, ValueType const& v,
Gradient const& grad_test)
{
return (v * grad_test);
}
template <class ValueType, class Gradient>
static auto fot(tag::scalar, tag::none, ValueType const& v,
Gradient const& grad_test)
{
return v * sum(grad_test);
}
template <clas