Skip to content

GitLab

Commit c3a20421 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

preconditioner-adapter and solveAdvanced

parent 7b196bfd
Hide whitespace changes
Inline Side-by-side
dune/amdis/Basic.hpp
View file @ c3a20421
......@@ -23,7 +23,15 @@ namespace AMDiS
template <class T>
using IdxList = std::map< int, std::list<std::shared_ptr<T> > >;
template <size_t I, class T, class A>
T const& get(std::vector<T,A> const& vec)
{
return vec[I];
}
namespace Impl
{
......@@ -63,6 +71,25 @@ namespace AMDiS
}
};
template <class Tuple>
struct FoldTuples
{
// add arg to repeated constructor argument list
template <size_t... Is, class... Tuples>
static Tuple make(Seq<Is...>, Tuples&&... tuples)
{
return Tuple{make_element(index_<Is>(), std::forward<Tuples>(tuples)...)...};
}
template <size_t I, class... Tuples>
static std::tuple_element_t<I, Tuple> make_element(index_<I>, Tuples&&... tuples)
{
using AMDiS::get;
return std::tuple_element_t<I, Tuple>{get<I>(std::forward<Tuples>(tuples))...};
}
};
} // end namespace Impl
// construct a tuple with each element constructed by the same argument arg.
......@@ -73,6 +100,15 @@ namespace AMDiS
std::forward<Arg>(arg));
}
template <class Tuple, class... Args>
Tuple fold_tuples(Args&&... args)
{
return Impl::FoldTuples<Tuple>::make(MakeSeq_t<std::tuple_size<Tuple>::value>(),
std::forward<Args>(args)...);
}
// -----------
template <template<class> class Base, class Tuple, class Indices> struct MakeTuple;
......
dune/amdis/DOFMatrix.hpp 0 → 100644
View file @ c3a20421
#pragma once
#include <string>
#include <memory>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/istl/bcrsmatrix.hh>
#include "ClonablePtr.hpp"
#include "Log.hpp"
namespace AMDiS
{
template <class RowFeSpaceType, class ColFeSpaceType, class ValueType = Dune::FieldMatrix<double,1,1>>
class DOFMatrix
{
public:
using RowFeSpace = RowFeSpaceType;
using ColFeSpace = ColFeSpaceType;
using BaseMatrix = Dune::BCRSMatrix<ValueType>;
using size_type = typename RowFeSpaceType::size_type;
using field_type = typename ValueType::field_type;
using value_type = ValueType;
/// Constructor. Constructs new BaseVector.
DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name)
: rowFeSpace(rowFeSpace)
, colFeSpace(colFeSpace)
, name(name)
, matrix(ClonablePtr<BaseMatrix>::make())
{}
/// Constructor. Takes pointer of data-vector.
DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name,
BaseMatrix& matrix_ref)
: rowFeSpace(rowFeSpace)
, colFeSpace(colFeSpace)
, name(name)
, matrix(matrix_ref)
{}
/// Return the row-basis \ref rowFeSpace of the matrix
RowFeSpace const& getRowFeSpace() const
{
return rowFeSpace;
}
/// Return the col-basis \ref colFeSpace of the matrix
ColFeSpace const& getColFeSpace() const
{
return colFeSpace;
}
/// Return the data-vector \ref vector
BaseMatrix const& getMatrix() const
{
return *matrix;
}
/// Return the data-vector \ref vector
BaseMatrix& getMatrix()
{
return *matrix;
}
/// Return the size of the \ref feSpace
size_type N() const
{
return rowFeSpace.size();
}
size_type M() const
{
return colFeSpace.size();
}
/// Return the \ref name of this vector
std::string getName() const
{
return name;
}
/// Access the entry \p i of the \ref vector with read-access.
value_type const& operator()(size_type r, size_type c) const
{
AMDIS_TEST_EXIT_DBG( r < matrix->N() && c < matrix->M() ,
"index (" << r << "," << c << ") out of range [0," << matrix->N() << ")x[0," << matrix->M() << ")" );
AMDIS_WARNING("Direct matrix-element access should not be used!");
return (*matrix)[r][c];
}
/// Access the entry \p i of the \ref vector with write-access.
value_type& operator()(size_type r, size_type c)
{
AMDIS_TEST_EXIT_DBG( r < matrix->N() && c < matrix->M() ,
"index (" << r << "," << c << ") out of range [0," << matrix->N() << ")x[0," << matrix->M() << ")" );
AMDIS_WARNING("Direct matrix-element access should not be used!");
return (*matrix)[r][c];
}
private:
RowFeSpace const& rowFeSpace;
ColFeSpace const& colFeSpace;
std::string name;
ClonablePtr<BaseMatrix> matrix;
// friend class declarations
template <class, class>
friend class SystemMatrix;
};
} // end namespace AMDiS
dune/amdis/DOFVector.hpp
View file @ c3a20421
......@@ -73,8 +73,7 @@ namespace AMDiS
/// Resize the \ref vector to the size of the \ref feSpace.
void compress()
{
if (vector->size() != getSize())
vector->resize(getSize());
vector->resize(getSize() / value_type::dimension);
}
......
dune/amdis/DirichletBC.hpp
View file @ c3a20421
......@@ -6,6 +6,9 @@
#include <dune/functions/common/functionconcepts.hh>
#include "Log.hpp"
#include "Traits.hpp"
namespace AMDiS
{
template <class WorldVector>
......@@ -21,22 +24,35 @@ namespace AMDiS
{}
template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
template <class RowFeSpace, class ColFeSpace, class Matrix, class Vector1, class Vector2>
void init(bool apply,
RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace,
DOFMatrix* matrix,
DOFVector* rhs,
DOFVector* solution);
Matrix* matrix,
Vector1* rhs,
Vector2* solution);
template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
template <class RowFeSpace, class ColFeSpace, class Matrix, class Vector1, class Vector2>
void finish(bool apply,
RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace,
DOFMatrix* matrix,
DOFVector* rhs,
DOFVector* solution);
Matrix* matrix,
Vector1* rhs,
Vector2* solution);
protected:
template <class Iter, class Category>
void setIdentity(bool condition, Iter cIt, Category)
{
AMDIS_ERROR_EXIT("Unknown block-type!");
}
template <class Iter>
void setIdentity(bool condition, Iter cIt, _scalar);
template <class Iter>
void setIdentity(bool condition, Iter cIt, _vector);
private:
std::function<bool(WorldVector)> predicate;
......
dune/amdis/DirichletBC.inc.hpp
View file @ c3a20421
......@@ -2,18 +2,17 @@
#include <dune/functions/functionspacebases/interpolate.hh>
#include "Log.hpp"
namespace AMDiS
{
template <class WorldVector>
template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
template <class RowFeSpace, class ColFeSpace, class Matrix, class Vector1, class Vector2>
void DirichletBC<WorldVector>::init(bool apply,
RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace,
DOFMatrix* matrix,
DOFVector* rhs,
DOFVector* solution)
Matrix* matrix,
Vector1* rhs,
Vector2* solution)
{
using Dune::Functions::interpolate;
......@@ -25,17 +24,18 @@ namespace AMDiS
template <class WorldVector>
template <class RowFeSpace, class ColFeSpace, class DOFMatrix, class DOFVector>
template <class RowFeSpace, class ColFeSpace, class Matrix, class Vector1, class Vector2>
void DirichletBC<WorldVector>::finish(bool apply,
RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace,
DOFMatrix* matrix,
DOFVector* rhs,
DOFVector* solution)
Matrix* matrix,
Vector1* rhs,
Vector2* solution)
{
using Dune::Functions::interpolate;
AMDIS_TEST_EXIT( initialized, "Boundary condition not initialized!" );
using _cat = ValueCategory_t<typename Vector2::block_type>;
// loop over the matrix rows
for (size_t i = 0; i < matrix->N(); ++i) {
......@@ -44,7 +44,7 @@ namespace AMDiS
auto cEndIt = (*matrix)[i].end();
// loop over nonzero matrix entries in current row
for (; cIt != cEndIt; ++cIt)
*cIt = (apply && i == cIt.index()) ? 1.0 : 0.0;
setIdentity(apply && i == cIt.index(), cIt, _cat{});
}
}
......@@ -53,5 +53,22 @@ namespace AMDiS
// interpolate(colFeSpace, *solution, values, dirichletNodes);
}
}
template <class WorldVector>
template <class Iter>
void DirichletBC<WorldVector>::setIdentity(bool condition, Iter cIt, _scalar)
{
*cIt = condition ? 1 : 0;
}
template <class WorldVector>
template <class Iter>
void DirichletBC<WorldVector>::setIdentity(bool condition, Iter cIt, _vector)
{
*cIt = condition ? Dune::ScaledIdentityMatrix<double, WorldVector::dimension>{} : 0;
}
} // end namespace AMDiS
dune/amdis/Operator.hpp
View file @ c3a20421
......@@ -6,6 +6,13 @@
namespace AMDiS
{
enum FirstOrderType
{
GRD_PSI,
GRD_PHI
};
template <class MeshView>
class Operator
{
......@@ -33,6 +40,29 @@ namespace AMDiS
return *this;
}
template <class Term>
Self& addFOT(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 Term, size_t I>
Self& addFOT(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 Term>
Self& addSOT(Term const& term)
{
......@@ -40,9 +70,17 @@ namespace AMDiS
return *this;
}
template <class Term, size_t I, size_t J>
Self& addSOT(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;
}
/// Calculates the needed quadrature degree for the given order.
int getQuadratureDegree(int order/*, FirstOrderType firstOrderType = GRD_PHI*/);
int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);
protected:
......@@ -55,6 +93,16 @@ namespace AMDiS
void assembleZeroOrderTerms(RowView const& rowView,
ElementVector& elementVector);
template <class RowView, class ColView, class ElementMatrix>
void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix);
template <class RowView, class ColView, class ElementMatrix>
void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix);
template <class RowView, class ColView, class ElementMatrix>
void assembleSecondOrderTerms(RowView const& rowView,
ColView const& colView,
......
dune/amdis/Operator.inc.hpp
View file @ c3a20421
......@@ -28,8 +28,14 @@ namespace AMDiS
ColView const& colView,
ElementMatrix& elementMatrix)
{
assembleZeroOrderTerms(rowView, colView, elementMatrix);
assembleSecondOrderTerms(rowView, colView, elementMatrix);
if (!zeroOrder.empty())
assembleZeroOrderTerms(rowView, colView, elementMatrix);
if (!firstOrderGrdPhi.empty())
assembleFirstOrderTermsGrdPhi(rowView, colView, elementMatrix);
if (!firstOrderGrdPsi.empty())
assembleFirstOrderTermsGrdPsi(rowView, colView, elementMatrix);
if (!secondOrder.empty())
assembleSecondOrderTerms(rowView, colView, elementMatrix);
}
......@@ -38,7 +44,8 @@ namespace AMDiS
void Operator<MeshView>::getElementVector(RowView const& rowView,
ElementVector& elementVector)
{
assembleZeroOrderTerms(rowView, elementVector);
if (!zeroOrder.empty())
assembleZeroOrderTerms(rowView, elementVector);
}
......@@ -59,9 +66,8 @@ namespace AMDiS
int order = getQuadratureDegree(0);
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
std::vector<double> functionValues(quad.size(), 0.0);
for (auto* operatorTerm : zeroOrder)
operatorTerm->evalAtQPs(element, quad, functionValues);
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
......@@ -82,7 +88,10 @@ namespace AMDiS
for (size_t j = 0; j < elementMatrix.M(); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
elementMatrix[local_i][local_j] += functionValues[iq] * (rowShapeValues[i] * colShapeValues[j]) * quad[iq].weight() * integrationElement;
for (auto* operatorTerm : zeroOrder)
elementMatrix[local_i][local_j]
+= operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j])
* quad[iq].weight() * integrationElement;
}
}
}
......@@ -105,9 +114,8 @@ namespace AMDiS
int order = getQuadratureDegree(0);
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
std::vector<double> functionValues(quad.size(), 0.0);
for (auto* operatorTerm : zeroOrder)
operatorTerm->evalAtQPs(element, quad, functionValues);
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
......@@ -121,7 +129,124 @@ namespace AMDiS
for (size_t i = 0; i < elementvector.size(); ++i) {
const int local_i = rowView.tree().localIndex(i);
elementvector[local_i] += functionValues[iq] * rowShapeValues[i] * quad[iq].weight() * integrationElement;
for (auto* operatorTerm : zeroOrder)
elementvector[local_i]
+= operatorTerm->evalZot(iq, rowShapeValues[i])
* quad[iq].weight() * integrationElement;
}
}
}
template <class MeshView>
template <class RowView, class ColView, class ElementMatrix>
void Operator<MeshView>::assembleFirstOrderTermsGrdPhi(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix)
{
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();
const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(2);
const auto& 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();
// 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);
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);
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dim> > colGradients(colReferenceGradients.size());
for (size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(colReferenceGradients[i][0], colGradients[i]);
for (size_t i = 0; i < elementMatrix.N(); ++i) {
for (size_t j = 0; j < elementMatrix.M(); ++j) {
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]
+= operatorTerm->evalFot1(iq, rowShapeValues[i], colGradients[j])
* quad[iq].weight() * integrationElement;
}
}
}
}
template <class MeshView>
template <class RowView, class ColView, class ElementMatrix>
void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix)
{
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();
const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(2);
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);
// 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]);
std::vector<Dune::FieldVector<double,1> > colShapeValues;
colLocalBasis.evaluateFunction(quadPos, colShapeValues);
for (size_t i = 0; i < elementMatrix.N(); ++i) {
for (size_t j = 0; j < elementMatrix.M(); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : firstOrderGrdPsi)
elementMatrix[local_i][local_j]
+= operatorTerm->evalFot2(iq, rowGradients[i], colShapeValues[j])
* quad[iq].weight() * integrationElement;
}
}
}
}
......@@ -144,9 +269,8 @@ namespace AMDiS
int order = getQuadratureDegree(2);
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
std::vector<double> functionValues(quad.size(), 0.0);
for (auto* operatorTerm : secondOrder)
operatorTerm->evalAtQPs(element, quad, functionValues);
operatorTerm->init(element, quad);
// TODO: currently only the implementation for equal fespaces
assert( psiDegree == phiDegree );
......@@ -175,7 +299,10 @@ namespace AMDiS
for (size_t j = 0; j < elementMatrix.M(); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
elementMatrix[local_i][local_j] += functionValues[iq] * (gradients[i] * gradients[j]) * quad[iq].weight() * integrationElement;
for (auto* operatorTerm : secondOrder)
elementMatrix[local_i][local_j]
+= operatorTerm->evalSot(iq, gradients[i], gradients[j])
* quad[iq].weight() * integrationElement;
}
}
}
......@@ -183,7 +310,7 @@ namespace AMDiS
template <class MeshView>
int Operator<MeshView>::getQuadratureDegree(int order/*, FirstOrderType firstOrderType*/)
int Operator<MeshView>::getQuadratureDegree(int order, FirstOrderType firstOrderType)
{
std::list<OperatorTermType*>* terms = NULL;