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

operator-terms for gradient adn partial derivative

parent 1635767b
cmake_minimum_required(VERSION 3.0)
project(dune-amdis CXX)
set(ALBERTA_ROOT /opt/software/alberta)
set(ALBERTA_INCLUDE_DIR /opt/software/alberta/include)
set(CMAKE_PREFIX_PATH /opt/software/dune/lib/cmake)
set(UG_DIR /opt/software/dune/lib/cmake/ug)
set(Vc_DIR /opt/software/dune/lib/cmake/Vc)
# set(ALBERTA_ROOT /opt/software/alberta)
# set(ALBERTA_INCLUDE_DIR /opt/software/alberta/include)
# set(CMAKE_PREFIX_PATH /opt/software/dune/lib/cmake)
# set(UG_DIR /opt/software/dune/lib/cmake/ug)
# set(Vc_DIR /opt/software/dune/lib/cmake/Vc)
set(MTL_DIR /opt/sources/amdis2/lib/mtl4)
set(MTL_DIR ${CMAKE_SOURCE_DIR}/install/share/mtl)
if(NOT (dune-common_DIR OR dune-common_ROOT OR
"${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
......
......@@ -3,6 +3,7 @@
#include <list>
#include "OperatorTerm.hpp"
#include "TermGenerator.hpp"
namespace AMDiS
{
......@@ -19,74 +20,78 @@ namespace AMDiS
using Self = Operator;
using OperatorTermType = OperatorTerm<MeshView>;
public:
template <class RowFeSpace, class ColFeSpace>
void init(RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace);
template <class RowView, class ColView, class ElementMatrix>
bool getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double* factor = NULL);
template <class RowView, class ElementVector>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
double* factor = NULL);
public:
/// Add coefficients for zero-order operator < c * u, v >.
/// The coefficient \p c must be a scalar expression.
template <class Term>
Self& addZOT(Term const& term)
Self& addZOT(Term const& c)
{
zeroOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
return *this;
return addZOTImpl(toTerm(c));
}
/// 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>
Self& addFOT(Term const& term, FirstOrderType firstOrderType)
Self& addFOT(Term const& b, FirstOrderType firstOrderType)
{
using OpTerm = GenericOperatorTerm<MeshView, Term>;
if (firstOrderType == GRD_PHI)
firstOrderGrdPhi.push_back(new OpTerm(term));
else
firstOrderGrdPsi.push_back(new OpTerm(term));
return *this;
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
/// \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>
Self& addFOT(Term const& term, const index_<I>, FirstOrderType firstOrderType)
Self& addFOT(Term const& b, const index_<I> _i, FirstOrderType firstOrderType)
{
using OpTerm = GenericOperatorTerm<MeshView, Term, VectorComponent<I>>;
if (firstOrderType == GRD_PHI)
firstOrderGrdPhi.push_back(new OpTerm(term));
else
firstOrderGrdPsi.push_back(new OpTerm(term));
return *this;
return addFOTImpl(toTerm(b), _i, firstOrderType);
}
/// Add coefficients for second-order operator < grad(psi), A * grad(phi) >,
/// where \p A can be a matrix or a scalar expression.
template <class Term>
Self& addSOT(Term const& term)
Self& addSOT(Term const& A)
{
secondOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
return *this;
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
/// second index \p _j to the derivative component of the trial function.
template <class Term, size_t I, size_t J>
Self& addSOT(Term const& term, const index_<I>, const index_<J>)
Self& addSOT(Term const& A, const index_<I> _i, const index_<J> _j)
{
using OpTerm = GenericOperatorTerm<MeshView, Term, MatrixComponent<I,J>>;
secondOrder.push_back(new OpTerm(term));
return *this;
return addSOTImpl(toTerm(A), _i, _j);
}
/// 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 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
/// and the polynomial degree of the coefficient functions
int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);
template <class RowView, class ColView, class ElementMatrix>
bool getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double* factor = NULL);
template <class RowView, class ElementVector>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
double* factor = NULL);
protected:
template <class RowView, class ColView, class ElementMatrix>
......@@ -112,12 +117,34 @@ namespace AMDiS
ElementMatrix& elementMatrix,
double factor);
template <class RowView, class ElementVector>
void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ElementVector& elementVector,
double factor);
template <class RowView, class ColView, class ElementMatrix>
void assembleSecondOrderTerms(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double factor);
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>,
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>);
private:
/// List of all second order terms
std::list<OperatorTermType*> secondOrder;
......
......@@ -24,6 +24,35 @@ namespace AMDiS
}
template <class MeshView>
int Operator<MeshView>::getQuadratureDegree(int order, FirstOrderType firstOrderType)
{
std::list<OperatorTermType*>* terms = NULL;
switch(order)
{
case 0:
terms = &zeroOrder;
break;
case 1:
if (firstOrderType == 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());
return psiDegree + phiDegree - order + maxTermDegree;
}
template <class MeshView>
template <class RowView, class ColView, class ElementMatrix>
bool Operator<MeshView>::getElementMatrix(RowView const& rowView,
......@@ -64,8 +93,13 @@ namespace AMDiS
if (fac == 0.0)
return false;
AMDIS_TEST_EXIT( firstOrderGrdPhi.empty() && secondOrder.empty(),
"Derivatives on ansatz-functions not allowed on the vector-side!");
if (!zeroOrder.empty())
assembleZeroOrderTerms(rowView, elementVector, fac);
if (!firstOrderGrdPsi.empty())
assembleFirstOrderTermsGrdPsi(rowView, elementVector, fac);
return true;
}
......@@ -286,6 +320,59 @@ namespace AMDiS
}
template <class MeshView>
template <class RowView, class ElementVector>
void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ElementVector& elementvector,
double factor)
{
AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPsi(elementvector)");
using Element = typename RowView::Element;
auto element = rowView.element();
const int dim = Element::dimension;
auto geometry = element.geometry();
const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(0);
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPsi)
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 transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
// The multiplicative factor in the integral transformation formula
const double integrationElement = geometry.integrationElement(quadPos) * factor;
// The gradients of the shape functions on the reference element
std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dim> > rowGradients(rowReferenceGradients.size());
for (size_t i = 0; i < rowGradients.size(); ++i)
jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
for (size_t i = 0; i < elementvector.size(); ++i) {
const int local_i = rowView.tree().localIndex(i);
for (auto* operatorTerm : firstOrderGrdPsi)
elementvector[local_i]
+= operatorTerm->evalFot2(iq, rowGradients[i])
* quad[iq].weight() * integrationElement;
}
}
}
template <class MeshView>
template <class RowView, class ColView, class ElementMatrix>
void Operator<MeshView>::assembleSecondOrderTerms(RowView const& rowView,
......@@ -345,33 +432,65 @@ namespace AMDiS
}
}
template <class MeshView>
int Operator<MeshView>::getQuadratureDegree(int order, FirstOrderType firstOrderType)
template <class Term>
Operator<MeshView>& Operator<MeshView>::addZOTImpl(Term const& term)
{
std::list<OperatorTermType*>* terms = NULL;
switch(order)
{
case 0:
terms = &zeroOrder;
break;
case 1:
if (firstOrderType == 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());
return psiDegree + phiDegree - order + maxTermDegree;
zeroOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
return *this;
}
template <class MeshView>
template <class Term>
Operator<MeshView>& Operator<MeshView>::addFOTImpl(Term const& term,
FirstOrderType firstOrderType)
{
using OpTerm = GenericOperatorTerm<MeshView, Term>;
if (firstOrderType == GRD_PHI)
firstOrderGrdPhi.push_back(new OpTerm(term));
else
firstOrderGrdPsi.push_back(new OpTerm(term));
return *this;
}
template <class MeshView>
template <class Term, size_t I>
Operator<MeshView>& Operator<MeshView>::addFOTImpl(Term const& term,
const index_<I>,
FirstOrderType firstOrderType)
{
using OpTerm = GenericOperatorTerm<MeshView, Term, VectorComponent<I>>;
if (firstOrderType == GRD_PHI)
firstOrderGrdPhi.push_back(new OpTerm(term));
else
firstOrderGrdPsi.push_back(new OpTerm(term));
return *this;
}
template <class MeshView>
template <class Term>
Operator<MeshView>& Operator<MeshView>::addSOTImpl(Term const& term)
{
secondOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
return *this;
}
template <class MeshView>
template <class Term, size_t I, size_t J>
Operator<MeshView>& Operator<MeshView>::addSOTImpl(Term const& term,
const index_<I>,
const index_<J>)
{
using OpTerm = GenericOperatorTerm<MeshView, Term, MatrixComponent<I,J>>;
secondOrder.push_back(new OpTerm(term));
return *this;
}
} // end namespace AMDiS
......@@ -298,6 +298,7 @@ namespace AMDiS
return {std::forward<DOFVectorType>(vector), factor};
}
// -----------------------------------------------------------------------------
......@@ -382,4 +383,177 @@ namespace AMDiS
}
// -----------------------------------------------------------------------------
/// An expression that evalues a DOFVector at a given element, either on the
/// dofs or on the quadrature points
template <class DOFVectorType>
class GradientTerm
{
using Basis = typename DOFVectorType::FeSpace::LocalView::Tree::FiniteElement::Traits::LocalBasisType;
using Jacobian = typename Basis::Traits::JacobianType;
public:
using value_type = Jacobian;
using field_type = typename DOFVectorType::value_type;
GradientTerm(DOFVectorType const& dofvector, double factor = 1.0)
: vector(dofvector.getVector())
, factor(factor)
, localView(dofvector.getFeSpace().localView())
, localIndexSet(dofvector.getFeSpace().localIndexSet())
{}
template <class Element, class PointList>
void init(Element const& element,
PointList const& points)
{
AMDIS_FUNCNAME("DOFVectorTerm::init()");
localView.bind(element);
localIndexSet.bind(localView);
const auto& localFiniteElem = localView.tree().finiteElement();
const size_t nBasisFct = localFiniteElem.size();
std::vector<Jacobian> shapeGradients(nBasisFct);
std::vector<field_type> localVec(nBasisFct);
for (size_t j = 0; j < nBasisFct; ++j) {
const auto global_idx = localIndexSet.index(j);
localVec[j] = factor * vector[global_idx];
}
values.resize(points.size());
for (size_t iq = 0; iq < points.size(); ++iq) {
localFiniteElem.localBasis().evaluateJacobian(points[iq].position(), shapeGradients);
value_type result = 0;
for (size_t j = 0; j < shapeGradients.size(); ++j)
result += localVec[j] * shapeGradients[j];
values[iq] = result;
}
}
value_type const& operator[](size_t iq) const
{
return values[iq];
}
int getDegree() const
{
return math::max(0, degree-1);
}
private:
typename DOFVectorType::BaseVector const& vector;
double factor;
using FeSpace = typename DOFVectorType::FeSpace;
typename FeSpace::LocalView localView;
typename FeSpace::LocalIndexSet localIndexSet;
int degree = GetDegree<FeSpace>::value;
std::vector<value_type> values;
};
/// generator function for \ref DOFVector expressions
template <class DOFVectorType>
GradientTerm<std::decay_t<DOFVectorType>>
gradientOf(DOFVectorType&& vector, double factor = 1.0)
{
return {std::forward<DOFVectorType>(vector), factor};
}
/// An expression that evalues a DOFVector at a given element, either on the
/// dofs or on the quadrature points
template <class DOFVectorType>
class DerivativeTerm
{
using Basis = typename DOFVectorType::FeSpace::LocalView::Tree::FiniteElement::Traits::LocalBasisType;
using Jacobian = typename Basis::Traits::JacobianType;
using Derivative = typename Basis::Traits::RangeType;
public:
using value_type = Derivative;
using field_type = typename DOFVectorType::value_type;
DerivativeTerm(DOFVectorType const& dofvector, int component, double factor = 1.0)
: vector(dofvector.getVector())
, component(component)
, factor(factor)
, localView(dofvector.getFeSpace().localView())
, localIndexSet(dofvector.getFeSpace().localIndexSet())
{
direction[0] = component;
}
template <class Element, class PointList>
void init(Element const& element,
PointList const& points)
{
AMDIS_FUNCNAME("DOFVectorTerm::init()");
localView.bind(element);
localIndexSet.bind(localView);
const auto& localFiniteElem = localView.tree().finiteElement();
const size_t nBasisFct = localFiniteElem.size();
// std::vector<Derivative> shapeDerivative(nBasisFct);
std::vector<Jacobian> shapeJacobian(nBasisFct);
std::vector<field_type> localVec(nBasisFct);
for (size_t j = 0; j < nBasisFct; ++j) {
const auto global_idx = localIndexSet.index(j);
localVec[j] = factor * vector[global_idx];
}
values.resize(points.size());
for (size_t iq = 0; iq < points.size(); ++iq) {
// localFiniteElem.localBasis().template evaluate<1>(direction, points[iq].position(), shapeDerivative);
localFiniteElem.localBasis().evaluateJacobian(points[iq].position(), shapeJacobian);
value_type result = 0;
for (size_t j = 0; j < nBasisFct; ++j)
result += localVec[j] * shapeJacobian[j][0][component];
values[iq] = result;
}
}
value_type const& operator[](size_t iq) const
{
return values[iq];
}
int getDegree() const
{
return math::max(0, degree-1);
}
private:
typename DOFVectorType::BaseVector const& vector;
int component;
double factor;
std::array<int, 1> direction;
using FeSpace = typename DOFVectorType::FeSpace;
typename FeSpace::LocalView localView;
typename FeSpace::LocalIndexSet localIndexSet;
int degree = GetDegree<FeSpace>::value;
std::vector<value_type> values;
};
/// generator function for \ref DOFVector expressions
template <class DOFVectorType>
DerivativeTerm<std::decay_t<DOFVectorType>>
derivativeOf(DOFVectorType&& vector, int component, double factor = 1.0)
{
return {std::forward<DOFVectorType>(vector), component, factor};
}
} // end namespace AMDiS
......@@ -35,12 +35,10 @@ namespace AMDiS
template <class Traits>
void ProblemInstat<Traits>::createUhOld()
{
if (oldSolution)
{
if (oldSolution) {
AMDIS_WARNING("oldSolution already created\n");
}
else
{
else {
const int size = problemStat.getNumComponents();
// create oldSolution
......
......@@ -167,6 +167,7 @@ namespace AMDiS
/// Return the \ref mesh
auto getMesh() { return mesh; }
/// Return the \ref meshView
auto getMeshView() { return meshView; }
/// Return the \ref feSpaces
......@@ -186,6 +187,7 @@ namespace AMDiS
return name;
}
/// Return a vector of names of the problem components
std::vector<std::string> getComponentNames() const
{
return componentNames;
......
/** \file TermGenerator.hpp */
#pragma once
// std c++ includes
#include <type_traits>