Skip to content
GitLab
Commit e653820b authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

partial derivative terms implemented

parent c64c2718
Expand all Hide whitespace changes
Inline Side-by-side
dune/amdis/Operator.hpp
View file @ e653820b
......@@ -3,6 +3,8 @@
#include <list>
#include <vector>
#include <dune/geometry/type.hh>
#include <dune/amdis/OperatorTermBase.hpp>
#include <dune/amdis/terms/TermGenerator.hpp>
......@@ -16,8 +18,8 @@ namespace AMDiS
GRD_PSI,
GRD_PHI
};
template <class MeshView>
class Operator
{
......@@ -26,10 +28,10 @@ namespace AMDiS
using ElementVector = Impl::ElementVector;
using ElementMatrix = Impl::ElementMatrix;
using IdType = typename MeshView::Grid::LocalIdSet::IdType;
public:
public:
/// Add coefficients for zero-order operator < c * u, v >.
/// The coefficient \p c must be a scalar expression.
template <class Term>
......@@ -37,7 +39,7 @@ namespace AMDiS
{
return addZOTImpl(toTerm(c));
}
template <class... Args>
static shared_ptr<Self> zot(Args&&... args)
{
......@@ -45,10 +47,10 @@ namespace AMDiS
op->addZOT(std::forward<Args>(args)...);
return op;
}
/// Add coefficients for first-order operator < psi, term * grad(phi) > or
/// < grad(psi), term * phi >, where the first corresponds to
/// Add coefficients for first-order operator < psi, term * grad(phi) > or
/// < grad(psi), term * phi >, where the first corresponds to
/// \p firstOrderType == GRD_PHI and the second one to \p firstOrderType == GRD_PSI.
/// The coefficient \p b must be a vector expression.
template <class Term>
......@@ -56,9 +58,9 @@ namespace AMDiS
{
return addFOTImpl(toTerm(b), firstOrderType);
}
/// Add coefficients for first-order operator < psi, b * d_i(phi) > or
/// < d_i(psi), b * phi >, where the first corresponds to
/// Add coefficients for first-order operator < psi, b * d_i(phi) > or
/// < d_i(psi), b * phi >, where the first corresponds to
/// \p firstOrderType == GRD_PHI and the second one to \p firstOrderType == GRD_PSI.
/// The coefficient \p b must be a scalar expression.
template <class Term, size_t I>
......@@ -66,7 +68,7 @@ namespace AMDiS
{
return addFOTImpl(toTerm(b), _i, firstOrderType);
}
template <class... Args>
static shared_ptr<Self> fot(Args&&... args)
{
......@@ -74,16 +76,16 @@ namespace AMDiS
op->addFOT(std::forward<Args>(args)...);
return op;
}
/// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
/// where \p A can be a matrix or a scalar expression.
/// where \p A can be a matrix or a scalar expression.
template <class Term>
Self& addSOT(Term const& A)
{
return addSOTImpl(toTerm(A));
}
/// Add coefficients for second-order operator < d_i(psi), A * d_j(phi) >,
/// where \p A can be a matrix or a scalar expression. The first index \p _i
/// corresponds to the derivative component of the test-function and the
......@@ -93,7 +95,7 @@ namespace AMDiS
{
return addSOTImpl(toTerm(A), _i, _j);
}
template <class... Args>
static shared_ptr<Self> sot(Args&&... args)
{
......@@ -101,110 +103,111 @@ namespace AMDiS
op->addSOT(std::forward<Args>(args)...);
return op;
}
/// Calls \ref zot(), \ref for() or \ref sot(), depending on template
/// Calls \ref zot(), \ref for() or \ref sot(), depending on template
/// parameter \p Order.
template <size_t Order, class... Args>
static shared_ptr<Self> create(Args&&... args)
{
return create(index_<Order>{}, std::forward<Args>(args)...);
}
/// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
template <class RowFeSpace, class ColFeSpace>
void init(RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace);
/// Calculates the needed quadrature degree for the given term-order \p order.
/// This takes into account the polynomial degree of trial and test functions
/// and the polynomial degree of the coefficient functions
int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);
template <class Element>
int getQuadratureDegree(Element const& element, int order, FirstOrderType firstOrderType = GRD_PHI);
template <class RowView, class ColView>
bool getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double* factor = NULL);
template <class RowView>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
double* factor = NULL);
protected:
template <class RowView, class ColView>
void assembleZeroOrderTerms(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix);
template <class RowView>
void assembleZeroOrderTerms(RowView const& rowView,
ElementVector& elementVector);
template <class RowView, class ColView>
void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix);
template <class RowView, class ColView>
void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix);
template <class RowView>
void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ElementVector& elementVector);
template <class RowView, class ColView>
void assembleSecondOrderTerms(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix);
private:
template <class Term>
Self& addZOTImpl(Term const& term);
template <class Term>
Self& addFOTImpl(Term const& term, FirstOrderType firstOrderType);
template <class Term, size_t I>
Self& addFOTImpl(Term const& term, const index_<I>,
Self& addFOTImpl(Term const& term, const index_<I>,
FirstOrderType firstOrderType);
template <class Term>
Self& addSOTImpl(Term const& term);
template <class Term, size_t I, size_t J>
Self& addSOTImpl(Term const& term, const index_<I>, const index_<J>);
template <class... Args>
static shared_ptr<Self> create(index_<0>, Args&&... args)
{
return zot(std::forward<Args>(args)...);
}
template <class... Args>
static shared_ptr<Self> create(index_<1>, Args&&... args)
{
return fot(std::forward<Args>(args)...);
}
template <class... Args>
static shared_ptr<Self> create(index_<2>, Args&&... args)
{
return sot(std::forward<Args>(args)...);
}
private:
template <class LocalView>
IdType getElementId(LocalView const& localView);
private:
/// List of all second order terms
std::list<OperatorTermType*> secondOrder;
......@@ -217,17 +220,17 @@ namespace AMDiS
/// List of all zero order terms
std::list<OperatorTermType*> zeroOrder;
int psiDegree = 1;
int phiDegree = 1;
IdType lastMatrixId = 0;
IdType lastVectorId = 0;
ElementMatrix cachedElementMatrix;
ElementVector cachedElementVector;
};
} // end namespace AMDiS
#include "Operator.inc.hpp"
dune/amdis/Operator.inc.hpp
View file @ e653820b
......@@ -8,29 +8,33 @@ namespace AMDiS
ColFeSpace const& colFeSpace)
{
AMDIS_FUNCNAME("Operator::init()");
using IdType = typename Operator<MeshView>::IdType;
lastMatrixId = std::numeric_limits<IdType>::max();
lastVectorId = std::numeric_limits<IdType>::max();
// auto const& rowFE = rowView.tree().finiteElement();
// auto const& colFE = colView.tree().finiteElement();
// psiDegree = rowFE.localBasis().order();
// phiDegree = colFE.localBasis().order();
psiDegree = getPolynomialDegree<RowFeSpace>;
phiDegree = getPolynomialDegree<ColFeSpace>;
// TODO: calc quadrature degree here.
}
template <class MeshView>
int Operator<MeshView>::getQuadratureDegree(int order, FirstOrderType firstOrderType)
template <class Element>
int Operator<MeshView>::getQuadratureDegree(Element const& element, int order, FirstOrderType firstOrderType)
{
std::list<OperatorTermType*>* terms = NULL;
auto const& t = element.type();
auto geometry = element.geometry();
switch(order)
{
case 0:
......@@ -49,24 +53,31 @@ namespace AMDiS
int maxTermDegree = 0;
for (OperatorTermType* term : *terms)
maxTermDegree = std::max(maxTermDegree, term->getDegree());
maxTermDegree = std::max(maxTermDegree, term->getDegree(t));
int degree = psiDegree + phiDegree + maxTermDegree;
if (t.isSimplex())
degree -= order;
if (!geometry.affine())
degree += 1; // oder += (order+1)
return psiDegree + phiDegree - order + maxTermDegree;
return degree;
}
template <class MeshView>
template <class LocalView>
typename Operator<MeshView>::IdType
typename Operator<MeshView>::IdType
Operator<MeshView>::getElementId(LocalView const& localView)
{
auto const& element = localView.element();
auto const& grid = localView.globalBasis().gridView().grid();
auto const& idSet = grid.localIdSet();
return idSet.id(element);
}
template <class MeshView>
template <class RowView, class ColView>
......@@ -74,28 +85,28 @@ namespace AMDiS
ColView const& colView,
ElementMatrix& elementMatrix,
double* factor)
{
{
AMDIS_FUNCNAME("Operator::getElementMatrix()");
double fac = factor ? *factor : 1.0;
if (fac == 0.0)
return false;
const auto nRows = rowView.tree().finiteElement().size();
const auto nCols = colView.tree().finiteElement().size();
auto currentId = getElementId(rowView);
bool useCachedMatrix = (lastMatrixId == currentId
&& num_rows(cachedElementMatrix) == nRows
&& num_rows(cachedElementMatrix) == nRows
&& num_cols(cachedElementMatrix) == nCols);
bool assign = true;
if (!useCachedMatrix) {
cachedElementMatrix.change_dim(nRows, nCols);
set_to_zero(cachedElementMatrix);
if (!zeroOrder.empty())
assembleZeroOrderTerms(rowView, colView, cachedElementMatrix);
if (!firstOrderGrdPhi.empty())
......@@ -104,63 +115,63 @@ namespace AMDiS
assembleFirstOrderTermsGrdPsi(rowView, colView, cachedElementMatrix);
if (!secondOrder.empty())
assembleSecondOrderTerms(rowView, colView, cachedElementMatrix);
assign = !zeroOrder.empty() || !firstOrderGrdPhi.empty() ||
!firstOrderGrdPsi.empty() || !secondOrder.empty();
}
}degree
if (assign)
elementMatrix += fac * cachedElementMatrix;
lastMatrixId = currentId;
return true;
}
template <class MeshView>
template <class RowView>
bool Operator<MeshView>::getElementVector(RowView const& rowView,
ElementVector& elementVector,
double* factor)
{
{
AMDIS_FUNCNAME("Operator::getElementVector()");
double fac = factor ? *factor : 1.0;
if (fac == 0.0)
return false;
AMDIS_TEST_EXIT( firstOrderGrdPhi.empty() && secondOrder.empty(),
"Derivatives on ansatz-functions not allowed on the vector-side!");
const auto nRows = rowView.tree().finiteElement().size();
const auto nRows = rowView.tree().finiteElement().size();
auto currentId = getElementId(rowView);
bool useCachedVector = (lastVectorId == currentId && size(cachedElementVector) == nRows);
bool assign = true;
if (!useCachedVector) {
if (!useCachedVector) {
cachedElementVector.change_dim(nRows);
set_to_zero(cachedElementVector);
if (!zeroOrder.empty())
assembleZeroOrderTerms(rowView, cachedElementVector);
if (!firstOrderGrdPsi.empty())
assembleFirstOrderTermsGrdPsi(rowView, cachedElementVector);
assign = !zeroOrder.empty() || !firstOrderGrdPsi.empty();
}
if (assign)
elementVector += fac * cachedElementVector;
lastVectorId = currentId;
return true;
}
template <class MeshView>
template <class RowView, class ColView>
void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
......@@ -168,31 +179,31 @@ namespace AMDiS
ElementMatrix& elementMatrix)
{
AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementMatrix)");
using Element = typename RowView::Element;
auto element = rowView.element();
const int dim = Element::dimension;
auto geometry = element.geometry();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(0);
int order = getQuadratureDegree(element, 0);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : zeroOrder)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
// The multiplicative factor in the integral transformation formula
const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
std::vector<Dune::FieldVector<double,1> > rowShapeValues, colShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
if (psiDegree == phiDegree)
colShapeValues = rowShapeValues;
else
......@@ -203,44 +214,42 @@ namespace AMDiS
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : zeroOrder)
elementMatrix[local_i][local_j]
elementMatrix[local_i][local_j]
+= operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j]) * factor;
}
}
}
}
template <class MeshView>
template <class RowView>
void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
ElementVector& elementvector)
{
AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementvector)");
using Element = typename RowView::Element;
auto element = rowView.element();
const int dim = Element::dimension;
auto geometry = element.geometry();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(0);
int order = getQuadratureDegree(element, 0);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : zeroOrder)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
// The multiplicative factor in the integral transformation formula
const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
std::vector<Dune::FieldVector<double,1> > rowShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
......@@ -251,7 +260,7 @@ namespace AMDiS
}
}
}
template <class MeshView>
template <class RowView, class ColView>
......@@ -260,23 +269,21 @@ namespace AMDiS
ElementMatrix& elementMatrix)
{
AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPhi(elementMatrix)");
using Element = typename RowView::Element;
auto element = rowView.element();
auto geometry = element.geometry();
const int dim = Element::dimension;
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(2);
int order = getQuadratureDegree(element, 1, GRD_PHI);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPhi)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
......@@ -289,7 +296,7 @@ namespace AMDiS
std::vector<Dune::FieldVector<double,1> > rowShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
// The gradients of the shape functions on the reference element
std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
colLocalBasis.evaluateJacobian(quadPos, colReferenceGradients);
......@@ -305,13 +312,13 @@ namespace AMDiS
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : firstOrderGrdPhi)
elementMatrix[local_i][local_j]
elementMatrix[local_i][local_j]
+= operatorTerm->evalFot1(iq, rowShapeValues[i], colGradients[j]) * factor;
}
}
}
}
template <class MeshView>
template <class RowView, class ColView>
......@@ -320,21 +327,21 @@ namespace AMDiS