diff --git a/doc/doxygen/Doxylocal b/doc/doxygen/Doxylocal
index c222d6b5c15445d47c50ab68cfca74d9a089e70b..46c6ffc56e3feacfe289982244995b54ae226876 100644
--- a/doc/doxygen/Doxylocal
+++ b/doc/doxygen/Doxylocal
@@ -1,27 +1,50 @@
# This file contains local changes to the doxygen configuration
# please us '+=' to add file/directories to the lists
+FILE_PATTERNS += *.hpp \
+ *.cpp
+
+HIDE_SCOPE_NAMES = YES
+HIDE_UNDOC_CLASSES = NO
+INTERNAL_DOCS = NO
+MARKDOWN_SUPPORT = YES
+
+EXCLUDE_SYMBOLS = AMDiS::Impl \
+ AMDiS::traits::Impl \
+ AMDiS::detail \
+ itl::details
+
+PREDEFINED += HAVE_UMFPACK \
+ HAVE_ALBERTA \
+ HAVE_UG \
+ AMDIS_BACKEND_MTL
+
# The INPUT tag can be used to specify the files and/or directories that contain
# documented source files. You may enter file names like "myfile.cpp" or
# directories like "/usr/src/myproject". Separate the files or directories
# with spaces.
-INPUT += @top_srcdir@/dune/
+INPUT += @top_srcdir@/dune/amdis \
+ @top_srcdir@/dune/amdis/common \
+ @top_srcdir@/dune/amdis/utility \
+ @top_srcdir@/dune/amdis/linear_algebra \
+ @top_srcdir@/dune/amdis/linear_algebra/mtl
# see e.g. dune-grid for the examples of mainpage and modules
-# INPUT += @srcdir@/mainpage \
-# @srcdir@/modules
+#INPUT += @srcdir@/mainpage \
+# @srcdir@/modules
# The EXCLUDE tag can be used to specify files and/or directories that should
# excluded from the INPUT source files. This way you can easily exclude a
# subdirectory from a directory tree whose root is specified with the INPUT tag.
-# EXCLUDE += @top_srcdir@/dune/amdis/test
+EXCLUDE += @top_srcdir@/dune/amdis/test \
+ @top_srcdir@/dune/amdis/linear_algebra/istl
# The EXAMPLE_PATH tag can be used to specify one or more files or
# directories that contain example code fragments that are included (see
# the \include command).
-# EXAMPLE_PATH += @top_srcdir@/src
+EXAMPLE_PATH += @top_srcdir@/src
# The IMAGE_PATH tag can be used to specify one or more files or
# directories that contain image that are included in the documentation (see
diff --git a/dune-amdis.pc.in b/dune-amdis.pc.in
index 81ec962ed91025de68e4548e5114de4f406781e2..03b2511183256c5e3590b1d13ce0771a9ee8f9a0 100644
--- a/dune-amdis.pc.in
+++ b/dune-amdis.pc.in
@@ -9,7 +9,7 @@ DEPENDENCIES=@REQUIRES@
Name: @PACKAGE_NAME@
Version: @VERSION@
Description: amdis module
-URL: http://dune-project.org/
+URL: http://www.github.com/spraetor
Requires: dune-common dune-geometry dune-localfunctions dune-istl dune-typetree dune-grid dune-functions
Libs: -L${libdir}
Cflags: -I${includedir}
diff --git a/dune/amdis/DirichletBC.hpp b/dune/amdis/DirichletBC.hpp
index 2769429469a117be547e2472a4495a36db1599f1..33148455794d132e9b6cc581b81580cc463c59fb 100644
--- a/dune/amdis/DirichletBC.hpp
+++ b/dune/amdis/DirichletBC.hpp
@@ -10,13 +10,28 @@
namespace AMDiS
{
+ /// Implements a boundary condition of Dirichlet-type.
+ /**
+ * By calling the methods \ref init() and \ref finish before and after
+ * assembling the system-matrix, respectively, dirichlet boundary conditions
+ * can be applied to the matrix and system vector. Therefore a predicate
+ * functions indicates the DOFs where values should be enforced and a second
+ * functor provided in the constructor is responsible for determining the
+ * values to be set at the DOFs.
+ *
+ * In the \ref finish method the matrix is called with \ref applyDirichletBC
+ * to erase the corresponding rows and columns for the DOF indices. This
+ * application of boundary conditions can be symmetric if the matrix does
+ * support this symmetric modification. As a result, this method returns a list
+ * of columns values, that should be subtracted from the rhs.
+ **/
template <class WorldVector>
class DirichletBC
{
public:
template <class Predicate, class Values,
- class = std::enable_if_t< concepts::Functor<Predicate, bool(WorldVector)> &&
- concepts::Functor<Values, double(WorldVector)> > >
+ class = std::enable_if_t< Concepts::Functor<Predicate, bool(WorldVector)> &&
+ Concepts::Functor<Values, double(WorldVector)> > >
DirichletBC(Predicate&& predicate, Values&& values)
: predicate(std::forward<Predicate>(predicate))
, values(std::forward<Values>(values))
diff --git a/dune/amdis/Math.hpp b/dune/amdis/Math.hpp
index bc6e353fbec7717321137f9d1a5e167d667a8eb5..962571e1594f8d51138a8a0805dd700fff8d6f31 100644
--- a/dune/amdis/Math.hpp
+++ b/dune/amdis/Math.hpp
@@ -1,6 +1,3 @@
-/** \defgroup Common Common
- */
-
#pragma once
#include <string>
@@ -8,26 +5,40 @@
#include <cmath>
#include <cfloat>
+#include <boost/math/special_functions/pow.hpp>
+
+#include <dune/amdis/common/ScalarTypes.hpp>
+
namespace AMDiS
{
- namespace math
+ namespace Math
{
- template <class T>
- std::enable_if_t<std::is_arithmetic<T>::value, T>
- constexpr abs(T a)
+ /// Implementation of the absolute value \f$|a|\f$ of arithmetic types.
+ template <class T,
+ class = std::enable_if_t<Concepts::Arithmetic<T>> >
+ constexpr T abs(T a)
{
return a >= 0 ? a : -a;
}
- template <class T>
- std::enable_if_t<std::is_arithmetic<T>::value, T>
- constexpr sqr(T a)
+ /// Implementation of the square \f$ a^2 \f$ of arithmetic types.
+ template <class T,
+ class = std::enable_if_t<Concepts::Arithmetic<T>> >
+ constexpr auto sqr(T a)
{
return a*a;
}
+ template <size_t p, class T,
+ class = std::enable_if_t<Concepts::Arithmetic<T>> >
+ constexpr auto pow(T v)
+ {
+ return boost::math::pow<p>(v);
+ }
+ /// Implementation of the minimum of two values \f$ min(a,b)\f$ of any type
+ /// supporting the `>` relation.
template <class T0, class T1>
constexpr auto min(T0 a, T1 b)
{
@@ -35,13 +46,15 @@ namespace AMDiS
}
+ /// Implementation of the maximum of two values \f$ max(a,b)\f$ of any type
+ /// supporting the `>` relation.
template <class T0, class T1>
constexpr auto max(T0 a, T1 b)
{
return a > b ? a : b;
}
- } // end namespace math
+ } // end namespace Math
template <class T> inline void nullify(T& a)
diff --git a/dune/amdis/Mesh.hpp b/dune/amdis/Mesh.hpp
index 53fe71de4ccbfc46da0710f5326763aed57834b0..2b7e3846688c9c9e51c625fe5b4cef1a64809a3e 100644
--- a/dune/amdis/Mesh.hpp
+++ b/dune/amdis/Mesh.hpp
@@ -129,6 +129,9 @@ namespace AMDiS
} else {
AMDIS_ERROR_EXIT("Construction of UGGrid without filename not yet implemented!");
}
+
+ AMDIS_ERROR_EXIT("No way to construct UG-Grid found");
+ return {};
}
};
#endif
diff --git a/dune/amdis/Operator.hpp b/dune/amdis/Operator.hpp
index ba8faf9c19a9b763b574b8b4c6387e40add1f7c8..64249eff44ac6d1f3629698c8237f287ebc24cdd 100644
--- a/dune/amdis/Operator.hpp
+++ b/dune/amdis/Operator.hpp
@@ -1,9 +1,13 @@
#pragma once
#include <list>
+#include <vector>
-#include "OperatorTerm.hpp"
-#include "TermGenerator.hpp"
+#include <dune/amdis/OperatorTermBase.hpp>
+#include <dune/amdis/terms/TermGenerator.hpp>
+
+#include <dune/amdis/common/TypeDefs.hpp>
+#include <dune/amdis/utility/GetDegree.hpp>
namespace AMDiS
{
@@ -19,6 +23,11 @@ namespace AMDiS
{
using Self = Operator;
using OperatorTermType = OperatorTerm<MeshView>;
+
+ using ElementVector = Impl::ElementVector;
+ using ElementMatrix = Impl::ElementMatrix;
+
+ using IdType = typename MeshView::Grid::LocalIdSet::IdType;
public:
/// Add coefficients for zero-order operator < c * u, v >.
@@ -115,52 +124,46 @@ namespace AMDiS
int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);
- template <class RowView, class ColView, class ElementMatrix>
+ template <class RowView, class ColView>
bool getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double* factor = NULL);
- template <class RowView, class ElementVector>
+ template <class RowView>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
double* factor = NULL);
protected:
- template <class RowView, class ColView, class ElementMatrix>
+ template <class RowView, class ColView>
void assembleZeroOrderTerms(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix,
- double factor);
+ ElementMatrix& elementMatrix);
- template <class RowView, class ElementVector>
+ template <class RowView>
void assembleZeroOrderTerms(RowView const& rowView,
- ElementVector& elementVector,
- double factor);
+ ElementVector& elementVector);
- template <class RowView, class ColView, class ElementMatrix>
+ template <class RowView, class ColView>
void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix,
- double factor);
+ ElementMatrix& elementMatrix);
- template <class RowView, class ColView, class ElementMatrix>
+ template <class RowView, class ColView>
void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix,
- double factor);
+ ElementMatrix& elementMatrix);
- template <class RowView, class ElementVector>
+ template <class RowView>
void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
- ElementVector& elementVector,
- double factor);
+ ElementVector& elementVector);
- template <class RowView, class ColView, class ElementMatrix>
+ template <class RowView, class ColView>
void assembleSecondOrderTerms(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix,
- double factor);
+ ElementMatrix& elementMatrix);
private:
template <class Term>
@@ -179,6 +182,7 @@ namespace AMDiS
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)
{
@@ -197,6 +201,9 @@ namespace AMDiS
return sot(std::forward<Args>(args)...);
}
+ private:
+ template <class LocalView>
+ IdType getElementId(LocalView const& localView);
private:
/// List of all second order terms
@@ -213,6 +220,12 @@ namespace AMDiS
int psiDegree = 1;
int phiDegree = 1;
+
+ IdType lastMatrixId = 0;
+ IdType lastVectorId = 0;
+
+ ElementMatrix cachedElementMatrix;
+ ElementVector cachedElementVector;
};
} // end namespace AMDiS
diff --git a/dune/amdis/Operator.inc.hpp b/dune/amdis/Operator.inc.hpp
index 5e6859f4f2425c34ff215cf872b3307de5fa3aa3..e7d4f8cc03545053fb05fe1d5fbedeb88a31423d 100644
--- a/dune/amdis/Operator.inc.hpp
+++ b/dune/amdis/Operator.inc.hpp
@@ -1,7 +1,5 @@
#pragma once
-#include <vector>
-
namespace AMDiS
{
template <class MeshView>
@@ -11,14 +9,18 @@ namespace AMDiS
{
AMDIS_FUNCNAME("Operator::init()");
-// const auto& rowFE = rowView.tree().finiteElement();
-// const auto& colFE = colView.tree().finiteElement();
+ 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 = GetDegree<RowFeSpace>::value;
- phiDegree = GetDegree<ColFeSpace>::value;
+ psiDegree = getPolynomialDegree<RowFeSpace>;
+ phiDegree = getPolynomialDegree<ColFeSpace>;
// TODO: calc quadrature degree here.
}
@@ -52,9 +54,22 @@ namespace AMDiS
return psiDegree + phiDegree - order + maxTermDegree;
}
+
+ template <class MeshView>
+ template <class LocalView>
+ 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, class ElementMatrix>
+ template <class RowView, class ColView>
bool Operator<MeshView>::getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
@@ -67,21 +82,44 @@ namespace AMDiS
if (fac == 0.0)
return false;
- if (!zeroOrder.empty())
- assembleZeroOrderTerms(rowView, colView, elementMatrix, fac);
- if (!firstOrderGrdPhi.empty())
- assembleFirstOrderTermsGrdPhi(rowView, colView, elementMatrix, fac);
- if (!firstOrderGrdPsi.empty())
- assembleFirstOrderTermsGrdPsi(rowView, colView, elementMatrix, fac);
- if (!secondOrder.empty())
- assembleSecondOrderTerms(rowView, colView, elementMatrix, fac);
+ 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_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())
+ assembleFirstOrderTermsGrdPhi(rowView, colView, cachedElementMatrix);
+ if (!firstOrderGrdPsi.empty())
+ assembleFirstOrderTermsGrdPsi(rowView, colView, cachedElementMatrix);
+ if (!secondOrder.empty())
+ assembleSecondOrderTerms(rowView, colView, cachedElementMatrix);
+
+ assign = !zeroOrder.empty() || !firstOrderGrdPhi.empty() ||
+ !firstOrderGrdPsi.empty() || !secondOrder.empty();
+ }
+
+ if (assign)
+ elementMatrix += fac * cachedElementMatrix;
+
+ lastMatrixId = currentId;
return true;
}
template <class MeshView>
- template <class RowView, class ElementVector>
+ template <class RowView>
bool Operator<MeshView>::getElementVector(RowView const& rowView,
ElementVector& elementVector,
double* factor)
@@ -96,21 +134,38 @@ namespace AMDiS
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);
+ const auto nRows = rowView.tree().finiteElement().size();
+
+ auto currentId = getElementId(rowView);
+ bool useCachedVector = (lastVectorId == currentId && size(cachedElementVector) == nRows);
+
+ bool assign = true;
+ 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, class ElementMatrix>
+ template <class RowView, class ColView>
void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix,
- double factor)
+ ElementMatrix& elementMatrix)
{
AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementMatrix)");
@@ -119,11 +174,11 @@ namespace AMDiS
const int dim = Element::dimension;
auto geometry = element.geometry();
- const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
- const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
+ auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+ auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(0);
- const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
+ auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : zeroOrder)
operatorTerm->init(element, quad);
@@ -133,7 +188,7 @@ namespace AMDiS
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
// The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos) * factor;
+ const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
std::vector<Dune::FieldVector<double,1> > rowShapeValues, colShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
@@ -143,14 +198,13 @@ namespace AMDiS
else
colLocalBasis.evaluateFunction(quadPos, colShapeValues);
- for (size_t i = 0; i < elementMatrix.N(); ++i) {
- for (size_t j = 0; j < elementMatrix.M(); ++j) {
+ for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
+ for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
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]
- += operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j])
- * quad[iq].weight() * integrationElement;
+ += operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j]) * factor;
}
}
}
@@ -159,10 +213,9 @@ namespace AMDiS
template <class MeshView>
- template <class RowView, class ElementVector>
+ template <class RowView>
void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
- ElementVector& elementvector,
- double factor)
+ ElementVector& elementvector)
{
AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementvector)");
@@ -171,41 +224,40 @@ namespace AMDiS
const int dim = Element::dimension;
auto geometry = element.geometry();
- const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+ auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(0);
- const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
+ 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 integrationElement = geometry.integrationElement(quadPos) * factor;
+ const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
std::vector<Dune::FieldVector<double,1> > rowShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
- for (size_t i = 0; i < elementvector.size(); ++i) {
+ for (size_t i = 0; i < size(elementvector); ++i) {
const int local_i = rowView.tree().localIndex(i);
for (auto* operatorTerm : zeroOrder)
- elementvector[local_i]
- += operatorTerm->evalZot(iq, rowShapeValues[i])
- * quad[iq].weight() * integrationElement;
+ elementvector[local_i] += operatorTerm->evalZot(iq, rowShapeValues[i]) * factor;
}
}
}
template <class MeshView>
- template <class RowView, class ColView, class ElementMatrix>
+ template <class RowView, class ColView>
void Operator<MeshView>::assembleFirstOrderTermsGrdPhi(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix,
- double factor)
+ ElementMatrix& elementMatrix)
{
AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPhi(elementMatrix)");
@@ -214,15 +266,17 @@ namespace AMDiS
auto geometry = element.geometry();
const int dim = Element::dimension;
- const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
- const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
+ auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+ auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(2);
- const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
+ 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();
@@ -231,7 +285,7 @@ namespace AMDiS
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
// The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos) * factor;
+ const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
std::vector<Dune::FieldVector<double,1> > rowShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
@@ -246,14 +300,13 @@ namespace AMDiS
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) {
+ for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
+ for (size_t j = 0; j < num_cols(elementMatrix); ++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;
+ += operatorTerm->evalFot1(iq, rowShapeValues[i], colGradients[j]) * factor;
}
}
}
@@ -261,11 +314,10 @@ namespace AMDiS
template <class MeshView>
- template <class RowView, class ColView, class ElementMatrix>
+ template <class RowView, class ColView>
void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix,
- double factor)
+ ElementMatrix& elementMatrix)
{
AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPsi(elementMatrix)");
@@ -274,11 +326,11 @@ namespace AMDiS
auto geometry = element.geometry();
const int dim = Element::dimension;
- const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
- const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
+ auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+ auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(2);
- const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
+ auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPsi)
operatorTerm->init(element, quad);
@@ -291,7 +343,7 @@ namespace AMDiS
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
// The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos) * factor;
+ const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
// The gradients of the shape functions on the reference element
std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
@@ -306,14 +358,13 @@ namespace AMDiS
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) {
+ for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
+ for (size_t j = 0; j < num_cols(elementMatrix); ++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;
+ += operatorTerm->evalFot2(iq, rowGradients[i], colShapeValues[j]) * factor;
}
}
}
@@ -322,10 +373,9 @@ namespace AMDiS
template <class MeshView>
- template <class RowView, class ElementVector>
+ template <class RowView>
void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
- ElementVector& elementvector,
- double factor)
+ ElementVector& elementvector)
{
AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPsi(elementvector)");
@@ -334,10 +384,10 @@ namespace AMDiS
const int dim = Element::dimension;
auto geometry = element.geometry();
- const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+ auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(0);
- const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
+ auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPsi)
operatorTerm->init(element, quad);
@@ -350,7 +400,7 @@ namespace AMDiS
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
// The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos) * factor;
+ const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
// The gradients of the shape functions on the reference element
std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
@@ -362,23 +412,20 @@ namespace AMDiS
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) {
+ for (size_t i = 0; i < size(elementvector); ++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;
+ elementvector[local_i] += operatorTerm->evalFot2(iq, rowGradients[i]) * factor;
}
}
}
template <class MeshView>
- template <class RowView, class ColView, class ElementMatrix>
+ template <class RowView, class ColView>
void Operator<MeshView>::assembleSecondOrderTerms(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix,
- double factor)
+ ElementMatrix& elementMatrix)
{
AMDIS_FUNCNAME("Operator::assembleSecondOrderTerms(elementMatrix)");
@@ -387,11 +434,11 @@ namespace AMDiS
const int dim = Element::dimension;
auto geometry = element.geometry();
- const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
-// const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
+ auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+// auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(2);
- const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
+ auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : secondOrder)
operatorTerm->init(element, quad);
@@ -407,7 +454,7 @@ namespace AMDiS
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
// The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos) * factor;
+ const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
// The gradients of the shape functions on the reference element
std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
@@ -419,14 +466,13 @@ namespace AMDiS
for (size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(referenceGradients[i][0], gradients[i]);
- for (size_t i = 0; i < elementMatrix.N(); ++i) {
- for (size_t j = 0; j < elementMatrix.M(); ++j) {
+ for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
+ for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : secondOrder)
elementMatrix[local_i][local_j]
- += operatorTerm->evalSot(iq, gradients[i], gradients[j])
- * quad[iq].weight() * integrationElement;
+ += operatorTerm->evalSot(iq, gradients[i], gradients[j]) * factor;
}
}
}
diff --git a/dune/amdis/OperatorTermBase.hpp b/dune/amdis/OperatorTermBase.hpp
index e6c8b3a4d32c367e05b3af4ff64d1110238582ac..6d12a92ac1eefea1287501e0a893cd31f387ec6c 100644
--- a/dune/amdis/OperatorTermBase.hpp
+++ b/dune/amdis/OperatorTermBase.hpp
@@ -1,5 +1,9 @@
#pragma once
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#include <vector>
#include <type_traits>
@@ -11,6 +15,7 @@
namespace AMDiS
{
+ /// Abstract base class of all operator terms
template <class MeshView>
class OperatorTerm
{
diff --git a/dune/amdis/ProblemStat.hpp b/dune/amdis/ProblemStat.hpp
index 0bddf047dbf511fcf97c5e66e0926531ecaf08d2..a8147fcb99c3f949a9214c5584e81fcc6cef422b 100644
--- a/dune/amdis/ProblemStat.hpp
+++ b/dune/amdis/ProblemStat.hpp
@@ -1,5 +1,9 @@
#pragma once
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#include <tuple>
#include <string>
#include <memory>
@@ -19,10 +23,12 @@
#include <dune/amdis/Mesh.hpp>
#include <dune/amdis/Operator.hpp>
#include <dune/amdis/ProblemStatBase.hpp>
+#include <dune/amdis/ProblemStatTraits.hpp>
#include <dune/amdis/StandardProblemIteration.hpp>
#include <dune/amdis/common/Timer.hpp>
#include <dune/amdis/common/TupleUtility.hpp>
+#include <dune/amdis/common/TypeDefs.hpp>
#include <dune/amdis/common/Utility.hpp>
namespace AMDiS
@@ -42,8 +48,8 @@ namespace AMDiS
using Codim0 = typename MeshView::template Codim<0>;
using Element = typename Codim0::Entity;
- using ElementVector = Dune::BlockVector<Dune::FieldVector<double,1> >;
- using ElementMatrix = Dune::Matrix<Dune::FieldMatrix<double,1,1> >;
+ using ElementVector = Impl::ElementVector;
+ using ElementMatrix = Impl::ElementMatrix;
/// Dimension of the mesh
@@ -53,7 +59,7 @@ namespace AMDiS
static constexpr int dow = Traits::dimworld;
/// Number of problem components
- static constexpr size_t nComponents = Traits::nComponents;
+ static constexpr int nComponents = Traits::nComponents;
template <size_t I>
@@ -96,7 +102,8 @@ namespace AMDiS
Flag adoptFlag = INIT_NOTHING);
- /// Adds an operator to \ref A.
+ /// Adds an operator to \ref A.
+ /** @{ */
void addMatrixOperator(OperatorType& op,
int i, int j,
double* factor = NULL,
@@ -107,9 +114,15 @@ namespace AMDiS
double* factor = NULL,
double* estFactor = NULL);
- void addMatrixOperator(std::map< std::pair<int,int>, shared_ptr<OperatorType> > ops);
+ void addMatrixOperator(std::map< std::pair<int,int>,
+ shared_ptr<OperatorType> > ops);
+ void addMatrixOperator(std::map< std::pair<int,int>,
+ std::pair<shared_ptr<OperatorType>, bool> > ops);
+ /** @} */
- /// Adds an operator to \ref rhs.
+
+ /// Adds an operator to \ref rhs.
+ /** @{ */
void addVectorOperator(OperatorType& op,
int i,
double* factor = NULL,
@@ -120,8 +133,11 @@ namespace AMDiS
double* factor = NULL,
double* estFactor = NULL);
- void addVectorOperator(std::map< int, shared_ptr<OperatorType> > ops);
-
+ void addVectorOperator(std::map< int,
+ shared_ptr<OperatorType> > ops);
+ void addVectorOperator(std::map< int,
+ std::pair<shared_ptr<OperatorType>, bool> > ops);
+ /** @} */
/// Adds a Dirichlet boundary condition
template <class Predicate, class Values>
@@ -149,7 +165,7 @@ namespace AMDiS
public: // get-methods
/// Returns nr of components \ref nComponents
- static constexpr size_t getNumComponents()
+ static constexpr int getNumComponents()
{
return nComponents;
}
@@ -264,8 +280,8 @@ namespace AMDiS
template <class Matrix, class Vector>
void assemble(std::pair<int, int> row_col,
- Matrix* matrix,
- Vector* rhs);
+ Matrix& matrix, bool asmMatrix,
+ Vector& rhs, bool asmVector);
template <class RowView, class ColView>
bool getElementMatrix(RowView const& rowView,
@@ -357,11 +373,15 @@ namespace AMDiS
/// each matrix block
MatrixEntries<shared_ptr<OperatorType>> matrixOperators;
MatrixEntries<double*> matrixFactors;
+ std::map< std::pair<int,int>, bool > matrixAssembled; // if false, do reassemble
+ std::map< std::pair<int,int>, bool > matrixChanging; // if true, or vectorAssembled false, do reassemble
/// A map (i) -> list<OperatorType> string the differential operators for
/// each rhs block
VectorEntries<shared_ptr<OperatorType>> vectorOperators;
VectorEntries<double*> vectorFactors;
+ std::map< int, bool > vectorAssembled; // if false, do reassemble
+ std::map< int, bool > vectorChanging; // if true, or vectorAssembled false, do reassemble
};
diff --git a/dune/amdis/ProblemStat.inc.hpp b/dune/amdis/ProblemStat.inc.hpp
index 339955c17e1b7bf5ff17a5a720c0b487588e4f47..73274b8e8d0045d685e5eb8023aec37fce46cbbe 100644
--- a/dune/amdis/ProblemStat.inc.hpp
+++ b/dune/amdis/ProblemStat.inc.hpp
@@ -117,6 +117,7 @@ namespace AMDiS
matrixOperators[{i,j}].push_back(op);
matrixFactors[{i,j}].push_back(factor);
+ matrixChanging[{i,j}] = true;
}
@@ -135,6 +136,17 @@ namespace AMDiS
for (auto op : ops)
addMatrixOperator(op.second, op.first.first, op.first.second);
}
+
+ template <class Traits>
+ void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::pair<shared_ptr<OperatorType>,bool> > ops)
+ {
+ for (auto op : ops) {
+ auto row_col = op.first;
+ matrixOperators[row_col].push_back(op.second.first);
+ matrixFactors[row_col].push_back(NULL);
+ matrixChanging[row_col] = matrixChanging[row_col] || op.second.second;
+ }
+ }
template <class Traits>
void ProblemStatSeq<Traits>::addVectorOperator(shared_ptr<OperatorType> op,
@@ -147,6 +159,7 @@ namespace AMDiS
vectorOperators[i].push_back(op);
vectorFactors[i].push_back(factor);
+ vectorChanging[i] = true;
}
template <class Traits>
@@ -157,7 +170,6 @@ namespace AMDiS
{
addVectorOperator(make_shared<OperatorType>(op), i, factor, estFactor);
}
-
template <class Traits>
void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, shared_ptr<OperatorType> > ops)
@@ -166,6 +178,16 @@ namespace AMDiS
addVectorOperator(op.second, op.first);
}
+ template <class Traits>
+ void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::pair<shared_ptr<OperatorType>, bool> > ops)
+ {
+ for (auto op : ops) {
+ vectorOperators[op.first].push_back(op.second.first);
+ vectorFactors[op.first].push_back(NULL);
+ vectorChanging[op.first] = vectorChanging[op.first] || op.second.second;
+ }
+ }
+
/// Adds a Dirichlet boundary condition
template <class Traits>
template <class Predicate, class Values>
@@ -173,10 +195,10 @@ namespace AMDiS
int row, int col,
Values const& values)
{
- static_assert(Dune::Functions::Concept::isCallable<Predicate, WorldVector>(),
+ static_assert( Concepts::Callable<Predicate, WorldVector>,
"Function passed to addDirichletBC for predicate does not model the Callable<WorldVector> concept");
- static_assert(Dune::Functions::Concept::isCallable<Values, WorldVector>(),
+ static_assert( Concepts::Callable<Values, WorldVector>,
"Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");
AMDIS_TEST_EXIT(row >= 0 && row < nComponents, "row number out of range: " << row);
@@ -232,7 +254,7 @@ namespace AMDiS
template <class Traits>
void ProblemStatSeq<Traits>::buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag,
- bool asmMatrix, bool asmVector)
+ bool asmMatrix_, bool asmVector_)
{
Timer t;
@@ -242,25 +264,29 @@ namespace AMDiS
size_t nnz = 0;
- For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector](auto const _r)
+ For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector_](auto const _r)
{
AMDIS_MSG(this->getFeSpace(_r).size() << " DOFs for FeSpace[" << _r << "]");
+ bool asmVector = asmVector_ && (!vectorAssembled[int(_r)] || vectorChanging[int(_r)]);
+
if (asmVector) {
rhs->compress(_r);
rhs->getVector(_r) = 0.0;
}
- For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector, _r](auto const _c)
+ For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector, _r](auto const _c)
{
// The DOFMatrix which should be assembled
auto& dofmatrix = (*systemMatrix)(_r, _c);
auto& solution_vec = (*solution)[_c];
auto& rhs_vec = (*rhs)[_r];
+ auto row_col = std::make_pair(int(_r), int(_c));
+ bool asmMatrix = asmMatrix_ && (!matrixAssembled[row_col] || matrixChanging[row_col]);
+
int r = 0, c = 0;
if (asmMatrix) {
-
// init boundary condition
for (auto bc_list : dirichletBc) {
std::tie(r, c) = bc_list.first;
@@ -269,10 +295,12 @@ namespace AMDiS
bc->init(c == int(_c), dofmatrix, solution_vec, rhs_vec);
}
}
+ }
- // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
- this->assemble({int(_r), int(_c)}, &dofmatrix, (_r == _c && asmVector ? &rhs_vec : NULL));
+ // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
+ this->assemble(row_col, dofmatrix, asmMatrix, rhs_vec, (_r == _c && asmVector));
+ if (asmMatrix) {
// finish boundary condition
for (auto bc_list : dirichletBc) {
std::tie(r, c) = bc_list.first;
@@ -305,53 +333,69 @@ namespace AMDiS
template <class Traits>
template <class Matrix, class Vector>
void ProblemStatSeq<Traits>::assemble(std::pair<int, int> row_col,
- Matrix* dofmatrix,
- Vector* rhs)
- {
- auto const& rowFeSpace = dofmatrix->getRowFeSpace();
- auto const& colFeSpace = dofmatrix->getColFeSpace();
+ Matrix& dofmatrix, bool asmMatrix,
+ Vector& rhs, bool asmVector)
+ {
+ auto const& rowFeSpace = dofmatrix.getRowFeSpace();
+ auto const& colFeSpace = dofmatrix.getColFeSpace();
- dofmatrix->init();
auto matrixOp = matrixOperators[row_col];
auto matrixOpFac = matrixFactors[row_col];
auto vectorOp = vectorOperators[row_col.first];
auto vectorOpFac = vectorFactors[row_col.first];
+ dofmatrix.init(asmMatrix);
+
+ // return if nothing to to.
+ if ((matrixOp.empty() || !asmMatrix) &&
+ (vectorOp.empty() || !asmVector))
+ {
+ dofmatrix.finish();
+ return;
+ }
+
for (auto op : matrixOp)
op->init(rowFeSpace, colFeSpace);
for (auto op : vectorOp)
op->init(rowFeSpace, colFeSpace);
-
+
auto rowLocalView = rowFeSpace.localView();
auto rowIndexSet = rowFeSpace.localIndexSet();
auto colLocalView = colFeSpace.localView();
auto colIndexSet = colFeSpace.localIndexSet();
- for (const auto& element : elements(*meshView)) {
+ for (auto const& element : elements(*meshView)) {
rowLocalView.bind(element);
rowIndexSet.bind(rowLocalView);
colLocalView.bind(element);
colIndexSet.bind(colLocalView);
- if (dofmatrix) {
+ if (asmMatrix) {
ElementMatrix elementMatrix;
bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix,
matrixOp, matrixOpFac);
if (add)
- addElementMatrix(*dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
+ addElementMatrix(dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
}
- if (rhs) {
+ if (asmVector) {
ElementVector elementVector;
bool add = getElementVector(rowLocalView, elementVector,
vectorOp, vectorOpFac);
if (add)
- addElementVector(*rhs, rowIndexSet, elementVector);
+ addElementVector(rhs, rowIndexSet, elementVector);
}
}
- dofmatrix->finish();
+
+ if (asmMatrix) {
+ dofmatrix.finish();
+ matrixAssembled[row_col] = true;
+ }
+ if (asmVector) {
+ vectorAssembled[row_col.first] = true;
+ }
}
@@ -363,15 +407,15 @@ namespace AMDiS
std::list<shared_ptr<OperatorType>>& operators,
std::list<double*> const& factors)
{
- const auto nRows = rowLocalView.tree().finiteElement().size();
- const auto nCols = colLocalView.tree().finiteElement().size();
+ auto const nRows = rowLocalView.tree().finiteElement().size();
+ auto const nCols = colLocalView.tree().finiteElement().size();
AMDIS_TEST_EXIT_DBG(operators.size() == factors.size(),
"Nr. of operators and factors must be consistent!");
// fills the entire matrix with zeroes
- elementMatrix.setSize(nRows, nCols);
- elementMatrix = 0.0;
+ elementMatrix.change_dim(nRows, nCols);
+ set_to_zero(elementMatrix);
bool add = false;
@@ -393,14 +437,14 @@ namespace AMDiS
std::list<shared_ptr<OperatorType>>& operators,
std::list<double*> const& factors)
{
- const auto nRows = rowLocalView.tree().finiteElement().size();
+ auto const nRows = rowLocalView.tree().finiteElement().size();
AMDIS_TEST_EXIT_DBG(operators.size() == factors.size(),
"Nr. of operators and factors must be consistent!");
// Set all entries to zero
- elementVector.resize(nRows);
- elementVector = 0.0;
+ elementVector.change_dim(nRows);
+ set_to_zero(elementVector);
bool add = false;
@@ -422,12 +466,12 @@ namespace AMDiS
ColIndexSet const& colIndexSet,
ElementMatrix const& elementMatrix)
{
- for (size_t i = 0; i < elementMatrix.N(); ++i) {
+ for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
// The global index of the i−th vertex of the element
- const auto row = rowIndexSet.index(i);
- for (size_t j = 0; j < elementMatrix.M(); ++j) {
+ auto const row = rowIndexSet.index(i);
+ for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
// The global index of the j−th vertex of the element
- const auto col = colIndexSet.index(j);
+ auto const col = colIndexSet.index(j);
dofmatrix(row, col) += elementMatrix[i][j];
}
}
@@ -440,9 +484,9 @@ namespace AMDiS
IndexSet const& indexSet,
ElementVector const& elementVector)
{
- for (size_t i = 0; i < elementVector.size(); ++i) {
+ for (size_t i = 0; i < size(elementVector); ++i) {
// The global index of the i-th vertex
- const auto row = indexSet.index(i);
+ auto const row = indexSet.index(i);
dofvector[row] += elementVector[i];
}
}
diff --git a/dune/amdis/ProblemStatBase.hpp b/dune/amdis/ProblemStatBase.hpp
index 29b8ce1a01ac1c1b7ee39a3201c87c676d48af37..82f43269838041c8fa779aaffba26d390a303897 100644
--- a/dune/amdis/ProblemStatBase.hpp
+++ b/dune/amdis/ProblemStatBase.hpp
@@ -6,14 +6,6 @@
#pragma once
#include <string>
-#include <tuple>
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/grid/albertagrid/agrid.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include "Flag.hpp"
@@ -124,29 +116,5 @@ namespace AMDiS
/// Returns the name of the problem.
virtual std::string getName() const = 0;
};
-
-
-
-
- template <int _dim, int _dimworld = _dim, int _degree0 = 1, int... _degrees>
- class ProblemStatTraits
- {
- template <class... Bs>
- using FeSpaceTuple = std::tuple<Bs...>;
-
- template <class Mesh, int deg>
- using Lagrange = Dune::Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
-
- public:
- static constexpr int dim = _dim;
- static constexpr int dimworld = _dimworld;
- static constexpr int nComponents = sizeof...(_degrees) + 1;
-
- // default values
- using Mesh = Dune::AlbertaGrid<dim, dimworld>;
- using FeSpaces = FeSpaceTuple<Lagrange<Mesh, _degree0>,
- Lagrange<Mesh, _degrees>...>;
- };
-
} // end namespace AMDiS
diff --git a/dune/amdis/ProblemStatTraits.hpp b/dune/amdis/ProblemStatTraits.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ca1798fb3c138a5dd4678223acba02fa28fbcaca
--- /dev/null
+++ b/dune/amdis/ProblemStatTraits.hpp
@@ -0,0 +1,59 @@
+#pragma once
+
+#include <tuple>
+
+#include <dune/amdis/Mesh.hpp>
+#include <dune/amdis/utility/GetDegree.hpp>
+
+namespace AMDiS
+{
+ namespace Impl
+ {
+ // by default use Lagrange basis
+ template <class Mesh, int deg>
+ struct DefaultFeSpaceFactory
+ {
+ using type = Dune::Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
+ };
+
+ template <class Mesh>
+ struct DefaultFeSpaceFactory<Mesh, 1>
+ {
+ using type = Dune::Functions::PQ1NodalBasis<typename Mesh::LeafGridView>;
+ };
+
+ template <class Mesh, int deg>
+ using DefaultFeSpaceFactory_t = typename DefaultFeSpaceFactory<Mesh, deg>::type;
+ }
+
+ /// The default traits, that can be used with variable mesh-type.
+ /** The traits class specifies dim, dow and the polynomial degrees of
+ * Lagrange basis functions. Together with a Grid type the necessary parameters
+ * for the \ref ProblemStat class are set.
+ **/
+ template <class MeshType, int _dim, int _dow = _dim, int... _degrees>
+ struct DefaultTraitsMesh
+ {
+ static constexpr int dim = _dim;
+ static constexpr int dimworld = _dow;
+ static constexpr int nComponents = sizeof...(_degrees);
+
+ static_assert( nComponents > 0, "Number of components > 0 required!" );
+
+ using Mesh = MeshType;
+ using FeSpaces = std::tuple<Impl::DefaultFeSpaceFactory_t<Mesh, _degrees>...>;
+ };
+
+
+ /// Specialization of \ref DefaultTraitsMesh for Grid type \ref Dune::YaspGrid.
+ template <int _dim, int _dow = _dim, int... _degrees>
+ using DefaultTraits =
+ DefaultTraitsMesh< Dune::YaspGrid<_dim>, _dim, _dow, _degrees...>;
+
+
+ // DEPRECATED!
+ template <int _dim, int _dimworld = _dim, int _degree0 = 1, int... _degrees>
+ using ProblemStatTraits =
+ DefaultTraits< _dim, _dimworld, _degree0, _degrees...>;
+
+} // end namespace AMDiS
diff --git a/dune/amdis/TermGenerator.hpp b/dune/amdis/TermGenerator.hpp
deleted file mode 100644
index 65a5fa25f216431360ba17657f1f36b7be442a6a..0000000000000000000000000000000000000000
--- a/dune/amdis/TermGenerator.hpp
+++ /dev/null
@@ -1,92 +0,0 @@
-/** \file TermGenerator.hpp */
-
-#pragma once
-
-// std c++ includes
-#include <type_traits>
-
-// AMDiS includes
-#include "OperatorTerm.hpp"
-
-// Dune includes
-#include <dune/functions/common/functionconcepts.hh>
-
-namespace AMDiS
-{
- namespace Impl
- {
- // named categories
- struct _constant {};
- struct _functor {};
- struct _default {};
-
- // TODO: generalize function concept to arbitrary dimensions!
- // TODO: add concept for terms
- // TODO: generalize constants, to accept also vectors
- template <class T>
- using TermCategory =
- std::conditional_t< std::is_arithmetic<T>::value,
- _constant,
- std::conditional_t< Dune::Functions::Concept::isFunction<T, double(Dune::FieldVector<double, 1>)>(),
- _functor,
- std::conditional_t< Dune::Functions::Concept::isFunction<T, double(Dune::FieldVector<double, 2>)>(),
- _functor,
- std::conditional_t< Dune::Functions::Concept::isFunction<T, double(Dune::FieldVector<double, 3>)>(),
- _functor,
- _default > > > >;
-
-
- /// Generator class that takes a type as argument
- /// and returns the corresponding term-type.
- /// By default all arguments are interpreted as constants.
- template <class T, class Category> struct ToTerm;
-
- template <class T>
- struct ToTerm<T, _constant>
- {
- using type = ConstantTerm<T>;
-
- static type get(T const& t)
- {
- return {t};
- }
- };
-
- // if the argument is a functor, that takes a WorldVector
- // and returns a double, see concept \ref CoordsFunctor.
- template <class F>
- struct ToTerm<F, _functor>
- {
- using type = CoordsTerm<F>;
-
- template <class F_>
- static type get(F_&& f)
- {
- return {std::forward<F_>(f)};
- }
- };
-
- // if the argument is already a term return it directly.
- template <class T>
- struct ToTerm<T, _default>
- {
- using type = T const&;
-
- static constexpr type get(T const& t)
- {
- return t;
- }
- };
-
- } // end namespace Impl
-
- template <class T>
- using ToTerm = Impl::ToTerm<std::decay_t<T>, Impl::TermCategory<std::decay_t<T>> >;
-
- template <class T>
- using ToTerm_t = std::decay_t<typename ToTerm<T>::type>;
-
- template <class T>
- inline typename ToTerm<T>::type toTerm(T&& t) { return ToTerm<T>::get(t); }
-
-} // end namespace AMDiS
diff --git a/dune/amdis/Terms.hpp b/dune/amdis/Terms.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d109b264fbf8346e676183d3e96dac31240c8e2f
--- /dev/null
+++ b/dune/amdis/Terms.hpp
@@ -0,0 +1,5 @@
+#pragma once
+
+#include <dune/amdis/terms/ConstantTerm.hpp>
+#include <dune/amdis/terms/CoordsTerm.hpp>
+#include <dune/amdis/terms/DOFVectorTerm.hpp>
diff --git a/dune/amdis/common/ConceptsBase.hpp b/dune/amdis/common/ConceptsBase.hpp
index 315f18c4a972356a96c9955fe38dc274e00b98a8..12f283017c6bbd4b41c388cd9cf07b0e90a13c8f 100644
--- a/dune/amdis/common/ConceptsBase.hpp
+++ b/dune/amdis/common/ConceptsBase.hpp
@@ -27,6 +27,8 @@
#define HAS_MEMBER(name) has_member_ ## name
+#include <dune/functions/common/functionconcepts.hh>
+
namespace AMDiS
{
namespace traits
@@ -65,51 +67,7 @@ namespace AMDiS
template <class T>
using HasResultType = decltype( Impl::HasResultTypeImpl<T>(int{}) );
-
- // Type traits to test whether a class is a functor, i.e. has a operator()
- // -------------------------------------------------------------------------
-
- namespace Impl
- {
- template <class F, class... Args>
- inline auto invoke_functor(F&& f, Args&&... args)
- {
- return std::forward<F>(f)(std::forward<Args>(args)...);
- }
-
- inline no_valid_type invoke_functor(...);
-
- } // end namespace Impl
- template <class, class>
- struct IsFunctor : std::false_type {};
-
- template <class F, class Return, class... Args>
- struct IsFunctor<F, Return(Args...)>
- : std::is_same< decltype( Impl::invoke_functor(std::declval<F>(), std::declval<Args>()...) ),
- Return > {};
-
- template <class, class>
- struct IsFunctorWeak : std::false_type {};
-
- template <class F, class Return, class... Args>
- struct IsFunctorWeak<F, Return(Args...)>
- : std::is_convertible< decltype( Impl::invoke_functor(std::declval<F>(), std::declval<Args>()...) ),
- Return > {};
-
-
- // -------------------------------------------------------------------------
-
-
- template <class>
- struct IsCallable : std::false_type {};
-
- template <class F, class... Args>
- struct IsCallable<F(Args...)>
- : not_< std::is_same< decltype( Impl::invoke_functor(std::declval<F>(), std::declval<Args>()...) ),
- no_valid_type > > {};
-
-
// -------------------------------------------------------------------------
@@ -117,9 +75,7 @@ namespace AMDiS
{
template <class F, class Arg>
inline auto invoke_vector(F&& f, Arg&& arg)
- {
- return std::forward<F>(f)[std::forward<Arg>(arg)];
- }
+ -> decltype( std::forward<F>(f)[std::forward<Arg>(arg)] );
inline no_valid_type invoke_vector(...);
@@ -130,7 +86,7 @@ namespace AMDiS
template <class F, class Return, class Arg>
struct IsVector<F, Return(Arg)>
- : std::is_same< decltype( Impl::invoke_vector(std::declval<F>(), std::declval<Arg>()) ),
+ : std::is_same< decltype( Impl::invoke_vector(std::declval<std::decay_t<F>>(), std::declval<Arg>()) ),
Return > {};
template <class, class>
@@ -138,7 +94,7 @@ namespace AMDiS
template <class F, class Return, class Arg>
struct IsVectorWeak<F, Return(Arg)>
- : std::is_convertible< decltype( Impl::invoke_vector(std::declval<F>(), std::declval<Arg>()) ),
+ : std::is_convertible< decltype( Impl::invoke_vector(std::declval<std::decay_t<F>>(), std::declval<Arg>()) ),
Return > {};
@@ -158,27 +114,61 @@ namespace AMDiS
/// Namespace that implements basic concept checks
- namespace concepts
+ namespace Concepts
{
- /// @defgroup concepts Concepts definitions, using decltype and other c++11 features.
-
- template <class F> // F(Arg)
- constexpr bool Callable = traits::IsCallable<F>::value;
+ /**
+ * \defgroup Concepts Concepts definitions
+ * @{
+ **/
+
+ /// A Collable is a function `F` that can be called with arguments of type `Args...`.
+ /**
+ * To be used as follows: `Concepts::Collable<F, Args...>`. Returns true, if
+ * an instance of `F` can be called by `operator()` with arguments of type
+ * `Args...`.
+ **/
+ template <class F, class... Args>
+ constexpr bool Callable = Dune::Functions::Concept::isCallable<F, Args...>();
+
- template <class F, class U> // F, U=Return(Arg)
- constexpr bool Functor = traits::IsFunctor<F,U>::value;
+ /// A Functor is a function `F` with signature `Signature`.
+ /**
+ * To be used as follows: `Concepts::Functor<F, R(Args...)>`. Returns true, if
+ * an instance of `F` can be called by `operator()` with arguments of type
+ * `Args...` and returns a value of type `R`, i.e. `Signature := R(Args...)`.
+ **/
+ template <class F, class Signature> // F, Signature=Return(Arg)
+ constexpr bool Functor = Dune::Functions::Concept::isFunction<F, Signature>();
+
+ /// A predicate is a function that returns a boolean.
+ template <class F, class... Args>
+ constexpr bool Predicate = Functor<F, bool(Args...)>;
+
+
+ /// Test whether a type `T` supports vector-access `vector[size_t]`.
+ /**
+ * To be used as follows: `Concepts::Vector<V>`. Returns true, if `V` has a
+ * value type, accessible by `Value_t<V>` and an instance of `V` can be called
+ * by square brackets `operator[]` with an index of type `size_t`. The return
+ * value of this element access should be convertible to `Value_t<V>`.
+ **/
template <class T>
constexpr bool Vector = traits::IsVectorWeak<T, Value_t<T>(size_t)>::value;
+
+ /// Test whether a type `T` supports matrix-access `matrix(size_t, size_t)`.
+ /**
+ * To be used as follows: `Concepts::Matrix<M>`. Returns true, if `M` has a
+ * value type, accessible by `Value_t<M>` and an instance of `M` can be called
+ * by round brackets `operator()` with two indices of type `size_t`. The return
+ * value of this element access should be convertible to `Value_t<V>`.
+ **/
template <class T>
- constexpr bool Matrix = traits::IsFunctorWeak<T, Value_t<T>(size_t, size_t)>::value;
-
- } // end namespace concepts
-
+ constexpr bool Matrix = Functor<T, Value_t<T>(size_t, size_t)>;
+
+ /** @} **/
- /// Namespace that implements basic concepts-checks using enable_if and
- /// type-traits from namespace \ref concepts
- namespace requires {}
+ } // end namespace Concepts
} // end namespace AMDiS
diff --git a/dune/amdis/common/ScalarTypes.hpp b/dune/amdis/common/ScalarTypes.hpp
index 039576530169a6c4c51f8b213efed2469b4e5a00..68b42393950cb7cc7224e8462d340163dc05ceda 100644
--- a/dune/amdis/common/ScalarTypes.hpp
+++ b/dune/amdis/common/ScalarTypes.hpp
@@ -15,14 +15,24 @@ namespace AMDiS
} // end namespace traits
- namespace concepts
+ namespace Concepts
{
+ /** \addtogroup Concepts
+ * @{
+ **/
+
+ /// \brief The types following the std type-trait \ref std::is_integral are
+ /// categorized as *integral types*.
template <class T>
constexpr bool Integral = traits::IsIntegral<T>::value;
+ /// \brief The types following the std type-trait \ref std::is_arithmetic are
+ /// categorized as *arithmetic types*.
template <class T>
constexpr bool Arithmetic = traits::IsArithmetic<T>::value;
- } // end namespace concepts
+ /** @} **/
+
+ } // end namespace Concepts
} // end namespace AMDiS
diff --git a/dune/amdis/common/Traits.hpp b/dune/amdis/common/Traits.hpp
index f2e2c8a87582eac808952ffae326831854f38e9b..0cefa62db8db4441d8ee2300551fc49b4f473cd8 100644
--- a/dune/amdis/common/Traits.hpp
+++ b/dune/amdis/common/Traits.hpp
@@ -32,10 +32,10 @@ namespace AMDiS
} // end namespace Impl
template <class T, class U = T>
- using IsAddable = IsCallable<Impl::PlusOp(T, U)>;
+ using IsAddable = bool_<Concepts::Callable<Impl::PlusOp, T, U>>;
template <class T, class U = T>
- using IsMultiplicable = IsCallable<Impl::MultipliesOp(T, U)>;
+ using IsMultiplicable = bool_<Concepts::Callable<Impl::MultipliesOp, T, U>>;
template <class T>
@@ -96,13 +96,17 @@ namespace AMDiS
} // end namespace traits
- namespace concepts
+ namespace Concepts
{
- /// Types support A + B
+ /** \addtogroup Concepts
+ * @{
+ **/
+
+ /// Types support \f$ A + B \f$
template <class A, class B>
constexpr bool Addable = traits::IsAddable<A, B>::value;
- /// Types support A * B
+ /// Types support \f$ A \cdot B \f$
template <class A, class B>
constexpr bool Multiplicable = traits::IsMultiplicable<A, B>::value;
@@ -110,6 +114,8 @@ namespace AMDiS
template <class A, class B>
constexpr bool Compatible = traits::IsCompatible<A, B>::value;
- } // end namespace concepts
+ /** @} **/
+
+ } // end namespace Concepts
} // end namespace AMDiS
diff --git a/dune/amdis/common/TypeDefs.hpp b/dune/amdis/common/TypeDefs.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b5077272314ee1f8f496015a1bd5dd6b067a3519
--- /dev/null
+++ b/dune/amdis/common/TypeDefs.hpp
@@ -0,0 +1,14 @@
+#pragma once
+
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+
+namespace AMDiS
+{
+ namespace Impl
+ {
+ using ElementVector = mtl::dense_vector<double>;
+ using ElementMatrix = mtl::dense2D<double>;
+
+ } // end namespace Impl
+} // end namespace AMDiS
diff --git a/dune/amdis/common/ValueCategory.hpp b/dune/amdis/common/ValueCategory.hpp
index 357d2bb17c6156b4ebfd062df7f9ca6bbdfb7205..0aa8fe04efd205b0cf01043b1d4ba5f8bc8b0d7d 100644
--- a/dune/amdis/common/ValueCategory.hpp
+++ b/dune/amdis/common/ValueCategory.hpp
@@ -14,12 +14,16 @@ namespace AMDiS
struct scalar {};
struct vector {};
struct matrix {};
+ struct unknown {};
} // end namespace tag
/// Category of type T, e.g. scalar, vector matrix, specified by a tag
template <class T, class = void>
- struct ValueCategory;
+ struct ValueCategory
+ {
+ using type = tag::unknown;
+ };
template <class T>
using ValueCategory_t = typename ValueCategory<T>::type;
@@ -53,5 +57,5 @@ namespace AMDiS
{
using type = tag::scalar;
};
-
+
} // end namespace AMDiS
diff --git a/dune/amdis/linear_algebra/istl/SystemMatrix.hpp b/dune/amdis/linear_algebra/istl/SystemMatrix.hpp
index 91efbea345ca5ca000606b11aba2dde6022d7956..118cda2324ea2c8e7385451ef4c6d4f864bdb1e2 100644
--- a/dune/amdis/linear_algebra/istl/SystemMatrix.hpp
+++ b/dune/amdis/linear_algebra/istl/SystemMatrix.hpp
@@ -14,7 +14,7 @@
#include <dune/amdis/linear_algebra/istl/DOFMatrix.hpp>
// for explicit instantiation
-#include <dune/amdis/ProblemStatBase.hpp>
+#include <dune/amdis/ProblemStatTraits.hpp>
namespace AMDiS
{
diff --git a/dune/amdis/linear_algebra/istl/SystemVector.hpp b/dune/amdis/linear_algebra/istl/SystemVector.hpp
index bccc91517cb1e248c46a1c84a90664db577139f0..9cb385b792faafa4994e7756eac0356814bd6c5a 100644
--- a/dune/amdis/linear_algebra/istl/SystemVector.hpp
+++ b/dune/amdis/linear_algebra/istl/SystemVector.hpp
@@ -15,7 +15,7 @@
#include <dune/amdis/linear_algebra/istl/DOFVector.hpp>
// for explicit instantiation
-#include <dune/amdis/ProblemStatBase.hpp>
+#include <dune/amdis/ProblemStatTraits.hpp>
namespace AMDiS
{
diff --git a/dune/amdis/linear_algebra/mtl/BITL_Solver.hpp b/dune/amdis/linear_algebra/mtl/BITL_Solver.hpp
index b17d0d260dc52d5a325272464598e74f2b004995..6be71e292601261af458e13393218b4c49f75f56 100644
--- a/dune/amdis/linear_algebra/mtl/BITL_Solver.hpp
+++ b/dune/amdis/linear_algebra/mtl/BITL_Solver.hpp
@@ -5,9 +5,22 @@
namespace AMDiS
{
-
+ /// Adds default creators for linear solvers based on `BlockMTLMatrix`.
+ /**
+ * Adds creators for block matrix aware solvers.
+ * - *cg*: conjugate gradient method, \see cg_solver_type
+ * - *cgs*: stabilized conjugate gradient mtheod, \see cgs_solver_type
+ * - *bcgs*: stabilized bi-conjugate gradient method, \see bicgstab_type
+ * - *tfqmr*: Transposed-Free Quasi-Minimal Residual method, \see tfqmr_solver_type
+ * - *bcgsl*: stabilized BiCG(ell) method, \see bicgstab_ell_type
+ * - *gmres*: Generalized minimal residula method, \see gmres_type
+ * - *fgmres*: Flexible GMRES, \see fgmres_type
+ * - *minres*: Minimal residul method, \see minres_solver_type
+ * - *gcr*: Generalized conjugate residual method, \see gcr_type
+ * - *umfpack*: external UMFPACK solver, \see UmfpackRunner
+ **/
template <class SubMatrix, size_t _N, size_t _M, class Vector>
- struct DefaultCreators<LinearSolverInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector, Vector>>
+ class DefaultCreators<LinearSolverInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector, Vector>>
{
using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
using SolverBase = LinearSolverInterface<Matrix, Vector, Vector>;
@@ -23,7 +36,7 @@ namespace AMDiS
using Map = CreatorMap<SolverBase>;
- /// Instantiate a new linear solver
+ public:
static void init()
{
auto cg = new SolverCreator<cg_solver_type>;
diff --git a/dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp b/dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp
index 77bd7040beaf4d84dd274ab5a5443245440326ae..10ab36a6c26c3c70aca56b6d8842d2be4cdbf033 100644
--- a/dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp
+++ b/dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp
@@ -44,7 +44,7 @@ namespace AMDiS
inline size_t num_rows(BlockMTLVector<MTLVector, _N> const& vec)
{
size_t nRows = 0;
- For<0, _N>::loop([&](const auto _i) {
+ for_<0, _N>([&](const auto _i) {
nRows += num_rows(std::get<_i>(vec));
});
return nRows;
@@ -68,7 +68,7 @@ namespace AMDiS
template <class MTLVector, size_t _N>
inline void set_to_zero(BlockMTLVector<MTLVector, _N>& vec)
{
- For<0, _N>::loop([&](const auto _i) {
+ for_<0, _N>([&](const auto _i) {
set_to_zero(std::get<_i>(vec));
});
}
@@ -80,7 +80,7 @@ namespace AMDiS
template <class BlockVector, class Vector, class NonConstBlockVector>
class BlockVectorWrapper;
-
+ // specialization for BlockMTLVector
template <class BlockVector, class Vector, class MTLVector, size_t _N>
class BlockVectorWrapper<BlockVector, Vector, BlockMTLVector<MTLVector, _N>>
{
@@ -132,7 +132,6 @@ namespace AMDiS
template <bool Copy>
std::enable_if_t<Copy> assignFrom(bool_<Copy>)
{
- AMDIS_MSG("assignFrom...");
size_t start = 0;
for (size_t r = 0; r < _N; ++r) {
size_t finish = start + num_rows(blockVector[r]);
diff --git a/dune/amdis/linear_algebra/mtl/DOFMatrix.hpp b/dune/amdis/linear_algebra/mtl/DOFMatrix.hpp
index a545d9e605e940f1d492b65d0a5a716278c378b1..99f6134fc830f6eab9d62779f5f69a682237779b 100644
--- a/dune/amdis/linear_algebra/mtl/DOFMatrix.hpp
+++ b/dune/amdis/linear_algebra/mtl/DOFMatrix.hpp
@@ -14,7 +14,7 @@
namespace AMDiS
{
/// \brief The basic container that stores a base matrix and a corresponding
- /// row/column feSpace
+ /// row/column feSpace.
template <class RowFeSpaceType, class ColFeSpaceType,
class ValueType = double>
class DOFMatrix
@@ -26,7 +26,7 @@ namespace AMDiS
/// The type of the finite element space / basis of the column
using ColFeSpace = ColFeSpaceType;
- /// The vector type of the underlying base matrix
+ /// The matrix type of the underlying base matrix
using BaseMatrix = mtl::compressed2D<ValueType>;
/// The index/size - type
@@ -117,15 +117,17 @@ namespace AMDiS
}
/// Create inserter. Assumes that no inserter is currently active on this matrix.
- void init()
+ void init(bool prepareForInsertion)
{
AMDIS_TEST_EXIT(!initialized && !inserter, "Matrix already in insertion mode!");
calculateSlotSize();
matrix->change_dim(rowFeSpace.size(), colFeSpace.size());
- set_to_zero(*matrix);
- inserter = new Inserter(*matrix, slot_size);
- initialized = true;
+ if (prepareForInsertion) {
+ set_to_zero(*matrix);
+ inserter = new Inserter(*matrix, slot_size);
+ initialized = true;
+ }
}
/// Delete inserter -> finish insertion. Must be called in order to fill the
diff --git a/dune/amdis/linear_algebra/mtl/DOFVector.hpp b/dune/amdis/linear_algebra/mtl/DOFVector.hpp
index 2d8b02e19827ba8cbe9c530d60729f30c40a66d5..91ffdf4f35d76205850d01c71029d9da723aa9ec 100644
--- a/dune/amdis/linear_algebra/mtl/DOFVector.hpp
+++ b/dune/amdis/linear_algebra/mtl/DOFVector.hpp
@@ -1,5 +1,9 @@
#pragma once
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#include <string>
#include <memory>
@@ -117,7 +121,7 @@ namespace AMDiS
/// Scale each DOFVector by the factor \p s.
template <class Scalar>
- std::enable_if_t< concepts::Arithmetic<Scalar>, Self&>
+ std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
operator*=(Scalar s)
{
(*vector) *= s;
@@ -126,7 +130,7 @@ namespace AMDiS
/// Sets each DOFVector to the scalar \p s.
template <class Scalar>
- std::enable_if_t< concepts::Arithmetic<Scalar>, Self&>
+ std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
operator=(Scalar s)
{
(*vector) = s;
diff --git a/dune/amdis/linear_algebra/mtl/ITL_Solver.hpp b/dune/amdis/linear_algebra/mtl/ITL_Solver.hpp
index d4f59262a2c9484f66627871090fa50f8a58cd6f..755080b22511a74523bcd446da7ad783d5a90532 100644
--- a/dune/amdis/linear_algebra/mtl/ITL_Solver.hpp
+++ b/dune/amdis/linear_algebra/mtl/ITL_Solver.hpp
@@ -1,9 +1,6 @@
/** \file ITL_Solver.h */
#pragma once
-#ifdef HAVE_CONFIG_H
-# include "config.h"
-#endif
// MTL4 includes
#include <boost/numeric/itl/itl.hpp>
@@ -37,32 +34,6 @@ namespace AMDiS
{
/**
* \ingroup Solver
- *
- * \brief
- * Creator for MTL4 itl-solvers.
- *
- * One of the following solvers can be chosen:
- * - @ref CGSolver "cg" (conjugate gradient method)
- * - @ref CGSSolver "cgs" (squared conjugate gradient method)
- * - @ref BiCGSolver "bicg" (biconjugate gradient method)
- * - @ref BiCGStabSolver "bicgstab" (stabilized BiCG method)
- * - @ref BiCGStab2Solver "bicgstab2" (stabilized BiCG(l) method with l=2)
- * - @ref QMRSolver "qmr" (Quasi-Minimal Residual method)
- * - @ref TFQMRSolver "tfqmr" (Transposed-Free Quasi-Minimal Residual method)
- * - @ref BiCGStabEllSolver "bicgstab_ell" (stabilized BiCG(l) method)
- * - @ref GMResSolver "gmres" (generalized minimal residual method)
- * - @ref IDRsSolver "idr_s" (Induced Dimension Reduction method)
- * - @ref MinResSolver "minres" (minimal residual method)
- * - @ref GcrSolver "gcr" (generalized conjugate residual method)
- * - @ref FGMResSolver "fgmres" (flexible GMRes method)
- * - @ref PreOnly "preonly" (solver that implements pure preconditioning applied to the rhs)
- */
-
- // ===========================================================================
-
- /**
- * \ingroup Solver
- * \class AMDiS::CGSolver
* \brief ITL_Solver <\ref cg_solver_type> implementation of conjugate gradient
* method \implements ITL_Solver
*
@@ -70,7 +41,6 @@ namespace AMDiS
* and can be used for symmetric positive definite system matrices.
* Right preconditioner is ignored.
*/
-
class cg_solver_type
{
public:
@@ -82,18 +52,15 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::CGSSolver
* \brief ITL_Solver <\ref cgs_solver_type> implementation of squared conjugate
* gradient method \implements ITL_Solver
*
* Solves a linear system \f$ Ax=b \f$ by the squared conjugate gradient method
* (CGS). Right preconditioner is ignored.
*/
-
class cgs_solver_type
{
public:
@@ -105,18 +72,15 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::BiCGSolver
* \brief ITL_Solver <\ref bicg_solver_type> implementation of bi-conjugate
* gradient method \implements ITL_Solver
*
* Solves a linear system \f$ Ax=b \f$ by a BiCG method and can be used for
* system matrices.
*/
-
class bicg_solver_type
{
public:
@@ -128,18 +92,15 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::BiCGStabSolver
* \brief ITL_Solver <\ref bicgstab_type> implementation of stabilized
* bi-conjugate gradient method \implements ITL_Solver
*
* Solves a linear system \f$ Ax=b \f$ by a stabilized BiCG method and can be
* used for system matrices.
*/
-
class bicgstab_type
{
public:
@@ -151,18 +112,15 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::BiCGStab2Solver
* \brief ITL_Solver <\ref bicgstab2_type> implementation of BiCGStab(l) method
* with l=2 \implements ITL_Solver
*
* Solves a linear system \f$ Ax=b \f$ by a stabilized BiCG(2) method and can
* be used for system matrices.
*/
-
class bicgstab2_type
{
public:
@@ -174,17 +132,14 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::QMRSolver
* \brief ITL_Solver <\ref qmr_solver_type> implementation of Quasi-Minimal
* Residual method \implements ITL_Solver
*
* Solves a linear system \f$ Ax=b \f$ by the Quasi-Minimal Residual method (QMR).
*/
-
class qmr_solver_type
{
public:
@@ -196,18 +151,15 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::TFQMRSolver
* \brief ITL_Solver <\ref tfqmr_solver_type> implementation of Transposed-Free
* Quasi-Minimal Residual method \implements ITL_Solver
*
* Solves a linear system by the Transposed-Free Quasi-Minimal Residual method
* (TFQMR). Does not use preconditioning currently.
*/
-
class tfqmr_solver_type
{
public:
@@ -219,18 +171,15 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::BiCGStabEllSolver
* \brief ITL_Solver <\ref bicgstab_ell_type> implementation of stabilized
* BiCG(ell) method \implements ITL_Solver
*
* Solves a linear system by a stabilized BiCG(ell) method and can be used for
* system matrices. The parameter ell [3] can be specified.
*/
-
class bicgstab_ell_type
{
int ell;
@@ -246,11 +195,9 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::GMResSolver
* \brief ITL_Solver <\ref gmres_type> implementation of generalized minimal
* residual method \implements ITL_Solver
*
@@ -258,7 +205,6 @@ namespace AMDiS
* The parameter restart [30] is the maximal number of orthogonalized vectors.
* The method is not preconditioned
*/
-
class gmres_type
{
@@ -292,11 +238,8 @@ namespace AMDiS
};
- // ===========================================================================
-
/**
* \ingroup Solver
- * \class AMDiS::IDRsSolver
* \brief ITL_Solver <\ref idr_s_type> implementation of Induced Dimension
* Reduction method \implements ITL_Solver
*
@@ -307,7 +250,6 @@ namespace AMDiS
* algorithms for solving large nonsymmetric linear systems.
* SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035-1062 (2008). (copyright SIAM)
*/
-
class idr_s_type
{
int s;
@@ -323,18 +265,15 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::MinResSolver
* \brief ITL_Solver <\ref minres_solver_type> implementation of minimal
* residual method \implements ITL_Solver
*
* Solves a linear system by the Minres method. Can be used for symmetric
* indefinite systems.
*/
-
class minres_solver_type
{
public:
@@ -346,11 +285,9 @@ namespace AMDiS
}
};
- // ===========================================================================
/**
* \ingroup Solver
- * \class AMDiS::GcrSolver
* \brief ITL_Solver <\ref gcr_type> implementation of generalized conjugate
* residual method \implements ITL_Solver
*
@@ -358,7 +295,6 @@ namespace AMDiS
* method. The parameter restart [30] is the maximal number of orthogonalized
* vectors.
*/
-
class gcr_type
{
int restart;
@@ -374,12 +310,10 @@ namespace AMDiS
return itl::gcr(A, x, b, l, r, iter, restart);
}
};
-
- // ===========================================================================
+
/**
* \ingroup Solver
- * \class AMDiS::FGMResSolver
* \brief ITL_Solver <\ref fgmres_type> implementation of flexible GMRes method
* \implements ITL_Solver
*
@@ -388,7 +322,6 @@ namespace AMDiS
* See reference "A Flexible Inner-Outer Preconditiones GMRES Algorithm",
* Youcef Saad, (1993)
*/
-
class fgmres_type
{
public:
@@ -418,14 +351,11 @@ namespace AMDiS
int restart;
int orthogonalization;
-
};
- // ===========================================================================
-
+
/**
* \ingroup Solver
- * \class AMDiS::PreOnly
* \brief ITL_Solver <\ref preonly_type> implementation of preconditioner as
* \implements ITL_Solver
*
@@ -442,12 +372,31 @@ namespace AMDiS
}
};
+
// ===========================================================================
-
+ /// Adds default creators for linear solvers based on `mtl::compressed2D`.
+ /**
+ * Adds creators for full-matrix aware solvers.
+ * - *cg*: conjugate gradient method, \see cg_solver_type
+ * - *cgs*: stabilized conjugate gradient mtheod, \see cgs_solver_type
+ * - *bicg*: BiCG mtheod, \see bicg_solver_type
+ * - *bcgs*: stabilized bi-conjugate gradient method, \see bicgstab_type
+ * - *bcgs2*: stabilized bi-conjugate gradient method, \see bicgstab2_type
+ * - *qmr*: Quasi-Minimal Residual method, \see qmr_solver_type
+ * - *tfqmr*: Transposed-Free Quasi-Minimal Residual method, \see tfqmr_solver_type
+ * - *bcgsl*: stabilized BiCG(ell) method, \see bicgstab_ell_type
+ * - *gmres*: Generalized minimal residula method, \see gmres_type
+ * - *fgmres*: Flexible GMRES, \see fgmres_type
+ * - *minres*: Minimal residul method, \see minres_solver_type
+ * - *idrs*: Induced Dimension Reduction method, \see idr_s_type
+ * - *gcr*: Generalized conjugate residual method, \see gcr_type
+ * - *preonly*: Just a pply a preconditioner, \see preonly_type
+ * - *umfpack*: external UMFPACK solver, \see UmfpackRunner
+ **/
template <class T, class Param, class Vector>
- struct DefaultCreators<LinearSolverInterface<mtl::compressed2D<T, Param>, Vector, Vector>>
+ class DefaultCreators<LinearSolverInterface<mtl::compressed2D<T, Param>, Vector, Vector>>
{
using Matrix = mtl::compressed2D<T, Param>;
using SolverBase = LinearSolverInterface<Matrix, Vector, Vector>;
@@ -463,7 +412,7 @@ namespace AMDiS
using Map = CreatorMap<SolverBase>;
- /// Instantiate a new linear solver
+ public:
static void init()
{
auto cg = new SolverCreator<cg_solver_type>;
diff --git a/dune/amdis/linear_algebra/mtl/SystemMatrix.hpp b/dune/amdis/linear_algebra/mtl/SystemMatrix.hpp
index f5f88bbd7aa360920030d6d9b415fd78c27446fb..cd70dc53c59616f27a3cb7328be38ba15b5082ea 100644
--- a/dune/amdis/linear_algebra/mtl/SystemMatrix.hpp
+++ b/dune/amdis/linear_algebra/mtl/SystemMatrix.hpp
@@ -15,7 +15,7 @@
#include <dune/amdis/linear_algebra/mtl/DOFMatrix.hpp>
// for explicit instantiation
-#include <dune/amdis/ProblemStatBase.hpp>
+#include <dune/amdis/ProblemStatTraits.hpp>
namespace AMDiS
{
diff --git a/dune/amdis/linear_algebra/mtl/SystemVector.hpp b/dune/amdis/linear_algebra/mtl/SystemVector.hpp
index e5e91a159e026f59725886a7bba0860b7f6a4696..b7f616932f4b0ca6bafef016032f1b876024688f 100644
--- a/dune/amdis/linear_algebra/mtl/SystemVector.hpp
+++ b/dune/amdis/linear_algebra/mtl/SystemVector.hpp
@@ -15,11 +15,10 @@
#include <dune/amdis/linear_algebra/mtl/DOFVector.hpp>
// for explicit instantiation
-#include <dune/amdis/ProblemStatBase.hpp>
+#include <dune/amdis/ProblemStatTraits.hpp>
namespace AMDiS
{
- /// \cond HIDDEN_SYMBOLS
namespace Impl
{
// DOFVectors = std::tuple< DOFVector<FeSpace, ValueType>... >
@@ -57,7 +56,6 @@ namespace AMDiS
}
} // end namespace Impl
- /// \endcond
/// \brief Container that repesents multiple data-Vectors of different value types.
@@ -200,7 +198,7 @@ namespace AMDiS
/// Scale each DOFVector by the factor \p s.
template <class Scalar>
- std::enable_if_t< concepts::Arithmetic<Scalar>, Self&>
+ std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
operator*=(Scalar s)
{
For<0, size()>::loop([this, s](const auto I) {
@@ -211,7 +209,7 @@ namespace AMDiS
/// Sets each DOFVector to the scalar \p s.
template <class Scalar>
- std::enable_if_t< concepts::Arithmetic<Scalar>, Self&>
+ std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
operator=(Scalar s)
{
For<0, size()>::loop([this, s](const auto I) {
diff --git a/dune/amdis/linear_algebra/mtl/UmfpackRunner.hpp b/dune/amdis/linear_algebra/mtl/UmfpackRunner.hpp
index 36ce56f4d79fa13c31d4bb11363ed86c5d94ba4f..83634a2c9221c7b6fc3e37ab570a6b1b6ab7737d 100644
--- a/dune/amdis/linear_algebra/mtl/UmfpackRunner.hpp
+++ b/dune/amdis/linear_algebra/mtl/UmfpackRunner.hpp
@@ -1,5 +1,9 @@
#pragma once
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#ifdef HAVE_UMFPACK
#include <iostream>
@@ -13,16 +17,19 @@
#include <dune/amdis/linear_algebra/mtl/Copy.hpp>
namespace AMDiS
-{
+{
/**
* \ingroup Solver
- * \class AMDiS::UmfPackSolver
- * \brief \implements LinearSolver
- * Wrapper for the external UMFPACK solver:
- * http://www.cise.ufl.edu/research/sparse/umfpack/
+ * \class AMDiS::UmfpackRunner
+ * \brief \implements RunnerInterface
+ * Wrapper for the external
+ * [UMFPACK solver](http://faculty.cse.tamu.edu/davis/suitesparse.html)
*
* This is a direct solver for large sparse matrices.
*/
+ template <class Matrix, class Vector>
+ class UmfpackRunner;
+
template <class Matrix, class Vector>
class UmfpackRunnerBase : public RunnerInterface<Matrix, Vector, Vector>
{
@@ -84,10 +91,6 @@ namespace AMDiS
double alloc_init;
};
-
- template <class Matrix, class Vector>
- class UmfpackRunner;
-
// specialization for block-matrices
template <class SubMatrix, size_t _N, size_t _M, class Vector>
class UmfpackRunner<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
diff --git a/dune/amdis/terms/ConstantTerm.hpp b/dune/amdis/terms/ConstantTerm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d43bfd6ba5234062855a2f6970fb627815b687b3
--- /dev/null
+++ b/dune/amdis/terms/ConstantTerm.hpp
@@ -0,0 +1,52 @@
+#pragma once
+
+#include <type_traits>
+
+#include <dune/amdis/OperatorTermBase.hpp>
+#include <dune/amdis/common/ScalarTypes.hpp>
+
+namespace AMDiS
+{
+ /**
+ * \defgroup OperatorTerms Expressions and terms to be used in GenericOperatorTerm.
+ * @{
+ **/
+
+ /// An expression representing a constant (arithmetic) value, \see constant
+ template <class ValueType>
+ class ConstantTerm
+ {
+ public:
+ // TODO: possibly convert to plain type, in case of expression templates.
+ using value_type = ValueType;
+
+ ConstantTerm(value_type value)
+ : value(value)
+ {}
+
+ template <class Element, class PointList>
+ void init(Element const& element,
+ PointList const& points) { /* do nothing */ }
+
+ value_type const& operator[](size_t iq) const
+ {
+ return value;
+ }
+
+ int getDegree() const
+ {
+ return 0;
+ }
+
+ private:
+ value_type value;
+ };
+
+
+ /// generator function for \ref ConstantTerm expressions
+ template <class T>
+ ConstantTerm<T> constant(T value) { return {value}; }
+
+ /** @} **/
+
+} // end namespace AMDiS
diff --git a/dune/amdis/terms/CoordsTerm.hpp b/dune/amdis/terms/CoordsTerm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..50f30e9b3f117cbc4cb891dcafaaac85f3b29efe
--- /dev/null
+++ b/dune/amdis/terms/CoordsTerm.hpp
@@ -0,0 +1,99 @@
+#pragma once
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <vector>
+#include <type_traits>
+
+#include <dune/istl/bvector.hh>
+
+#include <dune/amdis/OperatorTermBase.hpp>
+#include <dune/amdis/common/Traits.hpp>
+
+namespace AMDiS
+{
+ /**
+ * \addtogroup OperatorTerms
+ * @{
+ **/
+
+ namespace Impl
+ {
+ template <class F, int dow>
+ using CoordsFunctorResult
+ = std::result_of<F(Dune::FieldVector<double, dow>)>;
+
+ template <class F, int dow>
+ constexpr bool CallableDow =
+ Concepts::Callable<F, Dune::FieldVector<double, dow>>;
+
+#ifndef AMDIS_DOW
+ template <class F>
+ using CoordsFunctorResult_t = typename std::conditional_t<
+ CallableDow<F, 1>, CoordsFunctorResult<F, 1>, std::conditional_t<
+ CallableDow<F, 2>, CoordsFunctorResult<F, 2>, std::conditional_t<
+ CallableDow<F, 3>, CoordsFunctorResult<F, 3>, Id<void>> > >::type;
+#else
+ template <class F>
+ using CoordsFunctorResult_t = typename std::conditional_t<
+ CallableDow<F, AMDIS_DOW>, CoordsFunctorResult<F, AMDIS_DOW>, Id<void>>::type;
+#endif
+ }
+
+
+ /// \brief An expression that evaluates to the current coordinate of a dof or
+ /// quadrature point with given index. \see eval.
+ template <class Functor>
+ class CoordsTerm
+ {
+ public:
+ using value_type = Impl::CoordsFunctorResult_t<Functor>;
+ static_assert( !std::is_same<value_type, void>::value,
+ "Template `Functor` can not be used as a functor.");
+
+ template <class F,
+ class = std::enable_if_t<Concepts::Compatible<Functor, F>> >
+ CoordsTerm(F&& f, int degree = 1)
+ : fct(std::forward<F>(f))
+ , degree(degree)
+ {}
+
+ /// \brief Evaluates the functor \ref fct at all quadrature points and
+ /// stores the result in a vector \ref values.
+ template <class Element, class PointList>
+ void init(Element const& element,
+ PointList const& points)
+ {
+ AMDIS_FUNCNAME("CoordsTerm::init()");
+ values.resize(points.size());
+ for (size_t iq = 0; iq < points.size(); ++iq)
+ values[iq] = fct(element.geometry().global(points[iq].position()));
+ }
+
+ value_type const& operator[](size_t iq) const
+ {
+ return values[iq];
+ }
+
+ int getDegree() const
+ {
+ return degree;
+ }
+
+ private:
+ Functor fct;
+ int degree;
+
+ std::vector<value_type> values;
+ };
+
+
+ /// Generator function for \ref CoordsTerm expressions
+ template <class F>
+ CoordsTerm< std::decay_t<F> > eval(F&& f) { return {std::forward<F>(f)}; }
+
+ /** @} **/
+
+} // end namespace AMDiS
diff --git a/dune/amdis/OperatorTerm.hpp b/dune/amdis/terms/DOFVectorTerm.hpp
similarity index 70%
rename from dune/amdis/OperatorTerm.hpp
rename to dune/amdis/terms/DOFVectorTerm.hpp
index ef2bc3e30718ac7314d850ec021852b68b21274d..a0aabcbd980eceb841cd2b7f7dfdc8bc784bd925 100644
--- a/dune/amdis/OperatorTerm.hpp
+++ b/dune/amdis/terms/DOFVectorTerm.hpp
@@ -1,137 +1,29 @@
#pragma once
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
#include <algorithm>
#include <vector>
#include <type_traits>
#include <dune/istl/bvector.hh>
-#include <dune/functions/common/functionconcepts.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/amdis/OperatorTermBase.hpp>
#include <dune/amdis/common/Traits.hpp>
+#include <dune/amdis/utility/GetDegree.hpp>
namespace AMDiS
{
-// some example terms
-// -----------------------------------------------------------------------------
+ /**
+ * \addtogroup OperatorTerms
+ * @{
+ **/
- /// An expression representing a constant (arithmetic) value
- template <class ValueType>
- class ConstantTerm
- {
- public:
- using value_type = ValueType;
-
- ConstantTerm(value_type value)
- : value(value)
- {}
-
- template <class Element, class PointList>
- void init(Element const& element,
- PointList const& points) { /* do nothing */ }
-
- value_type operator[](size_t iq) const
- {
- return value;
- }
-
- int getDegree() const
- {
- return 0;
- }
-
- private:
- value_type value;
- };
-
-
- /// generator function for \ref ConstantTerm expressions
- template <class T>
- std::enable_if_t< concepts::Arithmetic<T>, ConstantTerm<T> >
- constant(T value) { return {value}; }
-
-// -----------------------------------------------------------------------------
-
- namespace Impl
- {
- template <class F, int dow>
- using CoordsFunctorResult
- = std::result_of<F(Dune::FieldVector<double, dow>)>;
-
- template <class F, int dow>
- static constexpr bool isCallableDow()
- { return concepts::Callable<F(Dune::FieldVector<double, dow>)>; }
-
- template <class F>
- using CoordsFunctorResult_t = typename std::conditional_t<
- isCallableDow<F, 1>(), CoordsFunctorResult<F, 1>, std::conditional_t<
- isCallableDow<F, 2>(), CoordsFunctorResult<F, 2>, std::conditional_t<
- isCallableDow<F, 3>(), CoordsFunctorResult<F, 3>, void> > >::type;
- }
-
-
-
- /// An expression that evaluates to the current coordinate of a dof or
- /// quadrature point with given index.
- template <class Functor>
- class CoordsTerm
- {
- public:
- using value_type = Impl::CoordsFunctorResult_t<Functor>;
-
- template <class F,
- class = std::enable_if_t<concepts::Compatible<Functor, F>> >
- CoordsTerm(F&& f, int degree = 1)
- : fct(std::forward<F>(f))
- , degree(degree)
- {}
-
- template <class Element, class PointList>
- void init(Element const& element,
- PointList const& points)
- {
- AMDIS_FUNCNAME("CoordsTerm::init()");
- values.resize(points.size());
- for (size_t iq = 0; iq < points.size(); ++iq)
- values[iq] = fct(element.geometry().global(points[iq].position()));
- }
-
- value_type const& operator[](size_t iq) const
- {
- return values[iq];
- }
-
- int getDegree() const
- {
- return degree;
- }
-
- private:
- Functor fct;
- int degree;
-
- std::vector<value_type> values;
- };
-
-
- /// generator function for \ref CoordsTerm expressions
- template <class F>
- CoordsTerm< std::decay_t<F> > eval(F&& f) { return {std::forward<F>(f)}; }
-
-
-// -----------------------------------------------------------------------------
-
-
- // helper class to extract the polynomial degree of a pqk nodal basis
- template <class FeSpace> struct GetDegree : int_<1> {};
- template <class GV, int k, class ST>
- struct GetDegree<Dune::Functions::PQkNodalBasis<GV, k, ST> > : int_<k> {};
-
-
- /// An expression that evalues a DOFVector at a given element, either on the
- /// dofs or on the quadrature points
+ /// \brief An expression that evalues a DOFVector at a given element, either on the
+ /// dofs or on the quadrature points. \see valueOf.
template <class DOFVectorType>
class DOFVectorTerm
{
@@ -192,13 +84,17 @@ namespace AMDiS
typename Basis::LocalView localView;
typename Basis::LocalIndexSet localIndexSet;
- int degree = GetDegree<Basis>::value;
+ int degree = getPolynomialDegree<Basis>;
std::vector<value_type> values;
};
- /// generator function for \ref DOFVector expressions
+ /// Generator function for \ref DOFVectorTerm expressions.
+ /** Generates a \ref DOFVectorTerm, that takes a \ref DOFVector \p vector,
+ * an optional factor \p factor the value of the DOFVector at the quadrature
+ * points are scaled with.
+ **/
template <class DOFVectorType>
DOFVectorTerm<std::decay_t<DOFVectorType>>
valueOf(DOFVectorType&& vector, double factor = 1.0)
@@ -210,8 +106,9 @@ namespace AMDiS
// -----------------------------------------------------------------------------
- /// An expression that evalues a DOFVector at a given element, either on the
- /// dofs or on the quadrature points, and applies a functor to the value
+ /// \brief An expression that evalues a DOFVector at a given element, either on the
+ /// dofs or on the quadrature points, and applies a functor to the value.
+ /// \see valueOfFunc.
template <class DOFVectorType, class Func>
class DOFVectorFuncTerm
{
@@ -275,14 +172,19 @@ namespace AMDiS
typename Basis::LocalView localView;
typename Basis::LocalIndexSet localIndexSet;
- int degree = GetDegree<Basis>::value;
+ int degree = getPolynomialDegree<Basis>;
int f_deg;
std::vector<value_type> values;
};
- /// generator function for \ref DOFVectorFuncTerm expressions
+ /// Generator function for \ref DOFVectorFuncTerm expressions.
+ /** Generates a \ref DOFVectorFuncTerm, that takes a \ref DOFVector \p vector,
+ * a functor \p f that operators on the value type of the DOFVector, and
+ * optionally a degree argument \p deg that represents an (approximate)
+ * polynomial degree of the functor, for choosing the right quadrature degree.
+ **/
template <class DOFVectorType, class F>
DOFVectorFuncTerm<std::decay_t<DOFVectorType>, std::decay_t<F>>
valueOfFunc(DOFVectorType&& vector, F&& f, int deg = 1)
@@ -294,8 +196,8 @@ namespace AMDiS
// -----------------------------------------------------------------------------
- /// An expression that evalues a DOFVector at a given element, either on the
- /// dofs or on the quadrature points
+ /// \brief An expression that evalues the gradient of a DOFVector at a given
+ /// element, either on the DOFs or on the quadrature points. \see gradientOf.
template <class DOFVectorType>
class GradientTerm
{
@@ -349,7 +251,7 @@ namespace AMDiS
int getDegree() const
{
- return math::max(0, degree-1);
+ return Math::max(0, degree-1);
}
private:
@@ -360,13 +262,17 @@ namespace AMDiS
typename FeSpace::LocalView localView;
typename FeSpace::LocalIndexSet localIndexSet;
- int degree = GetDegree<FeSpace>::value;
+ int degree = getPolynomialDegree<FeSpace>;
std::vector<value_type> values;
};
- /// generator function for an expression of gradients of \ref DOFVectors
+ /// Generator function for an expression of gradients of \ref DOFVectors
+ /** Generates a \ref GradientTerm, that takes a \ref DOFVector \p vector,
+ * and an optional factor \p factor the value of the gradient of the
+ * DOFVector at the quadrature points are scaled with.
+ **/
template <class DOFVectorType>
GradientTerm<std::decay_t<DOFVectorType>>
gradientOf(DOFVectorType&& vector, double factor = 1.0)
@@ -375,8 +281,8 @@ namespace AMDiS
}
- /// An expression that evalues a DOFVector at a given element, either on the
- /// dofs or on the quadrature points
+ /// \brief An expression that evalues a partial derivative of a DOFVector
+ /// at a given element, either on the DOFs or on the quadrature points. \see derivativeOf.
template <class DOFVectorType>
class DerivativeTerm
{
@@ -436,7 +342,7 @@ namespace AMDiS
int getDegree() const
{
- return math::max(0, degree-1);
+ return Math::max(0, degree-1);
}
private:
@@ -450,13 +356,18 @@ namespace AMDiS
typename FeSpace::LocalView localView;
typename FeSpace::LocalIndexSet localIndexSet;
- int degree = GetDegree<FeSpace>::value;
+ int degree = getPolynomialDegree<FeSpace>;
std::vector<value_type> values;
};
- /// generator function for an expression of partial derivatives of \ref DOFVectors
+ /// Generator function for an expression of partial derivatives of \ref DOFVectors
+ /** Generates a \ref DerivativeTerm, that takes a \ref DOFVector \p vector,
+ * a index representing the component of the gradient that should be returned,
+ * and an optional factor \p factor the value of the partial derivative of the
+ * DOFVector at the quadrature points are scaled with.
+ **/
template <class DOFVectorType>
DerivativeTerm<std::decay_t<DOFVectorType>>
derivativeOf(DOFVectorType&& vector, int component, double factor = 1.0)
@@ -464,4 +375,6 @@ namespace AMDiS
return {std::forward<DOFVectorType>(vector), component, factor};
}
+ /** @} **/
+
} // end namespace AMDiS
diff --git a/dune/amdis/terms/TermGenerator.hpp b/dune/amdis/terms/TermGenerator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..56169d29c9d9a1e9ea9fca74e385ecc062f41935
--- /dev/null
+++ b/dune/amdis/terms/TermGenerator.hpp
@@ -0,0 +1,120 @@
+/** \file TermGenerator.hpp */
+
+#pragma once
+
+// std c++ includes
+#include <type_traits>
+
+// AMDiS includes
+#include <dune/amdis/terms/ConstantTerm.hpp>
+#include <dune/amdis/terms/CoordsTerm.hpp>
+
+#include <dune/amdis/common/ConceptsBase.hpp>
+#include <dune/amdis/common/ValueCategory.hpp>
+
+namespace AMDiS
+{
+ namespace Impl
+ {
+ // annonymous namespace to hide tags
+ namespace
+ {
+ struct _constant {};
+ struct _functor {};
+ struct _default {};
+ }
+
+
+ template <class T>
+ constexpr bool ConstantCategory =
+ !std::is_same<ValueCategory_t<T>, tag::unknown>::value;
+
+
+ // TODO: generalize to functors with arbitrary return values
+#ifndef AMDIS_DOW
+ template <class T>
+ constexpr bool FunctorCategory =
+ Concepts::Functor<T, double(Dune::FieldVector<double, 1>)> ||
+ Concepts::Functor<T, double(Dune::FieldVector<double, 2>)> ||
+ Concepts::Functor<T, double(Dune::FieldVector<double, 3>)>;
+#else
+ template <class T>
+ constexpr bool FunctorCategory =
+ Concepts::Functor<T, double(Dune::FieldVector<double, AMDIS_DOW>)>;
+#endif
+
+
+ // TODO: add concept for terms
+ template <class T>
+ using TermCategory =
+ std::conditional_t< ConstantCategory<T>, _constant,
+ std::conditional_t< FunctorCategory<T>, _functor,
+ /* OperatorTerms */ _default > >;
+
+
+ // forward declaration
+ template <class T, class Tag> struct ToTerm;
+
+ template <class T>
+ struct ToTerm<T, _constant>
+ {
+ using type = ConstantTerm<T>;
+
+ static type get(T const& t)
+ {
+ return {t};
+ }
+ };
+
+ // if the argument is a functor, that takes a WorldVector
+ // and returns a double, see concept \ref CoordsFunctor.
+ template <class F>
+ struct ToTerm<F, _functor>
+ {
+ using type = CoordsTerm<F>;
+
+ template <class F_>
+ static type get(F_&& f)
+ {
+ return {std::forward<F_>(f)};
+ }
+ };
+
+ // if the argument is already a term return it directly.
+ template <class T>
+ struct ToTerm<T, _default>
+ {
+ using type = T const&;
+
+ static constexpr type get(T const& t)
+ {
+ return t;
+ }
+ };
+
+ } // end namespace Impl
+
+
+ /// Generator class that takes a type as argument and returns the corresponding term-type.
+ /** The arguments are transformed into terms in the following way:
+ * - If type `T` is arithmetic or any small vector/matrix type, generate a
+ * constant term, see \ref ConstantTerm
+ * - If the type `T` models a Functor concepts, that maps a WorldVector to a
+ * double type, construct a \ref CoordsTerm
+ * - In all other cases, the argument is retured without modification.
+ **/
+ template <class T>
+ using ToTerm = Impl::ToTerm<std::decay_t<T>, Impl::TermCategory<std::decay_t<T>> >;
+
+ template <class T>
+ using ToTerm_t = std::decay_t<typename ToTerm<T>::type>;
+
+
+ /// Generator function for term transformation, \see ToTerm
+ template <class T>
+ inline typename ToTerm<T>::type toTerm(T&& t)
+ {
+ return ToTerm<T>::get(std::forward<T>(t));
+ }
+
+} // end namespace AMDiS
diff --git a/dune/amdis/utility/GetDegree.hpp b/dune/amdis/utility/GetDegree.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0f2b419b1ab3d0e33d133dbe8034842831fdf5bb
--- /dev/null
+++ b/dune/amdis/utility/GetDegree.hpp
@@ -0,0 +1,40 @@
+#pragma once
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/functions/functionspacebases/lagrangedgbasis.hh>
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/pq1nodalbasis.hh>
+#include <dune/functions/functionspacebases/taylorhoodbasis.hh>
+
+#include <dune/amdis/common/Mpl.hpp>
+
+namespace AMDiS
+{
+ namespace Impl
+ {
+ // helper class to extract the polynomial degree of a pqk nodal basis
+ template <class FeSpace>
+ struct GetDegree : int_<1> {};
+
+ template <class GV, int k, class ST>
+ struct GetDegree<Dune::Functions::PQkNodalBasis<GV, k, ST> > : int_<k> {};
+
+ template <class GV, int k, class ST>
+ struct GetDegree<Dune::Functions::LagrangeDGBasis<GV, k, ST> > : int_<k> {};
+
+ template <class GV, class ST>
+ struct GetDegree<Dune::Functions::PQ1NodalBasis<GV, ST> > : int_<1> {};
+
+ template <class GV, class ST>
+ struct GetDegree<Dune::Functions::TaylorHoodBasis<GV, ST> > : int_<2> {};
+
+ } // end namespace Impl
+
+ /// \brief Return the polynomial degree of a given FESpace. Falls back to 1 if unknown.
+ template <class FeSpace>
+ constexpr int getPolynomialDegree = Impl::GetDegree<FeSpace>::value;
+
+} // end namespace AMDiS
diff --git a/init/pfc.json.2d b/init/pfc.json.2d
index 7591b4ea0fe1018b27ca7a6f6d3f952d0c363f77..2d38d18e96d8accb86e3acd2c6ad45fcfbe5a8c3 100644
--- a/init/pfc.json.2d
+++ b/init/pfc.json.2d
@@ -3,12 +3,12 @@
"pfcMesh": {
"macro file name": "./macro/macro.stand.2d",
- "global refinements": 4,
+ "global refinements": 5,
"min corner": "0,0",
"max corner": "1,1",
"num cells": "2,2",
- "dimension": "10,10"
+ "dimension": "50,50"
},
"pfc": {
@@ -20,6 +20,7 @@
"relative tolerance": 1e-6,
"info": 10,
"ell": 3,
+ "print cycle": 10,
"right precon": "pfc"
},
diff --git a/src/BlockPreconditioner.hpp b/src/BlockPreconditioner.hpp
index f0539df185be49baf837e127c04c46e7e22daaa5..19dac416b354d49a41b9bc34c5af453b4b746b0e 100644
--- a/src/BlockPreconditioner.hpp
+++ b/src/BlockPreconditioner.hpp
@@ -23,6 +23,7 @@ namespace AMDiS
A = &A_;
A_.getRanges(rows, cols);
+ initialized = true;
}
virtual void solve(Vector const& b, Vector& x) const override
@@ -35,18 +36,25 @@ namespace AMDiS
AMDIS_ERROR_EXIT("Must be implemented in derived classes!\n");
}
- mtl::irange const& getRowRange(size_t i) const
+ template <size_t I>
+ mtl::irange const& getRowRange(const index_<I> i = {}) const
{
- return rows[i];
+ AMDIS_TEST_EXIT_DBG(initialized,
+ "BlockPreconditioner not initialized. Call init() before.");
+ return std::get<I>(rows);
}
- mtl::irange const& getColRange(size_t i) const
+ template <size_t I>
+ mtl::irange const& getColRange(const index_<I> i = {}) const
{
- return cols[i];
+ AMDIS_TEST_EXIT_DBG(initialized,
+ "BlockPreconditioner not initialized. Call init() before.");
+ return std::get<I>(cols);
}
protected:
Matrix const* A = NULL;
+ bool initialized = false;
std::array<mtl::irange, Matrix::N()> rows;
std::array<mtl::irange, Matrix::M()> cols;
diff --git a/src/MTLPfcPrecon.hpp b/src/MTLPfcPrecon.hpp
index 16d9db3063e99032aacac387d6770b5081b1e0cb..d86c667249f5c7305c17734ab01d9c6f3254a688 100644
--- a/src/MTLPfcPrecon.hpp
+++ b/src/MTLPfcPrecon.hpp
@@ -42,13 +42,13 @@ namespace AMDiS
CreatorInterfaceName<LinearSolverType>* creator;
- creator = named(CreatorMap<LinearSolverType>::getCreator(solverNameM, "precon_pfc->M->solver"));
+ creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameM, "precon_pfc->M->solver") );
solverM = creator->create("precon_pfc->M");
- creator = named(CreatorMap<LinearSolverType>::getCreator(solverNameMpL, "precon_pfc->MpL->solver"));
+ creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL, "precon_pfc->MpL->solver") );
solverMpL = creator->create("precon_pfc->MpL");
- creator = named(CreatorMap<LinearSolverType>::getCreator(solverNameMpL2, "precon_pfc->MpL2->solver"));
+ creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL2, "precon_pfc->MpL2->solver") );
solverMpL2 = creator->create("precon_pfc->MpL2");
}
@@ -83,19 +83,19 @@ namespace AMDiS
void solveMpL(VectorX& x, const VectorB& b) const
{
SolverInfo solverInfo("precon_pfc->MpL");
- solverM->solve(MpL, x, b, solverInfo);
+ solverMpL->solve(MpL, x, b, solverInfo);
}
template <typename VectorX, typename VectorB>
void solveMpL2(VectorX& x, const VectorB& b) const
{
SolverInfo solverInfo("precon_pfc->MpL2");
- solverM->solve(MpL2, x, b, solverInfo);
+ solverMpL2->solve(MpL2, x, b, solverInfo);
}
- MTLMatrix const& getM() const { return matM ? *matM : (*Super::A)(2_c,2_c); } // < u, v >
- MTLMatrix const& getL0() const { return matL0 ? *matL0 : (*Super::A)(1_c,0_c); } // < M0*tau*grad(u), grad(v) >
- MTLMatrix const& getL() const { return matL ? *matL : (*Super::A)(2_c,1_c); } // < grad(u), grad(v) >
+ MTLMatrix const& getM() const { return matM ? *matM : (*Super::A)(2_c, 2_c); } // < u, v >
+ MTLMatrix const& getL0() const { return matL0 ? *matL0 : (*Super::A)(1_c, 0_c); } // < M0*tau*grad(u), grad(v) >
+ MTLMatrix const& getL() const { return matL ? *matL : (*Super::A)(2_c, 1_c); } // < grad(u), grad(v) >
double getTau() const { return *tau; }
diff --git a/src/MTLPfcPrecon.inc.hpp b/src/MTLPfcPrecon.inc.hpp
index 775553f94d66c8d1dd128b9c503ad5e48962b725..68cbdc7e9fe81848c3d8fefa9b097b0fb29ee4fc 100644
--- a/src/MTLPfcPrecon.inc.hpp
+++ b/src/MTLPfcPrecon.inc.hpp
@@ -4,17 +4,18 @@ namespace AMDiS
template <class Matrix, class Vector>
void MTLPfcPrecon<Matrix, Vector>::init(Matrix const& A)
{
- AMDIS_TEST_EXIT(tau != NULL, "Preconditioner not initialized!");
+ AMDIS_TEST_EXIT(tau, "Preconditioner not initialized!");
Super::init(A);
using size_type = typename MTLMatrix::size_type;
size_type n = num_rows(getM());
double delta = std::sqrt(M0 * getTau());
+ AMDIS_MSG("delta = " << delta);
// helper-matrix MpL = M + delta*L
MpL.change_dim(n, n);
- MpL = getM() + (1.0/delta)*getL0();
+ MpL = getM() + (1.0/delta) * getL0();
// helper-matrix MpL = M + sqrt(delta)*L
MpL2.change_dim(n, n);
@@ -32,13 +33,13 @@ namespace AMDiS
{
x.change_dim(num_rows(b));
- Vector x0(x[Super::getColRange(0)]);
- Vector x1(x[Super::getColRange(1)]);
- Vector x2(x[Super::getColRange(2)]);
+ Vector x0(x[Super::getColRange(0_c)]);
+ Vector x1(x[Super::getColRange(1_c)]);
+ Vector x2(x[Super::getColRange(2_c)]);
- const Vector b0(b[Super::getRowRange(0)]);
- const Vector b1(b[Super::getRowRange(1)]);
- const Vector b2(b[Super::getRowRange(2)]);
+ const Vector b0(b[Super::getRowRange(0_c)]);
+ const Vector b1(b[Super::getRowRange(1_c)]);
+ const Vector b2(b[Super::getRowRange(2_c)]);
double delta = std::sqrt(M0 * getTau());
@@ -54,7 +55,7 @@ namespace AMDiS
tmp = getM() * y1; // tmp := M*y1
solveMpL2(x1, tmp); // (M+eps*sqrt(delta)K) * x1 = tmp
- x0-= (1.0/delta)*x1; // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
+ x0-= (1.0/delta)*x1; // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
tmp = b2 - getL() * x1; // tmp := b2 - K*x1
solveM(x2, tmp);
diff --git a/src/amdis.cc b/src/amdis.cc
index bb002d7ee188cf731fa04d2779d2ee9df925f45d..a231300dd68b173fff60737b68b06c61aa9e2d96 100644
--- a/src/amdis.cc
+++ b/src/amdis.cc
@@ -5,6 +5,7 @@
#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/ProblemStat.hpp>
+#include <dune/amdis/Terms.hpp>
using namespace AMDiS;
diff --git a/src/ellipt.cc b/src/ellipt.cc
index 5019e89afb33b186add119114c5ef6c357407860..4260dc27f870efff043ed8b5e5e34bf3e1af1c1a 100644
--- a/src/ellipt.cc
+++ b/src/ellipt.cc
@@ -7,6 +7,7 @@
#include <dune/amdis/AdaptInstationary.hpp>
#include <dune/amdis/ProblemInstat.hpp>
#include <dune/amdis/ProblemStat.hpp>
+#include <dune/amdis/Terms.hpp>
#include <dune/amdis/common/Literals.hpp>
using namespace AMDiS;
diff --git a/src/heat.cc b/src/heat.cc
index e7af6ef07901ff9ffe58c7fdda043fdf19d26861..5f251aca474b791fd2c5fbbb5a29a0b357391e98 100644
--- a/src/heat.cc
+++ b/src/heat.cc
@@ -7,6 +7,7 @@
#include <dune/amdis/AdaptInstationary.hpp>
#include <dune/amdis/ProblemInstat.hpp>
#include <dune/amdis/ProblemStat.hpp>
+#include <dune/amdis/Terms.hpp>
#include <dune/amdis/common/Literals.hpp>
using namespace AMDiS;
diff --git a/src/navier_stokes.cc b/src/navier_stokes.cc
index aa4d01ac2c74ff117012b3d7e5aa34c20accbe39..32097eef327aa4192c77acb389e40289d349672f 100644
--- a/src/navier_stokes.cc
+++ b/src/navier_stokes.cc
@@ -6,6 +6,7 @@
#include <dune/amdis/AdaptInstationary.hpp>
#include <dune/amdis/ProblemStat.hpp>
#include <dune/amdis/ProblemInstat.hpp>
+#include <dune/amdis/Terms.hpp>
#ifdef DOW
#undef DOW
@@ -37,7 +38,7 @@ int main(int argc, char** argv)
Parameters::get("stokes->boundary velocity", vel);
using DOFVectorU = std::decay_t< decltype(prob.getSolution<0>()) >;
- using DOFVectorP = std::decay_t< decltype(prob.getSolution<DOW>()) >;
+// using DOFVectorP = std::decay_t< decltype(prob.getSolution<DOW>()) >;
std::array<DOFVectorU*, DOW> velocity;
For<0,DOW>::loop([&](auto const _i) { velocity[_i] = &prob.getSolution(_i); });
diff --git a/src/pfc.cc b/src/pfc.cc
index 596b2939eabe78308e2a66e066c2e663663a1bf7..f0b2414d1bde787378c9fcec9dc30f728afcd8b9 100644
--- a/src/pfc.cc
+++ b/src/pfc.cc
@@ -9,32 +9,17 @@
#include <dune/amdis/AdaptInstationary.hpp>
#include <dune/amdis/ProblemStat.hpp>
#include <dune/amdis/ProblemInstat.hpp>
+#include <dune/amdis/Terms.hpp>
#include "MTLPfcPrecon.hpp"
using namespace AMDiS;
+using PfcParam =
+ DefaultTraitsMesh<Dune::YaspGrid<AMDIS_DIM>, AMDIS_DIM, AMDIS_DOW, 2, 2, 2>;
-class PfcParam
-{
- template <class Mesh, int deg>
- using Lagrange = Dune::Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
-
- template <class Mesh, int... degs>
- using FeSpaceTuple = std::tuple<Lagrange<Mesh, degs>...>;
-
-public:
- static constexpr int dim = AMDIS_DIM;
- static constexpr int dimworld = AMDIS_DOW;
- static constexpr int nComponents = 3;
-
- // default values
- using Mesh = Dune::YaspGrid<dim>;
- using FeSpaces = FeSpaceTuple<Mesh, 2, 2, 2>;
-};
// 3 component with polynomial degree 1
-// using PfcParam = ProblemStatTraits<AMDIS_DIM, AMDIS_DOW, 1, 1, 1>;
using PfcProblem = ProblemStat<PfcParam>;
using PfcProblemInstat = ProblemInstat<PfcParam>;
@@ -67,19 +52,31 @@ int main(int argc, char** argv)
auto& Nu = prob.getSolution<2>();
using Op = PfcProblem::OperatorType;
- auto opLhs01 = std::make_shared<Op>(); // < [-(1+r) + 3*psi^2]*u, v > + < 2*grad(u), grad(v) >
- opLhs01->addZOT( -(1.0 + r) );
- opLhs01->addZOT( valueOfFunc(Psi, [](auto psi) { return 3.0 * std::pow(psi,2); }, 2) );
- opLhs01->addSOT( 2.0 );
+
+ // < [-(1+r) - 3*psi^2]*u, v > + < 2*grad(u), grad(v) >
+ auto opLhs01_ = Op::create<0>( -(1.0 + r) );
+ opLhs01_->addZOT( valueOfFunc(Psi, [](auto psi) { return -3.0 * Math::pow<2>(psi); }, 2) );
+ opLhs01_->addSOT( 2.0 );
+ auto opLhs01 = std::make_pair(opLhs01_, true);
- auto opRhs0 = Op::create<0>(valueOfFunc(Psi, [](auto psi) { return -2.0 * std::pow(psi,3); }, 3)); // < -2*psi^3, v >
+ // < -2*psi^3, v >
+ auto opRhs0 = Op::create<0>(valueOfFunc(Psi, [](auto psi) { return -2.0 * Math::pow<3>(psi); }, 3));
+
+ // < psi, v >
auto opRhs1 = Op::create<0>(valueOf(Psi));
- auto opM = Op::create<0>(1.0); // < u, v >
- auto opL = Op::create<2>(1.0); // < grad(u), grad(v) >
+ // < u, v >
+ auto opM = std::make_pair(Op::create<0>(1.0), false);
+ // < grad(u), grad(v) >
+ auto opL = std::make_pair(Op::create<2>(1.0), false);
+
+ // < M0*tau * grad(u), grad(v) >
double* tauPtr = adaptInfo.getTimestepPtr();
- auto opL0= Op::create<2>([M0, tauPtr](auto) { return M0 * (*tauPtr); }); // < M0*tau * grad(u), grad(v) >
+ auto opL0= std::make_pair(Op::create<2>([M0, tauPtr](auto) { return M0 * (*tauPtr); }), false);
+
+
+ // === Define the linear system ===
prob.addMatrixOperator({{ {{0,0}, opM}, {{0,1}, opLhs01}, {{0,2}, opL},
{{1,0}, opL0}, {{1,1}, opM},
@@ -88,12 +85,13 @@ int main(int argc, char** argv)
prob.addVectorOperator({{ {0, opRhs0},
{1, opRhs1} }});
- // initial solution
+ // === set the initial solution ===
Mu = 0.0;
Nu = 0.0;
std::srand( 12345 );
Psi.interpol([psi_mean](auto const& x) { return 0.5*(std::rand()/double(RAND_MAX) - 0.5) + psi_mean; });
+ // === configure preconditioner, if necessary ===
if (creator->precon) {
creator->precon->setData(adaptInfo.getTimestepPtr(), M0);
}
diff --git a/src/stokes.cc b/src/stokes.cc
index ed842bd7d5926d66880164901555f6bbc369d479..f0359c192f77764d56b64311511bc51568ca8f94 100644
--- a/src/stokes.cc
+++ b/src/stokes.cc
@@ -4,6 +4,7 @@
#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/ProblemStat.hpp>
+#include <dune/amdis/Terms.hpp>
#ifdef DOW
#undef DOW