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

corrected error in quadrature degree, restructuring of zero/first/secondOrderAssembler

parent 781a6094
...@@ -24,6 +24,7 @@ namespace AMDiS ...@@ -24,6 +24,7 @@ namespace AMDiS
using GridView = typename GlobalBasis::GridView; using GridView = typename GlobalBasis::GridView;
public: public:
/// Constructor, stores a shared-pointer to the feSpaces /// Constructor, stores a shared-pointer to the feSpaces
Assembler(GlobalBasis& globalBasis, Assembler(GlobalBasis& globalBasis,
MatrixOperators<GlobalBasis>& matrixOperators, MatrixOperators<GlobalBasis>& matrixOperators,
...@@ -48,7 +49,9 @@ namespace AMDiS ...@@ -48,7 +49,9 @@ namespace AMDiS
SystemVectorType& rhs, SystemVectorType& rhs,
bool asmMatrix, bool asmVector); bool asmMatrix, bool asmVector);
private: private:
/// Sets the system to zero and initializes all operators and boundary conditions /// Sets the system to zero and initializes all operators and boundary conditions
template <class SystemMatrixType, class SystemVectorType> template <class SystemMatrixType, class SystemVectorType>
void initMatrixVector( void initMatrixVector(
...@@ -58,11 +61,12 @@ namespace AMDiS ...@@ -58,11 +61,12 @@ namespace AMDiS
bool asmMatrix, bool asmVector) const; bool asmMatrix, bool asmVector) const;
template <class ElementContainer, class Container, class Operators, class... Bases> template <class ElementContainer, class Container, class Operators, class Geometry, class... Bases>
void assembleElementOperators( void assembleElementOperators(
ElementContainer& elementContainer, ElementContainer& elementContainer,
Container& container, Container& container,
Operators& operators, Operators& operators,
Geometry const& geometry,
Bases const&... subBases) const; Bases const&... subBases) const;
...@@ -90,22 +94,9 @@ namespace AMDiS ...@@ -90,22 +94,9 @@ namespace AMDiS
return globalBasis_.gridView(); return globalBasis_.gridView();
} }
public:
template <class RowNode, class ColNode>
Flag optimizationFlags(RowNode const& rowNode, ColNode const& colNode) const
{
Flag flag;
if (rowNode.treeIndex() == colNode.treeIndex())
flag |= EQUAL_BASES;
// NOTE: find a way to compare localBases directly
// if (rowNode.finiteElement().localBasis().order() == colNode.finiteElement().localBasis().order())
// flag |= EQUAL_LOCALBASES;
return flag;
}
private: private:
GlobalBasis& globalBasis_; GlobalBasis& globalBasis_;
MatrixOperators<GlobalBasis>& matrixOperators_; MatrixOperators<GlobalBasis>& matrixOperators_;
VectorOperators<GlobalBasis>& rhsOperators_; VectorOperators<GlobalBasis>& rhsOperators_;
......
...@@ -23,8 +23,8 @@ void Assembler<GlobalBasis>::assemble( ...@@ -23,8 +23,8 @@ void Assembler<GlobalBasis>::assemble(
// 2. create a local matrix and vector // 2. create a local matrix and vector
std::size_t localSize = localView.maxSize(); std::size_t localSize = localView.maxSize();
mtl::dense2D<double> elementMatrix(localSize, localSize); Impl::ElementMatrix elementMatrix(localSize, localSize);
mtl::dense_vector<double> elementVector(localSize); Impl::ElementVector elementVector(localSize);
// 3. traverse grid and assemble operators on the elements // 3. traverse grid and assemble operators on the elements
for (auto const& element : elements(globalBasis_.gridView())) for (auto const& element : elements(globalBasis_.gridView()))
...@@ -33,6 +33,7 @@ void Assembler<GlobalBasis>::assemble( ...@@ -33,6 +33,7 @@ void Assembler<GlobalBasis>::assemble(
set_to_zero(elementVector); set_to_zero(elementVector);
localView.bind(element); localView.bind(element);
auto geometry = element.geometry();
// traverse type-tree of global-basis // traverse type-tree of global-basis
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath) forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
...@@ -43,7 +44,7 @@ void Assembler<GlobalBasis>::assemble( ...@@ -43,7 +44,7 @@ void Assembler<GlobalBasis>::assemble(
auto& rhsOp = rhsOperators_[rowNode]; auto& rhsOp = rhsOperators_[rowNode];
if (rhsOp.assemble(asmVector) && !rhsOp.empty()) if (rhsOp.assemble(asmVector) && !rhsOp.empty())
assembleElementOperators(elementVector, rhs, rhsOp, rowLocalView); assembleElementOperators(elementVector, rhs, rhsOp, geometry, rowLocalView);
forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath) forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{ {
...@@ -53,8 +54,7 @@ void Assembler<GlobalBasis>::assemble( ...@@ -53,8 +54,7 @@ void Assembler<GlobalBasis>::assemble(
auto colLocalView = colBasis.localView(); auto colLocalView = colBasis.localView();
colLocalView.bind(element); // NOTE: Is this necessary colLocalView.bind(element); // NOTE: Is this necessary
assembleElementOperators(elementMatrix, matrix, matOp, rowLocalView, colLocalView, assembleElementOperators(elementMatrix, matrix, matOp, geometry, rowLocalView, colLocalView);
optimizationFlags(rowNode, colNode));
} }
}); });
}); });
...@@ -92,11 +92,12 @@ void Assembler<GlobalBasis>::assemble( ...@@ -92,11 +92,12 @@ void Assembler<GlobalBasis>::assemble(
template <class GlobalBasis> template <class GlobalBasis>
template <class ElementContainer, class Container, class Operators, class... LocalViews> template <class ElementContainer, class Container, class Operators, class Geometry, class... LocalViews>
void Assembler<GlobalBasis>::assembleElementOperators( void Assembler<GlobalBasis>::assembleElementOperators(
ElementContainer& elementContainer, ElementContainer& elementContainer,
Container& container, Container& container,
Operators& operators, Operators& operators,
Geometry const& geometry,
LocalViews const&... localViews) const LocalViews const&... localViews) const
{ {
auto const& element = getElement(localViews...); auto const& element = getElement(localViews...);
...@@ -106,7 +107,9 @@ void Assembler<GlobalBasis>::assembleElementOperators( ...@@ -106,7 +107,9 @@ void Assembler<GlobalBasis>::assembleElementOperators(
auto assemble_operators = [&](auto const& context, auto& operator_list) { auto assemble_operators = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list) { for (auto scaled : operator_list) {
bool add_op = scaled.op->assemble(context, elementContainer, localViews..., scaled.factor); scaled.op->bind(element, geometry);
bool add_op = scaled.op->assemble(context, elementContainer, localViews.tree()..., scaled.factor);
scaled.op->unbind();
add = add || add_op; add = add || add_op;
} }
}; };
...@@ -147,14 +150,14 @@ void Assembler<GlobalBasis>::initMatrixVector( ...@@ -147,14 +150,14 @@ void Assembler<GlobalBasis>::initMatrixVector(
msg(0, " DOFs for Basis[", to_string(rowTreePath), "]"); // TODO: add right values msg(0, " DOFs for Basis[", to_string(rowTreePath), "]"); // TODO: add right values
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath); auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
if (rhsOperators_[rowNode].assemble(asmVector)) //if (rhsOperators_[rowNode].assemble(asmVector))
rhsOperators_[rowNode].init(rowBasis); // rhsOperators_[rowNode].init(rowBasis);
forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath) forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{ {
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath); auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
if (matrixOperators_[rowNode][colNode].assemble(asmMatrix)) //if (matrixOperators_[rowNode][colNode].assemble(asmMatrix))
matrixOperators_[rowNode][colNode].init(rowBasis, colBasis); // matrixOperators_[rowNode][colNode].init(rowBasis, colBasis);
for (auto bc : constraints_[rowNode][colNode].scalar) for (auto bc : constraints_[rowNode][colNode].scalar)
bc->init(matrix, solution, rhs, rowBasis, colBasis); bc->init(matrix, solution, rhs, rowBasis, colBasis);
......
This diff is collapsed.
...@@ -18,19 +18,36 @@ namespace AMDiS ...@@ -18,19 +18,36 @@ namespace AMDiS
using LocalQuadratureRules = Dune::QuadratureRules<ctype, LocalContext::mydimension>; using LocalQuadratureRules = Dune::QuadratureRules<ctype, LocalContext::mydimension>;
using QuadratureRule = QuadratureRuleFactory_t<LocalContext, ctype, dim>; using QuadratureRule = QuadratureRuleFactory_t<LocalContext, ctype, dim>;
using Geometry = typename Impl::Get<LocalContext>::Geometry;
using Geometry = typename GridView::template Codim<0>::Entity::Geometry;
using LocalGeometry = typename Impl::Get<LocalContext>::Geometry;
public: public:
LocalAssembler(int degree, FirstOrderType type = GRD_PHI)
: degree_(degree)
, type_(type)
{}
/// Constructor, initializes the geometry of the element and a /// Constructor, initializes the geometry of the element and a
/// quadrature-rule wrapper /// quadrature-rule wrapper
template <class Operator> template <class Operator>
LocalAssembler(Operator& op, LocalContext const& element, int degree, FirstOrderType type = GRD_PHI) void bind(Operator& op, LocalContext const& localContext, Geometry const& geometry, LocalGeometry const& localGeometry)
: geometry(get_geometry(element))
{ {
int order = op.getQuadratureDegree(geometry.type(), geometry, degree, type); localContext_ = &localContext;
auto const& quad = LocalQuadratureRules::rule(geometry.type(), order); geometry_ = &geometry;
localGeometry_ = &localGeometry;
int order = op.getQuadratureDegree(localGeometry.type(), localGeometry, degree_, type_);
auto const& quad = LocalQuadratureRules::rule(localGeometry.type(), order);
quadrature_.reset(new QuadratureRule(localContext, quad));
}
quadrature.reset(new QuadratureRule(element, quad)); void unbind()
{
localContext_ = nullptr;
geometry_ = nullptr;
localGeometry_ = nullptr;
} }
/// \brief Returns reference to the transformed quadrature rule /// \brief Returns reference to the transformed quadrature rule
...@@ -42,34 +59,44 @@ namespace AMDiS ...@@ -42,34 +59,44 @@ namespace AMDiS
**/ **/
QuadratureRule const& getQuadrature() const QuadratureRule const& getQuadrature() const
{ {
return *quadrature; assert( quadrature_ );
return *quadrature_;
}
LocalContext const& getLocalContext() const
{
assert( localContext_ );
return *localContext_;
} }
/// Return the geometry of the Object
/**
* The geometry is either the element geometry, or the geometry of the
* inside-element for intersections.
**/
Geometry const& getGeometry() const Geometry const& getGeometry() const
{ {
return geometry; assert( geometry_ );
return *geometry_;
} }
protected: LocalGeometry const& getLocalGeometry() const
{
assert( localGeometry_ );
return *localGeometry_;
}
// transform coords from intersection coords to element coords
template <class Coordinate> template <class Coordinate>
decltype(auto) map(Coordinate const& p) const decltype(auto) getLocal(Coordinate const& p) const
{ {
return get_position<LocalContext>(geometry, p); return get_position<LocalContext>(getLocalGeometry(), p);
} }
private: private:
/// transformed quadrature rule int degree_;
std::shared_ptr<QuadratureRule const> quadrature; FirstOrderType type_;
/// geometry() for entities or geometryInInside() for intersections LocalContext const* localContext_ = nullptr;
Geometry geometry; Geometry const* geometry_ = nullptr;
LocalGeometry const* localGeometry_ = nullptr;
/// transformed quadrature rule
std::shared_ptr<QuadratureRule const> quadrature_;
}; };
} // end namespace AMDiS } // end namespace AMDiS
...@@ -27,11 +27,22 @@ namespace AMDiS ...@@ -27,11 +27,22 @@ namespace AMDiS
using Self = Operator; using Self = Operator;
using OperatorTermType = OperatorTerm<LocalContext>; using OperatorTermType = OperatorTerm<LocalContext>;
using Element = typename GridView::template Codim<0>::Entity;
using Geometry = typename Element::Geometry;
using LocalGeometry = typename Impl::Get<LocalContext>::Geometry;
using ElementVector = Impl::ElementVector;
using ElementMatrix = Impl::ElementMatrix;
using QuadratureRule = QuadratureRuleFactory_t<LocalContext, typename GridView::ctype, GridView::dimension>;
template <class,class> friend class ZeroOrderAssembler; template <class,class> friend class ZeroOrderAssembler;
template <class,class, FirstOrderType> friend class FirstOrderAssembler; template <class,class, FirstOrderType> friend class FirstOrderAssembler;
template <class,class> friend class SecondOrderAssembler; template <class,class> friend class SecondOrderAssembler;
public: public:
/// \brief Add coefficients for zero-order operator < C(phi), psi >. /// \brief Add coefficients for zero-order operator < C(phi), psi >.
/** /**
* If `c` is a scalar expression, then C(phi) := c * phi, * If `c` is a scalar expression, then C(phi) := c * phi,
...@@ -113,38 +124,77 @@ namespace AMDiS ...@@ -113,38 +124,77 @@ namespace AMDiS
public: // initialize and assemble operator on element public: // initialize and assemble operator on element
/// Extract the polynomial degree from \p rowBasis and \p colBasis. /// Extract the polynomial degree from \p rowBasis and \p colBasis.
template <class RowBasis, class ColBasis> template <class RowNode, class ColNode>
void init(RowBasis const& rowBasis, void init(RowNode const& rowNode, ColNode const& colNode);
ColBasis const& colBasis);
template <class RowNode>
void init(RowNode const& node);
void bind(Element const& element, Geometry const& geometry)
{
element_ = &element;
geometry_ = &geometry;
}
template <class RowBasis> void unbind()
void init(RowBasis const& rowBasis); {
element_ = nullptr;
geometry_ = nullptr;
}
Element const& getElement() const
{
assert( element_ );
return *element_;
}
Geometry const& getGeometry() const
{
assert( geometry_ );
return *geometry_;
}
decltype(auto) getLocalGeometry(LocalContext const& context) const
{
return getLocalGeometry(context, std::is_same<LocalContext, Element>{});
}
private: // implementation detail
Geometry const& getLocalGeometry(Element const& element, std::true_type) const
{
return *geometry_;
}
LocalGeometry getLocalGeometry(LocalContext const& intersection, std::false_type) const
{
return intersection.geometryInInside();
}
public:
/// Calculates the needed quadrature degree for the given term-order \p order. /// Calculates the needed quadrature degree for the given term-order \p order.
/// This takes into account the polynomial degree of trial and test functions /// This takes into account the polynomial degree of trial and test functions
/// and the polynomial degree of the coefficient functions /// and the polynomial degree of the coefficient functions
template <class Geometry> int getQuadratureDegree(Dune::GeometryType t,
int getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, LocalGeometry const& geometry, int order,
FirstOrderType firstOrderType = GRD_PHI); FirstOrderType firstOrderType = GRD_PHI);
// assemble matrix operator // assemble matrix operator
template <class RowView, class ColView, class ElementMatrix> template <class RowNode, class ColNode>
bool assemble( bool assemble(
LocalContext const& element, LocalContext const& element,
ElementMatrix& elementMatrix, ElementMatrix& elementMatrix,
RowView const& rowView, RowNode const& rowNode,
ColView const& colView, ColNode const& colNode,
Flag const& optimizationFlag,
double* factor = NULL); double* factor = NULL);
// assemble vector operator // assemble vector operator
template <class LocalView, class ElementVector> template <class RowNode>
bool assemble( bool assemble(
LocalContext const& element, LocalContext const& element,
ElementVector& elementVector, ElementVector& elementVector,
LocalView const& localView, RowNode const& node,
double* factor = NULL); double* factor = NULL);
...@@ -168,8 +218,8 @@ namespace AMDiS ...@@ -168,8 +218,8 @@ namespace AMDiS
protected: // sum up constribution from all operator-terms protected: // sum up constribution from all operator-terms
template <class OrderTag, class E, class Quadrature> template <class OrderTag>
void init(OrderTag, E const& element, Quadrature const& quad) void init(OrderTag, LocalContext const& element, QuadratureRule const& quad)
{ {
for (auto* operatorTerm : getOperatorTerms(OrderTag{})) for (auto* operatorTerm : getOperatorTerms(OrderTag{}))
operatorTerm->init(element, quad); operatorTerm->init(element, quad);
...@@ -182,6 +232,7 @@ namespace AMDiS ...@@ -182,6 +232,7 @@ namespace AMDiS
result += operatorTerm->eval(iq, args...) * factor; result += operatorTerm->eval(iq, args...) * factor;
} }
private: private:
std::list<OperatorTermType*>& getOperatorTerms(tag::zot) { return zeroOrder; } std::list<OperatorTermType*>& getOperatorTerms(tag::zot) { return zeroOrder; }
...@@ -189,12 +240,6 @@ namespace AMDiS ...@@ -189,12 +240,6 @@ namespace AMDiS
std::list<OperatorTermType*>& getOperatorTerms(tag::fot_grd_psi) { return firstOrderGrdPsi; } std::list<OperatorTermType*>& getOperatorTerms(tag::fot_grd_psi) { return firstOrderGrdPsi; }
std::list<OperatorTermType*>& getOperatorTerms(tag::sot) { return secondOrder; } std::list<OperatorTermType*>& getOperatorTerms(tag::sot) { return secondOrder; }
template <class LocalView>
std::size_t getElementIndex(GridView const& gridView, LocalView const& localView) const
{
return gridView.indexSet().index(localView.element());
}
private: private:
...@@ -210,11 +255,16 @@ namespace AMDiS ...@@ -210,11 +255,16 @@ namespace AMDiS
/// List of all zero order terms /// List of all zero order terms
std::list<OperatorTermType*> zeroOrder; std::list<OperatorTermType*> zeroOrder;
int psiDegree = 1; int psiDegree_ = 1;
int phiDegree = 1; int phiDegree_ = 1;
Element const* element_ = nullptr;
Geometry const* geometry_ = nullptr;
std::size_t lastMatrixIndex = 0; ZeroOrderAssembler<GridView, LocalContext> zeroOrderAssembler;
std::size_t lastVectorIndex = 0; FirstOrderAssembler<GridView, LocalContext, GRD_PHI> firstOrderGrdPhiAssembler;
FirstOrderAssembler<GridView, LocalContext, GRD_PSI> firstOrderGrdPsiAssembler;
SecondOrderAssembler<GridView, LocalContext> secondOrderAssembler;
}; };
......
...@@ -2,41 +2,28 @@ ...@@ -2,41 +2,28 @@
namespace AMDiS { namespace AMDiS {
template <class GridView, class Element> template <class GridView, class LocalContext>
template <class RowBasis, class ColBasis> template <class RowNode, class ColNode>
void Operator<GridView, Element>:: void Operator<GridView, LocalContext>::init(
init(RowBasis const& rowBasis, ColBasis const& colBasis) RowNode const& rowNode,
ColNode const& colNode)
{ {
lastMatrixIndex = std::numeric_limits<std::size_t>::max(); psiDegree_ = getPolynomialDegree(rowNode);
lastVectorIndex = std::numeric_limits<std::size_t>::max(); phiDegree_ = getPolynomialDegree(colNode);
// auto const& rowFE = rowView.tree().finiteElement();
// auto const& colFE = colView.tree().finiteElement();
// psiDegree = rowFE.localBasis().order();
// phiDegree = colFE.localBasis().order();
psiDegree = getPolynomialDegree<RowBasis>;
phiDegree = getPolynomialDegree<ColBasis>;
// TODO: calc quadrature degree here.
} }
template <class GridView, class Element> template <class GridView, class LocalContext>
template <class RowBasis> template <class RowNode>
void Operator<GridView, Element>:: void Operator<GridView, LocalContext>::init(RowNode const& node)
init(RowBasis const& rowBasis)
{ {
lastVectorIndex = std::numeric_limits<std::size_t>::max(); psiDegree_ = getPolynomialDegree(node);
psiDegree = getPolynomialDegree<RowBasis>;
} }
template <class GridView, class Element> template <class GridView, class LocalContext>
template <class Geometry> int Operator<GridView, LocalContext>::
int Operator<GridView, Element>:: getQuadratureDegree(Dune::GeometryType t, LocalGeometry const& geometry, int order, FirstOrderType type)
getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, FirstOrderType type)
{ {
std::list<OperatorTermType*>* terms = NULL; std::list<OperatorTermType*>* terms = NULL;
...@@ -60,7 +47,7 @@ getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, F ...@@ -60,7 +47,7 @@ getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, F
for (OperatorTermType* term : *terms) for (OperatorTermType* term : *terms)
maxTermDegree = std::max(maxTermDegree, term->getDegree(t)); maxTermDegree = std::max(maxTermDegree, term->getDegree(t));
int degree = psiDegree + phiDegree + maxTermDegree; int degree = psiDegree_ + phiDegree_ + maxTermDegree;
if (t.isSimplex()) if (t.isSimplex())
degree -= order; degree -= order;
...@@ -71,14 +58,13 @@ getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, F ...@@ -71,14 +58,13 @@ getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, F
} }
template <class GridView, class Element> template <class GridView, class LocalContext>
template <class RowView, class ColView, class ElementMatrix> template <class RowNode, class ColNode>
bool Operator<GridView, Element>::assemble( bool Operator<GridView, LocalContext>::assemble(
Element const& element, LocalContext const& context,
ElementMatrix& elementMatrix, ElementMatrix& elementMatrix,
RowView const& rowView, RowNode const& rowNode,
ColView const& colView, ColNode const& colNode,
Flag const& /*optimizationFlag*/,