Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • feature/blockgssteps_autoCopy
  • feature/cmakelists-sources-target
  • feature/incomplete-cholesky-rebased
  • feature/istl-preconditioners
  • feature/optional-ignore
  • feature/update-buildsystem
  • feature/update-to-clang-7
  • feature/use-smart-ptr-ignorenodes
  • feature/whitespace-fix
  • fix/error-norm
  • fix_linking_module
  • flexible-loopsolver-max
  • generalized-blockgsstep-rebased
  • implement-overlappingblockgsstep
  • make-getiterationstep-return-shared-ptr
  • master
  • more-features-for-cholmodsolver
  • new_interface
  • releases/2.0-1
  • releases/2.1-1
  • releases/2.10
  • releases/2.2-1
  • releases/2.3-1
  • releases/2.4-1
  • releases/2.5-1
  • releases/2.6-1
  • releases/2.7
  • releases/2.8
  • releases/2.9
  • subversion->git
30 results

Target

Select target project
  • lisa_julia.nebel_at_tu-dresden.de/dune-solvers
  • patrick.jaap_at_tu-dresden.de/dune-solvers
  • burchardt_at_igpm.rwth-aachen.de/dune-solvers
  • agnumpde/dune-solvers
4 results
Select Git revision
  • feature/cmakelists-sources-target
  • feature/incomplete-cholesky-rebased
  • feature/istl-preconditioners
  • feature/optional-ignore
  • fix/dune_deprecated_macro
  • fix_linking_module
  • flexible-loopsolver-max
  • generalized-blockgsstep-rebased
  • make-getiterationstep-return-shared-ptr
  • master
  • new_interface
  • releases/2.0-1
  • releases/2.1-1
  • releases/2.2-1
  • releases/2.3-1
  • releases/2.4-1
  • releases/2.5-1
  • releases/2.6-1
  • use-keyword-signature-of-target_link_libraries
  • subversion->git
20 results
Show changes
Showing
with 800 additions and 1065 deletions
...@@ -3,8 +3,13 @@ add_subdirectory("iterationsteps") ...@@ -3,8 +3,13 @@ add_subdirectory("iterationsteps")
add_subdirectory("norms") add_subdirectory("norms")
add_subdirectory("operators") add_subdirectory("operators")
add_subdirectory("solvers") add_subdirectory("solvers")
add_subdirectory("test")
add_subdirectory("transferoperators") add_subdirectory("transferoperators")
target_sources(dunesolvers PRIVATE
iterationsteps/blockgssteps.cc
solvers/criterion.cc)
install(FILES install(FILES
computeenergy.hh computeenergy.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/solvers) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/solvers)
SUBDIRS = common iterationsteps norms operators solvers test transferoperators
dune_solversdir = $(includedir)/dune/solvers
dune_solvers_HEADERS = computeenergy.hh
EXTRA_DIST = CMakeLists.txt
include $(top_srcdir)/am/global-rules
install(FILES install(FILES
arithmetic.hh
boxconstraint.hh boxconstraint.hh
canignore.hh canignore.hh
copyorreference.hh
defaultbitvector.hh
genericvectortools.hh genericvectortools.hh
interval.hh interval.hh
numproc.hh numproc.hh
permutationmanager.hh permutationmanager.hh
preconditioner.hh preconditioner.hh
resize.hh
staticmatrixtools.hh staticmatrixtools.hh
tuplevector.hh
wrapownshare.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/solvers/common) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/solvers/common)
SUBDIRS =
commondir = $(includedir)/dune/solvers/common
common_HEADERS = arithmetic.hh \
boxconstraint.hh \
genericvectortools.hh \
canignore.hh \
interval.hh \
numproc.hh \
permutationmanager.hh \
preconditioner.hh \
staticmatrixtools.hh
EXTRA_DIST = CMakeLists.txt
include $(top_srcdir)/am/global-rules
#ifndef ARITHMETIC_HH
#define ARITHMETIC_HH
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/diagonalmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/scaledidmatrix.hh>
#include <dune/common/typetraits.hh>
/** \brief Namespace containing helper classes and functions for arithmetic operations
*
* Everything in this namespace is experimental and might change
* in the near future. Especially the naming of namespace, structs,
* and functions is not final.
*/
namespace Arithmetic
{
/** \brief Class to identify scalar types
*
* Specialize this class for all types that can be used
* like scalar quantities.
*/
template<class T>
struct ScalarTraits
{
enum { isScalar=false };
};
template<>
struct ScalarTraits<float>
{
enum { isScalar=true };
};
template<>
struct ScalarTraits<double>
{
enum { isScalar=true };
};
template<class T>
struct ScalarTraits<Dune::FieldVector<T,1> >
{
enum { isScalar=true };
};
template<class T>
struct ScalarTraits<Dune::FieldMatrix<T,1,1> >
{
enum { isScalar=true };
};
template<class T>
struct ScalarTraits<Dune::DiagonalMatrix<T,1> >
{
enum { isScalar=true };
};
template<class T>
struct ScalarTraits<Dune::ScaledIdentityMatrix<T,1> >
{
enum { isScalar=true };
};
/** \brief Class to identify matrix types and extract information
*
* Specialize this class for all types that can be used like a matrix.
*/
template<class T>
struct MatrixTraits
{
enum { isMatrix=false };
enum { rows=-1};
enum { cols=-1};
};
template<class T, int n, int m>
struct MatrixTraits<Dune::FieldMatrix<T,n,m> >
{
enum { isMatrix=true };
enum { rows=n};
enum { cols=m};
};
template<class T, int n>
struct MatrixTraits<Dune::DiagonalMatrix<T,n> >
{
enum { isMatrix=true };
enum { rows=n};
enum { cols=n};
};
template<class T, int n>
struct MatrixTraits<Dune::ScaledIdentityMatrix<T,n> >
{
enum { isMatrix=true };
enum { rows=n};
enum { cols=n};
};
template<class T>
struct MatrixTraits<Dune::BCRSMatrix<T> >
{
enum { isMatrix=true };
};
/** \brief Internal helper class for product operations
*
*/
template<class A, class B, class C, bool AisScalar, bool BisScalar, bool CisScalar>
struct ProductHelper
{
static void addProduct(A& a, const B& b, const C& c)
{
b.umv(c, a);
}
static void subtractProduct(A& a, const B& b, const C& c)
{
b.mmv(c, a);
}
};
// Internal helper class for scaled product operations (i.e., b is always a scalar)
template<class A, class B, class C, class D, bool AisScalar, bool CisScalar, bool DisScalar>
struct ScaledProductHelper
{
static void addProduct(A& a, const B& b, const C& c, const D& d)
{
c.usmv(b, d, a);
}
static void subtractProduct(A& a, const B& b, const C& c, const D& d)
{
c.usmv(-b, d, a);
}
};
template<class T, int n, int m, int k>
struct ProductHelper<Dune::FieldMatrix<T,n,k>, Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,m,k>, false, false, false>
{
typedef Dune::FieldMatrix<T,n,k> A;
typedef Dune::FieldMatrix<T,n,m> B;
typedef Dune::FieldMatrix<T,m,k> C;
static void addProduct(A& a, const B& b, const C& c)
{
for (size_t row = 0; row < n; ++row)
{
for (size_t col = 0 ; col < k; ++col)
{
for (size_t i = 0; i < m; ++i)
a[row][col] += b[row][i]*c[i][col];
}
}
}
static void subtractProduct(A& a, const B& b, const C& c)
{
for (size_t row = 0; row < n; ++row)
{
for (size_t col = 0 ; col < k; ++col)
{
for (size_t i = 0; i < m; ++i)
a[row][col] -= b[row][i]*c[i][col];
}
}
}
};
template<class T, int n, int m, int k>
struct ScaledProductHelper<Dune::FieldMatrix<T,n,k>, T, Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,m,k>, false, false, false>
{
typedef Dune::FieldMatrix<T,n,k> A;
typedef Dune::FieldMatrix<T,n,m> C;
typedef Dune::FieldMatrix<T,m,k> D;
static void addProduct(A& a, const T& b, const C& c, const D& d)
{
for (size_t row = 0; row < n; ++row) {
for (size_t col = 0 ; col < k; ++col) {
for (size_t i = 0; i < m; ++i)
a[row][col] += b * c[row][i]*d[i][col];
}
}
}
static void subtractProduct(A& a, const T& b, const C& c, const D& d)
{
for (size_t row = 0; row < n; ++row) {
for (size_t col = 0 ; col < k; ++col) {
for (size_t i = 0; i < m; ++i)
a[row][col] -= b * c[row][i]*d[i][col];
}
}
}
};
template<class T, int n>
struct ProductHelper<Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
{
typedef Dune::DiagonalMatrix<T,n> A;
typedef A B;
typedef B C;
static void addProduct(A& a, const B& b, const C& c)
{
for (size_t i=0; i<n; ++i)
a.diagonal(i) += b.diagonal(i)*c.diagonal(i);
}
static void subtractProduct(A& a, const B& b, const C& c)
{
for (size_t i=0; i<n; ++i)
a.diagonal(i) -= b.diagonal(i)*c.diagonal(i);
}
};
template<class T, int n>
struct ScaledProductHelper<Dune::DiagonalMatrix<T,n>, T, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
{
typedef Dune::DiagonalMatrix<T,n> A;
typedef A C;
typedef C D;
static void addProduct(A& a, const T& b, const C& c, const D& d)
{
for (size_t i=0; i<n; ++i)
a.diagonal(i) += b * c.diagonal(i)*d.diagonal(i);
}
static void subtractProduct(A& a, const T& b, const C& c, const D& d)
{
for (size_t i=0; i<n; ++i)
a.diagonal(i) -= b * c.diagonal(i)*d.diagonal(i);
}
};
template<class T, int n>
struct ProductHelper<Dune::FieldMatrix<T,n,n>, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
{
typedef Dune::FieldMatrix<T,n,n> A;
typedef Dune::DiagonalMatrix<T,n> B;
typedef B C;
static void addProduct(A& a, const B& b, const C& c)
{
for (size_t i=0; i<n; ++i)
a[i][i] += b.diagonal(i)*c.diagonal(i);
}
static void subtractProduct(A& a, const B& b, const C& c)
{
for (size_t i=0; i<n; ++i)
a[i][i] -= b.diagonal(i)*c.diagonal(i);
}
};
template<class T, int n>
struct ScaledProductHelper<Dune::FieldMatrix<T,n,n>, T, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
{
typedef Dune::FieldMatrix<T,n,n> A;
typedef Dune::DiagonalMatrix<T,n> C;
typedef C D;
static void addProduct(A& a, const T& b, const C& c, const D& d)
{
for (size_t i=0; i<n; ++i)
a[i][i] += b * c.diagonal(i)*d.diagonal(i);
}
static void subtractProduct(A& a, const T& b, const C& c, const D& d)
{
for (size_t i=0; i<n; ++i)
a[i][i] -= b * c.diagonal(i)*d.diagonal(i);
}
};
/** \brief Specialization for b being a scalar type
*/
template<class A, class B, class C, bool AisScalar, bool CisScalar>
struct ProductHelper<A, B, C, AisScalar, true, CisScalar>
{
static void addProduct(A& a, const B& b, const C& c)
{
a.axpy(b, c);
}
static void subtractProduct(A& a, const B& b, const C& c)
{
a.axpy(-b, c);
}
};
/** \brief Specialization for c being a scalar type
*/
template<class A, class B, class C, class D, bool AisScalar, bool DisScalar>
struct ScaledProductHelper<A, B, C, D, AisScalar, true, DisScalar>
{
static void addProduct(A& a, const B& b, const C& c, const D& d)
{
a.axpy(b * c, d);
}
static void subtractProduct(A& a, const B& b, const C& c, const D& d)
{
a.axpy(-b * c, d);
}
};
template<class A, class B, class C>
struct ProductHelper<A, B, C, true, true, true>
{
static void addProduct(A& a, const B& b, const C& c)
{
a += b*c;
}
static void subtractProduct(A& a, const B& b, const C& c)
{
a -= b*c;
}
};
template<class A, class B, class C, class D>
struct ScaledProductHelper<A, B, C, D, true, true, true>
{
static void addProduct(A& a, const B& b, const C& c, const D& d)
{
a += b*c*d;
}
static void subtractProduct(A& a, const B& b, const C& c, const D& d)
{
a -= b*c*d;
}
};
//! Helper class for computing the cross product
template <class T, int n>
struct CrossProductHelper
{
static Dune::FieldVector<T,n> crossProduct(const Dune::FieldVector<T,n>& a, const Dune::FieldVector<T,n>& b)
{
DUNE_THROW(Dune::Exception, "You can only call crossProduct with dim==3");
}
};
//! Specialisation for n=3
template <class T>
struct CrossProductHelper<T,3>
{
static Dune::FieldVector<T,3> crossProduct(const Dune::FieldVector<T,3>& a, const Dune::FieldVector<T,3>& b)
{
Dune::FieldVector<T,3> r;
r[0] = a[1]*b[2] - a[2]*b[1];
r[1] = a[2]*b[0] - a[0]*b[2];
r[2] = a[0]*b[1] - a[1]*b[0];
return r;
}
};
template<class A, class TransposedA>
struct TransposeHelper
{
static void transpose(const A& a, TransposedA& aT)
{
DUNE_THROW(Dune::Exception, "Not implemented for general matrix types!");
}
};
//! Specialization for Dune::FieldMatrix
template<class T, int n, int m>
struct TransposeHelper<Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,m,n> >
{
static void transpose(const Dune::FieldMatrix<T,n,m>& a, Dune::FieldMatrix<T,m,n>& aT)
{
for (int row = 0; row < m; ++row)
for (int col = 0 ; col < n; ++col)
aT[row][col] = a[col][row];
}
};
//! Specialization for Dune::BCRSMatrix Type
template<class A, class TransposedA>
struct TransposeHelper<Dune::BCRSMatrix<A>, Dune::BCRSMatrix<TransposedA> >
{
static void transpose(const Dune::BCRSMatrix<A>& a, Dune::BCRSMatrix<TransposedA>& aT)
{
Dune::MatrixIndexSet idxSetaT(a.M(),a.N());
typedef typename Dune::BCRSMatrix<A>::ConstColIterator ColIterator;
// add indices into transposed matrix
for (int row = 0; row < a.N(); ++row)
{
ColIterator col = a[row].begin();
ColIterator end = a[row].end();
for ( ; col != end; ++col)
idxSetaT.add(col.index(),row);
}
idxSetaT.exportIdx(aT);
for (int row = 0; row < a.N(); ++row)
{
ColIterator col = a[row].begin();
ColIterator end = a[row].end();
for ( ; col != end; ++col)
TransposeHelper<A,TransposedA>::transpose(*col, aT[col.index()][row]);
}
}
};
/** \brief Add a product to some matrix or vector
*
* This function computes a+=b*c.
*
* This function should tolerate all meaningful
* combinations of scalars, vectors, and matrices.
*
* a,b,c could be matrices with appropriate
* dimensions. But b can also always be a scalar
* represented by a 1-dim vector or a 1 by 1 matrix.
*/
template<class A, class B, class C>
void addProduct(A& a, const B& b, const C& c)
{
ProductHelper<A,B,C,ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar, ScalarTraits<C>::isScalar>::addProduct(a,b,c);
}
/** \brief Subtract a product from some matrix or vector
*
* This function computes a-=b*c.
*
* This function should tolerate all meaningful
* combinations of scalars, vectors, and matrices.
*
* a,b,c could be matrices with appropriate
* dimensions. But b can also always be a scalar
* represented by a 1-dim vector or a 1 by 1 matrix.
*/
template<class A, class B, class C>
void subtractProduct(A& a, const B& b, const C& c)
{
ProductHelper<A,B,C,ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar, ScalarTraits<C>::isScalar>::subtractProduct(a,b,c);
}
//! Compute the cross product of two vectors. Only works for n==3
template<class T, int n>
Dune::FieldVector<T,n> crossProduct(const Dune::FieldVector<T,n>& a, const Dune::FieldVector<T,n>& b)
{
return CrossProductHelper<T,n>::crossProduct(a,b);
}
//! Compute the transposed of a matrix
template <class A, class TransposedA>
void transpose(const A& a, TransposedA& aT)
{
TransposeHelper<A,TransposedA>::transpose(a,aT);
}
/** \brief Add a scaled product to some matrix or vector
*
* This function computes a+=b*c*d.
*
* This function should tolerate all meaningful
* combinations of scalars, vectors, and matrices.
*
* a,c,d could be matrices with appropriate dimensions. But b must
* (currently) and c can also always be a scalar represented by a
* 1-dim vector or a 1 by 1 matrix.
*/
template<class A, class B, class C, class D>
typename Dune::enable_if<ScalarTraits<B>::isScalar, void>::type
addProduct(A& a, const B& b, const C& c, const D& d)
{
ScaledProductHelper<A,B,C,D,ScalarTraits<A>::isScalar, ScalarTraits<C>::isScalar, ScalarTraits<D>::isScalar>::addProduct(a,b,c,d);
}
/** \brief Subtract a scaled product from some matrix or vector
*
* This function computes a-=b*c*d.
*
* This function should tolerate all meaningful
* combinations of scalars, vectors, and matrices.
*
* a,c,d could be matrices with appropriate dimensions. But b must
* (currently) and c can also always be a scalar represented by a
* 1-dim vector or a 1 by 1 matrix.
*/
template<class A, class B, class C, class D>
typename Dune::enable_if<ScalarTraits<B>::isScalar, void>::type
subtractProduct(A& a, const B& b, const C& c, const D& d)
{
ScaledProductHelper<A,B,C,D,ScalarTraits<A>::isScalar, ScalarTraits<C>::isScalar, ScalarTraits<D>::isScalar>::subtractProduct(a,b,c,d);
}
/** \brief Internal helper class for Matrix operations
*
*/
template<class OperatorType, bool isMatrix>
struct OperatorHelper
{
template <class VectorType>
static typename VectorType::field_type
Axy(const OperatorType &A,
const VectorType &x,
const VectorType &y)
{
VectorType tmp(y.size());
tmp = 0.0;
addProduct(tmp, A, x);
return tmp * y;
}
template <class VectorType>
static typename VectorType::field_type
bmAxy(const OperatorType &A, const VectorType &b,
const VectorType &x, const VectorType &y)
{
VectorType tmp = b;
subtractProduct(tmp, A, x);
return tmp * y;
}
};
template<class MatrixType>
struct OperatorHelper<MatrixType, true>
{
template <class VectorType>
static typename VectorType::field_type
Axy(const MatrixType &A,
const VectorType &x,
const VectorType &y)
{
assert(x.N() == A.M());
assert(y.N() == A.N());
typename VectorType::field_type outer = 0.0;
typename VectorType::block_type inner;
typename MatrixType::ConstIterator endi=A.end();
for (typename MatrixType::ConstIterator i=A.begin(); i!=endi; ++i) {
inner = 0.0;
const size_t iindex = i.index();
typename MatrixType::row_type::ConstIterator endj = i->end();
for (typename MatrixType::row_type::ConstIterator j=i->begin();
j!=endj; ++j)
addProduct(inner, *j, x[j.index()]);
outer += inner * y[iindex];
}
return outer;
}
template <class VectorType>
static typename VectorType::field_type
bmAxy(const MatrixType &A, const VectorType &b,
const VectorType &x, const VectorType &y)
{
assert(x.N() == A.M());
assert(y.N() == A.N());
assert(y.N() == b.N());
typename VectorType::field_type outer = 0.0;
typename VectorType::block_type inner;
typename MatrixType::ConstIterator endi=A.end();
for (typename MatrixType::ConstIterator i=A.begin(); i!=endi; ++i) {
const size_t iindex = i.index();
inner = b[iindex];
typename MatrixType::row_type::ConstIterator endj = i->end();
for (typename MatrixType::row_type::ConstIterator j=i->begin();
j!=endj; ++j)
subtractProduct(inner, *j, x[j.index()]);
outer += inner * y[iindex];
}
return outer;
}
};
//! Compute \f$(Ax,y)\f$
template <class OperatorType, class VectorType>
typename VectorType::field_type
Axy(const OperatorType &A,
const VectorType &x,
const VectorType &y)
{
return OperatorHelper<OperatorType,
MatrixTraits<OperatorType>::isMatrix>::Axy(A, x, y);
}
//! Compute \f$(b-Ax,y)\f$
template <class OperatorType, class VectorType>
typename VectorType::field_type
bmAxy(const OperatorType &A,
const VectorType &b,
const VectorType &x,
const VectorType &y)
{
return OperatorHelper<OperatorType,
MatrixTraits<OperatorType>::isMatrix>::bmAxy(A, b, x, y);
}
};
#endif
#ifndef DUNE_BOX_CONSTRAINT_HH // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
#define DUNE_BOX_CONSTRAINT_HH // vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_SOLVERS_COMMON_BOX_CONSTRAINT_HH
#define DUNE_SOLVERS_COMMON_BOX_CONSTRAINT_HH
#include <array>
#include <limits> #include <limits>
#include <dune/common/array.hh>
#include <dune/common/fvector.hh> #include <dune/common/fvector.hh>
#include "interval.hh"
template <class T, int dim> template <class T, int dim>
class BoxConstraint { class BoxConstraint {
public: public:
using Interval = typename Dune::Solvers::Interval<T>;
/** \brief Default constructor, makes the obstacle as large as possible */ /** \brief Default constructor, makes the obstacle as large as possible */
BoxConstraint() { BoxConstraint() {
for (int i=0; i<dim; i++) { clear();
val[2*i] = -std::numeric_limits<T>::max();
val[2*i+1] = std::numeric_limits<T>::max();
}
} }
/** \brief Set obstacle values to +/- std::numeric_limits<T>::max() */ /** \brief Set obstacle values to +/- std::numeric_limits<T>::max() */
void clear() { void clear() {
for (int i=0; i<dim; i++) { for (int i=0; i<dim; i++)
val[2*i] = -std::numeric_limits<T>::max(); val[i] = { -std::numeric_limits<T>::max(),
val[2*i+1] = std::numeric_limits<T>::max(); std::numeric_limits<T>::max() };
}
} }
//! Subtract vector from box //! Subtract vector from box
BoxConstraint<T,dim>& operator-= (const Dune::FieldVector<T, dim>& v) BoxConstraint<T,dim>& operator-= (const Dune::FieldVector<T, dim>& v)
{ {
for (int i=0; i<dim; i++) { for (int i=0; i<dim; i++) {
val[2*i] -= v[i]; val[i][0] -= v[i];
val[2*i+1] -= v[i]; val[i][1] -= v[i];
} }
return *this; return *this;
} }
Interval const &operator[](size_t i) const {
return val[i];
}
Interval &operator[](size_t i) {
return val[i];
}
//! Access to the lower obstacles //! Access to the lower obstacles
T& lower(size_t i) {return val[2*i];} T& lower(size_t i) {return val[i][0];}
//! Const access to the lower obstacles //! Const access to the lower obstacles
const T& lower(size_t i) const {return val[2*i];} const T& lower(size_t i) const {return val[i][0];}
//! Access to the upper obstacles //! Access to the upper obstacles
T& upper(size_t i) {return val[2*i+1];} T& upper(size_t i) {return val[i][1];}
//! Const access to the upper obstacles //! Const access to the upper obstacles
const T& upper(size_t i) const {return val[2*i+1];} const T& upper(size_t i) const {return val[i][1];}
//! Send box constraint to an output stream //! Send box constraint to an output stream
friend friend
std::ostream& operator<< (std::ostream& s, const BoxConstraint<T,dim>& v) std::ostream& operator<< (std::ostream& s, const BoxConstraint<T,dim>& v)
{ {
for (int i=0; i<dim; i++) for (int i=0; i<dim; i++)
s << "Direction: " << i <<", val[0] = " << v.val[2*i] s << "Direction: " << i <<", val[0] = " << v[i][0]
<< ", val[1] = " << v.val[2*i+1] << std::endl; << ", val[1] = " << v[i][1] << std::endl;
return s; return s;
} }
protected: protected:
std::array<Interval, dim> val;
Dune::array<T, 2*dim> val;
}; };
#endif #endif
#ifndef CAN_IGNORE_HH // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
#define CAN_IGNORE_HH // vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_SOLVERS_COMMON_CANIGNORE_HH
#define DUNE_SOLVERS_COMMON_CANIGNORE_HH
#include <cassert>
/** \brief Abstract base class for solvers and iterationsteps that are /** \brief Abstract base class for solvers and iterationsteps that are
* able to ignore degrees of freedom described by a bit field. * able to ignore degrees of freedom described by a bit field.
...@@ -9,19 +13,60 @@ class CanIgnore ...@@ -9,19 +13,60 @@ class CanIgnore
{ {
public: public:
/** \brief Default constructor */ using BitVector = BitVectorType;
CanIgnore()
: ignoreNodes_(0)
{}
/** \brief Virtual destructor. Does NOT delete the bitfield! */ /**
virtual ~CanIgnore() * \brief Default constructor
*/
CanIgnore() = default;
/**
* \brief Constructor from bit vector
*
* This class stores a non-owning pointer to the bit vector.
*/
CanIgnore(const BitVector& i)
: ignoreNodes_(&i)
{} {}
/** \brief A flag for each degree of freedom stating whether the /**
dof should be ignored by the solver */ * \brief Destructor
const BitVectorType* ignoreNodes_; *
* Does NOT delete the bitfield!
*/
~CanIgnore() = default;
/**
* \brief Set bit vector
*
* This class stores a non-owning pointer to the bit vector.
*/
void setIgnore(const BitVector& i)
{
ignoreNodes_ = &i;
}
/**
* \brief Only if this returns true can ignore() safely be called.
*/
bool hasIgnore() const {
return ignoreNodes_ != nullptr;
}
/**
* \brief Const access to bit vector
*/
const BitVector& ignore() const
{
assert(hasIgnore());
return *ignoreNodes_;
}
/**
* \brief A flag for each degree of freedom stating whether the dof should be ignored by the solver
*/
const BitVectorType* ignoreNodes_ = nullptr;
}; };
#endif #endif // DUNE_SOLVERS_COMMON_CANIGNORE_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_SOLVERS_COMMON_COPYORREFERENCE_HH
#define DUNE_SOLVERS_COMMON_COPYORREFERENCE_HH
#include <dune/common/indices.hh>
#include <dune/common/concept.hh>
#include <dune/solvers/common/defaultbitvector.hh>
namespace Dune {
namespace Solvers {
/**
* \brief A wrapper storing either a copy or reference
*
* \tparam T Type to be stored
*
* If T is of reference type, the wrapper stores a reference
* to the passed object. Otherwise it stores a copy.
*
* Internally, for reference type T, this wrapper stores a pointer, which
* results in CopyOrReference being copyable (etc.). This sets this wrapper
* apart from having a pure member type T in a class in particular.
* Furthermore, using this wrapper provides gclarity for the reader that type T
* may be both, of raw type or of reference type.
*
* Usage of this wrapper comes at the price of get()-semantics.
* As an alternative with pointer semantics consider using wrap_own_share but
* mind the dynamic overhead.
*/
template<class T>
class CopyOrReference
{
public:
using Type = std::decay_t<T>;
constexpr CopyOrReference(const Type& other) :
t_(other)
{}
constexpr CopyOrReference(Type&& other) :
t_(std::move(other))
{}
constexpr const Type& get() const
{
return t_;
}
Type& get()
{
return t_;
}
constexpr std::true_type hasCopy() const
{
return {};
}
private:
Type t_;
};
template<class T>
class CopyOrReference<T&>
{
public:
using Type = std::decay_t<T>;
constexpr CopyOrReference(Type& other) :
t_(&other)
{}
constexpr const Type& get() const
{
return *t_;
}
Type& get()
{
return *t_;
}
constexpr std::false_type hasCopy() const
{
return {};
}
private:
Type* t_;
};
template<class T>
class CopyOrReference<const T&>
{
public:
using Type = std::decay_t<T>;
constexpr CopyOrReference(const Type& other) :
t_(&other)
{}
constexpr const Type& get() const
{
return *t_;
}
constexpr std::false_type hasCopy() const
{
return {};
}
private:
const Type* t_;
};
/**
* \brief A wrapper storing either a const copy or const reference
*
* \tparam T Type to be stored
*
* This is an alias for CopyOrReference<TT> where TT
* is either 'const T' if T is a raw type or 'const S&'
* if T is a 'S&'. Notice that in the latter case this
* is different from 'const T' which would be the same
* as T.
*/
template<class T>
using ConstCopyOrReference =
CopyOrReference<
std::conditional_t<
std::is_reference<T>::value,
const std::remove_reference_t<T>&,
const T>>;
} // end namespace Solvers
} // end namespace Dune
#endif // DUNE_SOLVERS_COMMON_COPYORREFERENCE_HH
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=8 sw=2 sts=2:
#ifndef DUNE_SOLVERS_COMMON_DEFAULTBITVECTOR_HH
#define DUNE_SOLVERS_COMMON_DEFAULTBITVECTOR_HH
#include <dune/common/bitsetvector.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/multitypeblockvector.hh>
#include <dune/common/tuplevector.hh>
namespace Dune {
namespace Solvers {
namespace Imp {
/** \brief Construct a bitvector type with a nesting that matches the one of a given vector
*
* The default implementation implements the case where the vector is a number type
* (to end recursion), and returns 'char' in such cases.
*
* \tparam Vector The vector types whose nesting shall be matched
*/
template<class Vector>
struct DefaultBitVector
{
static_assert(IsNumber<Vector>::value, "No DefaultBitVector known for this type.");
using type = char;
};
template<class T, int i>
struct DefaultBitVector<FieldVector<T,i>>
{
using type = std::bitset<i>;
};
template<class T, class A>
struct DefaultBitVector<BlockVector<T,A>>
{
using type = typename std::vector<typename DefaultBitVector<T>::type>;
};
template<class T, class A>
struct DefaultBitVector<std::vector<T,A>>
{
using type = typename std::vector<typename DefaultBitVector<T>::type>;
};
template<class T, int i, class A>
struct DefaultBitVector<BlockVector<FieldVector<T,i>, A>>
{
using type = BitSetVector<i>;
};
template<typename... Args>
struct DefaultBitVector<MultiTypeBlockVector<Args...> >
{
using type = TupleVector<typename DefaultBitVector<Args>::type...>;
};
} // end namespace Imp
/**
* \brief Type alias for a bit vector that fits to T
*
* \tparam T A nested container type
*
* This alias provides a default choice of a bit vector
* that has the same nested structure as T.
*/
template<class T>
using DefaultBitVector_t = typename Imp::DefaultBitVector<T>::type;
} // end namespace Solvers
} // end namespace Dune
#endif // DUNE_SOLVERS_COMMON_DEFAULTBITVECTOR_HH
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef GENERIC_VECTOR_TOOL_HH #ifndef GENERIC_VECTOR_TOOL_HH
#define GENERIC_VECTOR_TOOL_HH #define GENERIC_VECTOR_TOOL_HH
...@@ -11,13 +13,17 @@ ...@@ -11,13 +13,17 @@
#include <dune/common/diagonalmatrix.hh> #include <dune/common/diagonalmatrix.hh>
#include <dune/common/fvector.hh> #include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh> #include <dune/common/fmatrix.hh>
#include <dune/common/static_assert.hh>
#include <dune/common/diagonalmatrix.hh> #include <dune/common/diagonalmatrix.hh>
#include <dune/common/indices.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/istl/scaledidmatrix.hh> #include <dune/istl/scaledidmatrix.hh>
#include <dune/istl/bcrsmatrix.hh> #include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh> #include <dune/istl/bvector.hh>
#include <dune/matrix-vector/genericvectortools.hh>
#include <dune/matrix-vector/resize.hh>
#include <dune/solvers/operators/sumoperator.hh> #include <dune/solvers/operators/sumoperator.hh>
/** \brief Various tools for working with istl vectors of arbitrary nesting depth /** \brief Various tools for working with istl vectors of arbitrary nesting depth
...@@ -25,114 +31,42 @@ ...@@ -25,114 +31,42 @@
struct GenericVector struct GenericVector
{ {
//! Write vector to given stream //! Write vector to given stream
template <class VectorType> template <class VectorType>
[[deprecated("Please use Dune::MatrixVector::Generic::writeBinary instead.")]]
static void writeBinary(std::ostream& s, const VectorType& v) static void writeBinary(std::ostream& s, const VectorType& v)
{ {
typename VectorType::const_iterator it = v.begin(); Dune::MatrixVector::Generic::writeBinary(s, v);
typename VectorType::const_iterator end = v.end();
for(; it!=end; ++it)
GenericVector::writeBinary(s, *it);
}
template <class field_type, int n>
static void writeBinary(std::ostream& s, const Dune::FieldVector<field_type,n>& v)
{
typedef typename Dune::FieldVector<field_type,n> VectorType;
typename VectorType::const_iterator it = v.begin();
typename VectorType::const_iterator end = v.end();
for(; it!=end; ++it)
s.write(reinterpret_cast<const char*>(&(*it)), sizeof(field_type));
} }
//! Read vector from a given stream //! Read vector from a given stream
template <class VectorType> template <class VectorType>
[[deprecated("Please use Dune::MatrixVector::Generic::readBinary instead.")]]
static void readBinary(std::istream& s, VectorType& v) static void readBinary(std::istream& s, VectorType& v)
{ {
typename VectorType::iterator it = v.begin(); Dune::MatrixVector::Generic::readBinary(s, v);
typename VectorType::iterator end = v.end();
for(; it!=end; ++it)
GenericVector::readBinary(s, *it);
} }
template <class field_type, int n> //! Resize vector recursively to size of given vector/matrix
static void readBinary(std::istream& s, Dune::FieldVector<field_type,n>& v)
{
typedef typename Dune::FieldVector<field_type,n> VectorType;
typename VectorType::iterator it = v.begin();
typename VectorType::iterator end = v.end();
for(; it!=end; ++it)
s.read(reinterpret_cast<char*>(&(*it)), sizeof(field_type));
}
//! Resize vector recursivly to size of given vector/matrix
template <class VectorTypeA, class VectorTypeB> template <class VectorTypeA, class VectorTypeB>
[[deprecated("Please use Dune::MatrixVector::resize instead.")]]
static void resize(VectorTypeA& a, const VectorTypeB& b) static void resize(VectorTypeA& a, const VectorTypeB& b)
{ {
a.resize(b.size()); Dune::MatrixVector::resize(a, b);
typename VectorTypeB::const_iterator it = b.begin();
typename VectorTypeB::const_iterator end = b.end();
for(; it!=end; ++it)
GenericVector::resize(a[it.index()], *it);
}
template <class VectorTypeA, class T, class A>
static void resize(VectorTypeA& a, const Dune::BCRSMatrix<T,A>& b)
{
a.resize(b.N());
typedef typename Dune::BCRSMatrix<T,A> VectorTypeB;
typename VectorTypeB::const_iterator it = b.begin();
typename VectorTypeB::const_iterator end = b.end();
for(; it!=end; ++it)
GenericVector::resize(a[it.index()], *(it->begin()));
} }
template <class VectorTypeA, class SparseMatrixType, class LowRankMatrixType> //! Set vector to zero at indices that are true in bitvector recursively
static void resize(VectorTypeA& a, const SumOperator<SparseMatrixType, LowRankMatrixType>& b)
{
a.resize(b.N());
typename SparseMatrixType::const_iterator it = b.sparseMatrix().begin();
typename SparseMatrixType::const_iterator end = b.sparseMatrix().end();
for(; it!=end; ++it)
GenericVector::resize(a[it.index()], *(it->begin()));
}
template <class field_type, int n, class VectorTypeB>
static void resize(Dune::FieldVector<field_type,n>& a, const VectorTypeB& b)
{}
template <int n, class VectorTypeB>
static void resize(std::bitset<n>& a, const VectorTypeB& b)
{}
//! Set vector to zero at indices that are true in bitvector recursivly
template <class VectorType, class BitVectorType> template <class VectorType, class BitVectorType>
[[deprecated("Please use Dune::MatrixVector::Generic::truncate instead.")]]
static void truncate(VectorType& v, const BitVectorType& tr) static void truncate(VectorType& v, const BitVectorType& tr)
{ {
typename VectorType::iterator it = v.begin(); Dune::MatrixVector::Generic::truncate(v, tr);
typename VectorType::iterator end = v.end();
for(; it!=end; ++it)
GenericVector::truncate(*it, tr[it.index()]);
}
template <class field_type, int n, class BitVectorType>
static void truncate(Dune::FieldVector<field_type,n>& v, const BitVectorType& tr)
{
typedef typename Dune::FieldVector<field_type,n> VectorType;
typename VectorType::iterator it = v.begin();
typename VectorType::iterator end = v.end();
for(; it!=end; ++it)
if (tr[it.index()])
*it = 0;
} }
//! weave two vector blocks into each other e.g. [[u1...un][w1...wn]] --> [[u1 w1]...[un wn]] //! weave two vector blocks into each other e.g. [[u1...un][w1...wn]] --> [[u1 w1]...[un wn]]
template <class LVectorType, class RVectorType> template <class LVectorType, class RVectorType>
static void interlace(const LVectorType& lvec, RVectorType& rvec) static void interlace([[maybe_unused]] const LVectorType& lvec, [[maybe_unused]] RVectorType& rvec)
{ {
DUNE_THROW(Dune::NotImplemented,"GenericVector::interlace not implemented for given VectorTypes"); DUNE_THROW(Dune::NotImplemented,"GenericVector::interlace not implemented for given VectorTypes");
} }
...@@ -183,7 +117,7 @@ struct GenericVector ...@@ -183,7 +117,7 @@ struct GenericVector
template <class LVectorBlock, class RVectorBlock> template <class LVectorBlock, class RVectorBlock>
static void interlace(const Dune::BlockVector<Dune::BlockVector<LVectorBlock> >& lvec, Dune::BlockVector<RVectorBlock>& rvec) static void interlace(const Dune::BlockVector<Dune::BlockVector<LVectorBlock> >& lvec, Dune::BlockVector<RVectorBlock>& rvec)
{ {
dune_static_assert(RVectorBlock::dimension % LVectorBlock::dimension == 0, "Block sizes don't match."); static_assert(RVectorBlock::dimension % LVectorBlock::dimension == 0, "Block sizes don't match.");
rvec.resize(lvec[0].size()); rvec.resize(lvec[0].size());
rvec = 0.0; rvec = 0.0;
...@@ -198,7 +132,7 @@ struct GenericVector ...@@ -198,7 +132,7 @@ struct GenericVector
//! unweave two vectors previously interlaced e.g. [[u1 w1]...[un wn]] --> [[u1...un][w1...wn]] //! unweave two vectors previously interlaced e.g. [[u1 w1]...[un wn]] --> [[u1...un][w1...wn]]
template <class LVectorType, class RVectorType> template <class LVectorType, class RVectorType>
static void deinterlace(const LVectorType& lvec, RVectorType& rvec) static void deinterlace([[maybe_unused]] const LVectorType& lvec, [[maybe_unused]] RVectorType& rvec)
{ {
DUNE_THROW(Dune::NotImplemented,"GenericVector::deinterlace not implemented for given VectorTypes"); DUNE_THROW(Dune::NotImplemented,"GenericVector::deinterlace not implemented for given VectorTypes");
} }
...@@ -206,7 +140,7 @@ struct GenericVector ...@@ -206,7 +140,7 @@ struct GenericVector
template <class LVectorBlock, class RVectorBlock> template <class LVectorBlock, class RVectorBlock>
static void deinterlace(const Dune::BlockVector<LVectorBlock>& lvec, Dune::BlockVector<Dune::BlockVector<RVectorBlock> >& rvec) static void deinterlace(const Dune::BlockVector<LVectorBlock>& lvec, Dune::BlockVector<Dune::BlockVector<RVectorBlock> >& rvec)
{ {
dune_static_assert(LVectorBlock::dimension % RVectorBlock::dimension == 0, "Block sizes don't match."); static_assert(LVectorBlock::dimension % RVectorBlock::dimension == 0, "Block sizes don't match.");
rvec.resize(LVectorBlock::dimension / RVectorBlock::dimension); rvec.resize(LVectorBlock::dimension / RVectorBlock::dimension);
for (size_t k=0; k<rvec.size(); ++k) for (size_t k=0; k<rvec.size(); ++k)
...@@ -249,7 +183,7 @@ struct GenericVector ...@@ -249,7 +183,7 @@ struct GenericVector
struct GenericMatrix struct GenericMatrix
{ {
//! Set matrix to zero at indices that are true in bitvector recursivly //! Set matrix to zero at indices that are true in bitvector recursively
template <class MatrixType, class BitVectorTypeR, class BitVectorTypeC> template <class MatrixType, class BitVectorTypeR, class BitVectorTypeC>
static void truncate(MatrixType& mat, const BitVectorTypeR& trows, const BitVectorTypeC& tcols, bool setTruncatedDiagonalOne) static void truncate(MatrixType& mat, const BitVectorTypeR& trows, const BitVectorTypeC& tcols, bool setTruncatedDiagonalOne)
{ {
...@@ -274,7 +208,7 @@ struct GenericMatrix ...@@ -274,7 +208,7 @@ struct GenericMatrix
template <class field_type, int n, int m, class BitVectorTypeR, class BitVectorTypeC> template <class field_type, int n, int m, class BitVectorTypeR, class BitVectorTypeC>
static void truncate(Dune::FieldMatrix<field_type,n,m>& mat, const BitVectorTypeR& trows, const BitVectorTypeC& tcols, bool setTruncatedDiagonalOne) static void truncate(Dune::FieldMatrix<field_type,n,m>& mat, const BitVectorTypeR& trows, const BitVectorTypeC& tcols, bool setTruncatedDiagonalOne)
{ {
// dune_static_assert(BitVectorTypeR.size()==n and BitVectorTypeC.size()==m,"BitVector length doesn't match matrix size."); // static_assert(BitVectorTypeR.size()==n and BitVectorTypeC.size()==m,"BitVector length doesn't match matrix size.");
assert(trows.size()==n and tcols.size()==m); assert(trows.size()==n and tcols.size()==m);
for(size_t row=0; row<n; ++row) for(size_t row=0; row<n; ++row)
......
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_SOLVERS_COMMON_INTERVAL_HH #ifndef DUNE_SOLVERS_COMMON_INTERVAL_HH
#define DUNE_SOLVERS_COMMON_INTERVAL_HH #define DUNE_SOLVERS_COMMON_INTERVAL_HH
#include <array>
#include <iostream> #include <iostream>
#include <algorithm> #include <algorithm>
#include <dune/common/array.hh>
namespace Dune { namespace Dune {
namespace Solvers { namespace Solvers {
...@@ -19,27 +18,69 @@ template <class field_type> ...@@ -19,27 +18,69 @@ template <class field_type>
class Interval class Interval
{ {
public: public:
/** \brief Default constructor */
Interval()
{}
/** \brief Construct from an initializer list */
Interval(std::initializer_list<field_type> const &input)
{
std::copy(input.begin(), input.end(), data_.begin());
}
/** \brief Construct from a scalar */
Interval(const field_type& s)
{
data_[0] = data_[1] = s;
}
/** \brief Array-like access /** \brief Array-like access
*/ */
field_type& operator[](int i) field_type& operator[](int i)
{ {
return data_[i]; return data_[i];
} }
/** \brief Const array-like access /** \brief Const array-like access
*/ */
const field_type& operator[](int i) const const field_type& operator[](int i) const
{ {
return data_[i]; return data_[i];
} }
/** \brief Addition
*/
Interval& operator+=(const Interval& other)
{
data_[0] += other.data_[0];
data_[1] += other.data_[1];
return *this;
}
/** \brief Scalar multiplication
*/
Interval& operator*=(const field_type& c)
{
data_[0] *= c;
data_[1] *= c;
// swap if multiplication with negative number
if ( c < 0.0 )
std::swap(data_[0],data_[1]);
return *this;
}
/** \brief Project a scalar onto the interval /** \brief Project a scalar onto the interval
*/ */
field_type projectIn(const field_type& x) const field_type projectIn(const field_type& x) const
{ {
// NB: C++17 has std::clamp(x, data_[0], data_[1]).
// Note, however, that we need to guarantee neither
// data_[0] nor data_[1] is NaN then.
return std::max(std::min(x,data_[1]), data_[0]); return std::max(std::min(x,data_[1]), data_[0]);
} }
/** \brief Fast projection onto the interval if you know that your value /** \brief Fast projection onto the interval if you know that your value
* is smaller than your upper bound * is smaller than your upper bound
*/ */
...@@ -65,11 +106,27 @@ public: ...@@ -65,11 +106,27 @@ public:
}; };
private: private:
/** \brief The actual data */ /** \brief The actual data */
Dune::array<field_type,2> data_; std::array<field_type,2> data_;
}; };
//! Scalar multiplication from the left
template <class field_type>
Interval<field_type> operator*(const field_type& c, const Interval<field_type>& i)
{
auto ret = i;
return ret *= c;
}
//! Scalar multiplication from the right
template <class field_type>
Interval<field_type> operator*(const Interval<field_type>& i, const field_type& c)
{
return c * i;
}
//! Output operator for Interval //! Output operator for Interval
template <class field_type> template <class field_type>
inline std::ostream& operator<< (std::ostream& s, const Interval<field_type>& interval) inline std::ostream& operator<< (std::ostream& s, const Interval<field_type>& interval)
...@@ -81,5 +138,5 @@ inline std::ostream& operator<< (std::ostream& s, const Interval<field_type>& in ...@@ -81,5 +138,5 @@ inline std::ostream& operator<< (std::ostream& s, const Interval<field_type>& in
} // namespace Solvers } // namespace Solvers
} // namespace Dune } // namespace Dune
#endif #endif
#ifndef DUNE_NUMPROC_HH // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
#define DUNE_NUMPROC_HH // vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_SOLVERS_COMMON_NUMPROC_HH
#define DUNE_SOLVERS_COMMON_NUMPROC_HH
#include <dune/common/exceptions.hh> #include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
/** \brief Exception thrown by solvers */ /** \brief Exception thrown by solvers */
class SolverError : public Dune::Exception {}; class SolverError : public Dune::Exception {};
/** \brief Base class for numerical procedures */ /** \brief Base class for numerical procedures */
class NumProc class NumProc
{ {
public: public:
/** \brief Different levels of verbosity */
enum VerbosityMode {QUIET, REDUCED, FULL};
/** \brief Different levels of verbosity */ /** \brief Constructor, per default uses FULL. */
enum VerbosityMode {QUIET, REDUCED, FULL}; NumProc(VerbosityMode verbosity = FULL)
: verbosity_(verbosity)
{}
NumProc() : verbosity_(FULL) {} /** \brief Construct from parameter file. */
NumProc(const Dune::ParameterTree& config)
{
setVerbosity(config);
}
NumProc(VerbosityMode verbosity) /** \brief Set verbosity from parameter file. */
: verbosity_(verbosity) void setVerbosity(const Dune::ParameterTree& config) {
{} verbosity_ = config.get<VerbosityMode>("verbosity");
}
/** \brief Controls the amount of screen output of a numproc */ /** \brief Set the verbosity level */
VerbosityMode verbosity_; void setVerbosity(VerbosityMode verbosity) { verbosity_ = verbosity; }
}; /** \brief Get the verbosity level */
const VerbosityMode& getVerbosity() const { return verbosity_; }
/** \brief Check if verbosity is quiet. */
bool isQuiet() const { return verbosity_ == QUIET; }
/** \brief Check if verbosity is reduced. */
bool isReduced() const { return verbosity_ == REDUCED; }
/** \brief Check if verbosity is full. */
bool isFull() const { return verbosity_ == FULL; }
protected:
/** \brief Controls the amount of screen output of a numproc */
VerbosityMode verbosity_;
};
inline
std::istream& operator>>(std::istream &lhs, NumProc::VerbosityMode &e)
{
std::string s;
lhs >> s;
if (s == "full" || s == "2")
e = NumProc::FULL;
else if (s == "reduced" || s == "1")
e = NumProc::REDUCED;
else if (s == "quiet" || s == "0")
e = NumProc::QUIET;
else
lhs.setstate(std::ios_base::failbit);
return lhs;
}
#endif #endif
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_PERMUTATIONMANAGER_HH #ifndef DUNE_PERMUTATIONMANAGER_HH
#define DUNE_PERMUTATIONMANAGER_HH #define DUNE_PERMUTATIONMANAGER_HH
...@@ -13,7 +15,7 @@ ...@@ -13,7 +15,7 @@
//! Idea of implementation: define an ordering for the nodes (with respect to the axis of anisotropy) //! Idea of implementation: define an ordering for the nodes (with respect to the axis of anisotropy)
//! and then use a merge sort algorithm to renumerate the nodes along the lines which are parallel to this axis. //! and then use a merge sort algorithm to renumerate the nodes along the lines which are parallel to this axis.
//! to do so, we introduce a 'NodeList' class. //! to do so, we introduce a 'NodeList' class.
//! everything is handeled by the PermutationManager //! everything is handled by the PermutationManager
// type of nodes: // type of nodes:
...@@ -83,7 +85,7 @@ public: ...@@ -83,7 +85,7 @@ public:
}; };
template <class CoordinateImp, int dim_domain> template <class CoordinateImp, int dim_domain>
class AxisComparison class AxisComparison
{ {
...@@ -114,7 +116,7 @@ public: ...@@ -114,7 +116,7 @@ public:
while ( current_axis >= 0 ) while ( current_axis >= 0 )
{ {
if ( current_axis == primaryAxis_ ) if ( current_axis == primaryAxis_ )
current_axis -= 1; current_axis -= 1;
if ( current_axis < 0 ) if ( current_axis < 0 )
break; break;
...@@ -149,12 +151,12 @@ public: ...@@ -149,12 +151,12 @@ public:
// the information about the axis of anisotropy is in NodeList, since this list contains a methode "less_equal", which describes the ordering of the nodes. // the information about the axis of anisotropy is in NodeList, since this list contains a method "less_equal", which describes the ordering of the nodes.
// this ordering depends on the axis of anisotropy // this ordering depends on the axis of anisotropy
template< class GridViewImp > template< class GridViewImp >
class PermutationManager{ class PermutationManager{
public: public:
...@@ -166,15 +168,16 @@ public: ...@@ -166,15 +168,16 @@ public:
#endif #endif
static const int dim = GridType::dimension; static const int dim = GridType::dimension;
static const int dimworld = GridType::dimensionworld;
typedef typename GridViewType :: template Codim<dim> :: Iterator NodeIterator; typedef typename GridViewType :: template Codim<dim> :: Iterator NodeIterator;
typedef typename GridType :: template Codim<dim> :: Geometry NodeGeometry; typedef typename GridType :: template Codim<dim> :: Geometry NodeGeometry;
typedef Dune::FieldVector<double, dim> DomainPointType; // type of point in the domain typedef Dune::FieldVector<double, dimworld> DomainPointType; // type of point in the domain
typedef Node< DomainPointType , dim > NodeType; // type of nodes typedef Dune::FieldVector<typename GridType::ctype, dimworld> CoordinateType;
typedef Node< CoordinateType, dim > NodeType; // type of nodes
typedef std::vector< NodeType > NodeListType; // array type list of nodes typedef std::vector< NodeType > NodeListType; // array type list of nodes
typedef Dune::FieldVector<typename GridType::ctype, GridType::dimensionworld> CoordinateType;
const GridViewType gridView_; const GridViewType gridView_;
...@@ -224,7 +227,7 @@ private: ...@@ -224,7 +227,7 @@ private:
{ {
const NodeGeometry& node_geometry = node->geometry(); const NodeGeometry& node_geometry = node->geometry();
DomainPointType global_coordinate_node = node_geometry.corner(0); CoordinateType global_coordinate_node = node_geometry.corner(0);
node_list_[idSet.index(*node)].set_node_index( idSet.index(*node) ); node_list_[idSet.index(*node)].set_node_index( idSet.index(*node) );
node_list_[idSet.index(*node)].set_permuted_node_index( idSet.index(*node) ); node_list_[idSet.index(*node)].set_permuted_node_index( idSet.index(*node) );
...@@ -366,7 +369,7 @@ public: ...@@ -366,7 +369,7 @@ public:
if ( block_number < number_of_blocks_ ) if ( block_number < number_of_blocks_ )
{ return block_sizes[block_number]; } { return block_sizes[block_number]; }
else else
{ std :: cout << "Request for block number " << block_number << ", but there are only " << number_of_blocks_ << " blocks." << std :: endl; { std :: cout << "Request for block number " << block_number << ", but there are only " << number_of_blocks_ << " blocks." << std :: endl;
std :: cout << "Note that the numeration of the blocks starts with 0." << std :: endl; std :: cout << "Note that the numeration of the blocks starts with 0." << std :: endl;
abort(); abort();
} }
...@@ -427,9 +430,9 @@ public: ...@@ -427,9 +430,9 @@ public:
for( int i = 0; i < node_list_.size(); ++i ) for( int i = 0; i < node_list_.size(); ++i )
{ {
std::cout << "Original node[" << get_inverse_permuted_index( i ) << "] = (" << node_list_[i].get_global_coordinate() << ") recieved permuted index " << i << std::endl; std::cout << "Original node[" << get_inverse_permuted_index( i ) << "] = (" << node_list_[i].get_global_coordinate() << ") received permuted index " << i << std::endl;
// alternative: // alternative:
//std::cout << "Original node[" << get_inverse_permuted_index( i ) << "] = (" << original_node_list_[get_inverse_permuted_index( i )].get_global_coordinate() << ") recieved permuted index " << i << std::endl; //std::cout << "Original node[" << get_inverse_permuted_index( i ) << "] = (" << original_node_list_[get_inverse_permuted_index( i )].get_global_coordinate() << ") received permuted index " << i << std::endl;
} }
#if 0 #if 0
......
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef PRECONDITIONER_HH #ifndef PRECONDITIONER_HH
#define PRECONDITIONER_HH #define PRECONDITIONER_HH
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_SOLVERS_COMMON_RESIZE_HH
#define DUNE_SOLVERS_COMMON_RESIZE_HH
#include <dune/common/indices.hh>
#include <dune/common/concept.hh>
#include <dune/common/typetraits.hh>
#include <dune/solvers/common/defaultbitvector.hh>
namespace Dune {
namespace Solvers {
namespace Concept {
struct HasResize
{
template<class C>
auto require(C&& c) -> decltype(
c.resize(0)
);
};
} // namespace Concept
namespace Impl {
template<int size, class Alloc, class Vector>
void resizeInitialize(Dune::BitSetVector<size, Alloc>& x, const Vector& y, bool value)
{
x.resize(y.size());
if (value)
x.setAll();
else
x.unsetAll();
}
template<class TargetVector, class Vector, class Value,
std::enable_if_t<not Dune::IsIndexable<TargetVector, std::integral_constant<std::size_t,0> >::value, int> = 0>
void resizeInitialize(TargetVector& x, const Vector& y, const Value& value)
{
x = value;
}
template<class TargetVector, class Vector, class Value,
std::enable_if_t<Dune::IsIndexable<TargetVector, std::integral_constant<std::size_t,0> >::value, int> = 0>
void resizeInitialize(TargetVector& x, const Vector& y, const Value& value)
{
namespace H = Dune::Hybrid;
auto size = H::size(y);
H::ifElse(models<Concept::HasResize, TargetVector>(), [&](auto&& id) {
id(x).resize(size);
}, [&](auto&& id) {
if (H::size(x) != size)
DUNE_THROW(RangeError, "Can't resize statically sized vector of size " << H::size(x) << " to size " << size);
});
H::forEach(H::integralRange(size), [&](auto&& i) {
resizeInitialize(x[i], y[i], value);
});
}
} // namespace Impl
/**
* \brief Resize and initialization vector to match size of given vector
*
* \param x Vector to resize
* \param y Model for resizing
* \param value Value to use for initialization
*
* This will resize the given vector x to match
* the size of the given vector and assign the given
* value to it.
*/
template<class TargetVector, class Vector, class Value>
void resizeInitialize(TargetVector& x, const Vector& y, Value&& value)
{
Impl::resizeInitialize(x, y, value);
}
/**
* \brief Resize and initialization vector to match size of given vector
*
* \param x Vector to resize
* \param y Model for resizing
*
* This will resize the given vector x to match
* the size of the given vector and initialize
* all entries with zero.
*/
template<class TargetVector, class Vector>
void resizeInitializeZero(TargetVector& x, const Vector& y)
{
resizeInitialize(x, y, 0);
}
} // end namespace Solvers
} // end namespace Dune
#endif // DUNE_SOLVERS_COMMON_RESIZE_HH
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef STATIC_MATRIX_TOOL_HH #ifndef STATIC_MATRIX_TOOL_HH
#define STATIC_MATRIX_TOOL_HH #define STATIC_MATRIX_TOOL_HH
#include "dune/common/diagonalmatrix.hh"
#include "dune/common/fmatrix.hh" #include "dune/common/fmatrix.hh"
#include "dune/common/diagonalmatrix.hh"
#include "dune/istl/scaledidmatrix.hh" #include "dune/istl/scaledidmatrix.hh"
#include "dune/istl/bcrsmatrix.hh" #include "dune/istl/bcrsmatrix.hh"
#include "dune/istl/matrixindexset.hh" #include "dune/istl/matrixindexset.hh"
class StaticMatrix #include <dune/matrix-vector/addtodiagonal.hh>
{ #include <dune/matrix-vector/axpy.hh>
public: #include <dune/matrix-vector/ldlt.hh>
#include <dune/matrix-vector/promote.hh>
#include <dune/matrix-vector/scalartraits.hh>
#include <dune/matrix-vector/transformmatrix.hh>
namespace StaticMatrix
{
// type promotion (smallest matrix that can hold the sum of two matrices ******************* // type promotion (smallest matrix that can hold the sum of two matrices *******************
template<class MatrixA, class MatrixB> template<class MatrixA, class MatrixB>
struct Promote struct Promote
{ {
typedef Dune::FieldMatrix<typename MatrixA::field_type, MatrixA::rows, MatrixA::cols> Type; typedef typename Dune::MatrixVector::Promote<MatrixA, MatrixB>::Type Type;
};
template<class Matrix>
struct Promote<Matrix, Matrix>
{
typedef Matrix Type;
};
template<typename FieldType, int n>
struct Promote<Dune::FieldMatrix<FieldType,n,n>, Dune::DiagonalMatrix<FieldType,n> >
{
typedef Dune::FieldMatrix<FieldType,n,n> Type;
}; };
template<typename FieldType, int n>
struct Promote<Dune::DiagonalMatrix<FieldType,n>, Dune::FieldMatrix<FieldType,n,n> >
{
typedef Dune::FieldMatrix<FieldType,n,n> Type;
};
template<typename FieldType, int n>
struct Promote<Dune::DiagonalMatrix<FieldType,n>, Dune::ScaledIdentityMatrix<FieldType,n> >
{
typedef Dune::DiagonalMatrix<FieldType,n> Type;
};
template<typename FieldType, int n>
struct Promote<Dune::ScaledIdentityMatrix<FieldType,n>, Dune::DiagonalMatrix<FieldType,n> >
{
typedef Dune::DiagonalMatrix<FieldType,n> Type;
};
// add scalar to matrix diagonal *********************************************************** // add scalar to matrix diagonal ***********************************************************
template<class Matrix> template<class Matrix>
static void addToDiagonal(Matrix& x, const typename Matrix::field_type a) static void addToDiagonal(Matrix& x, const typename Matrix::field_type a)
{ {
for(int i=0; i<Matrix::rows; ++i) Dune::MatrixVector::addToDiagonal(x, a);
x[i][i] += a;
}
template <typename FieldType, int n>
// static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const FieldType a)
// the commented line should NOT be used as it may lead to ambiguities in which case the general method above will be used.
// an example is the line taken from massassembler.hh (line 83 at the time):
// StaticMatrix::addToDiagonal(localMatrix[i][i], values[i] * zi);
// where values[i] is of type FieldVector<FieldType,1> and zi is a double. This call then does not exactly fit this specialization (without the implicit cast of FV<FT,1>)
// and hence some wild template voodoo decides which of the templates is to be taken - in this case with a gcc 4.4.5 that was the one above leading to a wrong factor n
// in the diagonal value
static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const typename Dune::ScaledIdentityMatrix<FieldType,n>::field_type a)
{
x.scalar() += a;
}
// add scalar times matrix to matrix *******************************************************
template<class MatrixA, class MatrixB>
static void axpy(MatrixA& A, const typename MatrixA::field_type a, const MatrixB& B)
{
for(size_t i=0; i<B.N(); ++i)
{
typename MatrixB::row_type::ConstIterator it=B[i].begin();
typename MatrixB::row_type::ConstIterator end=B[i].end();
for(; it!=end; ++it)
A[i][it.index()] += a * (*it);
}
}
// template<class MatrixA, class MatrixB>
// static void axpy(MatrixA& x, const typename MatrixA::field_type a, const MatrixB& y)
// {
// for(int i=0; i<MatrixA::rows; ++i)
// for(int j=0; j<MatrixA::cols; ++j)
// x[i][j] += a * y[i][j];
// }
template <typename FieldType, int n, class MatrixB>
static void axpy(Dune::ScaledIdentityMatrix<FieldType,n>& x, const typename Dune::ScaledIdentityMatrix<FieldType,n>::field_type a, const MatrixB& y)
{
x.scalar() += a * y[0][0];
}
template <typename FieldType, int n, class MatrixB>
static void axpy(Dune::DiagonalMatrix<FieldType,n>& x, const typename Dune::DiagonalMatrix<FieldType,n>::field_type a, const MatrixB& y)
{
for(int i=0; i<n; ++i)
x.diagonal()[i] += a * y[i][i];
}
template <typename FieldType, int n, typename Scalar>
static void axpy(Dune::FieldVector<FieldType,n>& x, const Scalar a, const Dune::FieldVector<FieldType,n>& y)
{
x.axpy(a, y);
} }
// add transformed matrix A += T1^t*B*T2 ******************************************************
// add transformed matrix X = A^T*Y*B ******************************************************
template<class MatrixA, class MatrixB, class TransformationMatrix> template<class MatrixA, class MatrixB, class TransformationMatrix>
static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2) void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{
for (size_t i=0; i<A.N(); ++i)
{
typename MatrixA::row_type::ConstIterator jIt=A[i].begin();
typename MatrixA::row_type::ConstIterator jEnd=A[i].end();
for(; jIt!=jEnd; ++jIt)
{
for (size_t k=0; k<B.N(); ++k)
{
if (T1[k][i]!=0)
{
typename MatrixB::row_type::ConstIterator lIt=B[k].begin();
typename MatrixB::row_type::ConstIterator lEnd=B[k].end();
for(; lIt!=lEnd; ++lIt)
A[i][jIt.index()] += T1[k][i] * B[k][lIt.index()] * T2[lIt.index()][jIt.index()];
}
}
}
}
}
template<class K1, class K2, class K3, int n, int m>
static void addTransformedMatrix(
typename Dune::FieldMatrix<K1, m, m>& A,
const typename Dune::FieldMatrix<K2, n, m>& T1,
const typename Dune::FieldMatrix<K3, n, n>& B,
const typename Dune::FieldMatrix<K2, n, m>& T2)
{
typename Dune::FieldMatrix<K1, m, n> T2transposedB;
T2transposedB = 0;
for (size_t i=0; i<A.N(); ++i)
{
for (size_t k=0; k<B.N(); ++k)
{
if (T1[k][i]!=0)
for (size_t l=0; l<B.M(); ++l)
T2transposedB[i][l] += T1[k][i] * B[k][l];
}
}
for (size_t k=0; k<T2.N(); ++k)
{
for (size_t l=0; l<T2.M(); ++l)
{
if (T2[k][l]!=0)
for (size_t i=0; i<A.N(); ++i)
A[i][l] += T2transposedB[i][k] * T2[k][l];
}
}
}
template<class MatrixA, class MatrixB>
static void addTransformedMatrix(MatrixA& X, const double& A, const MatrixB& Y, const double& B)
{ {
axpy(X, A*B, Y); Dune::MatrixVector::addTransformedMatrix(A,T1,B,T2);
} }
template<class MatrixA, typename FieldType, int n, class TransformationMatrix> template<class MatrixA, class MatrixB, class TransformationMatrix>
static void addTransformedMatrix(MatrixA& X, const TransformationMatrix& A, const Dune::DiagonalMatrix<FieldType,n>& Y, const TransformationMatrix& B) void transformMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{
for (int i=0; i<MatrixA::rows; ++i)
for (int j=0; j<MatrixA::cols; ++j)
for (int k=0; k<n; k++)
X[i][j] += A[k][i] * Y.diagonal()[k] * B[k][j];
}
template<class MatrixA, typename FieldType, int n, class TransformationMatrix>
static void addTransformedMatrix(MatrixA& X, const TransformationMatrix& A, const Dune::ScaledIdentityMatrix<FieldType,n>& Y, const TransformationMatrix& B)
{
for (int i=0; i<MatrixA::rows; ++i)
for (int j=0; j<MatrixA::cols; ++j)
for (int k=0; k<n; k++)
X[i][j] += A[k][i] * Y.scalar() * B[k][j];
}
template <typename FieldType, int n, class TransformationMatrix>
static void addTransformedMatrix(Dune::DiagonalMatrix<FieldType,n>& X, const TransformationMatrix& A, const Dune::DiagonalMatrix<FieldType,n>& Y, const TransformationMatrix& B)
{ {
for (int i=0; i<n; ++i) Dune::MatrixVector::transformMatrix(A,T1,B,T2);
for (int k=0; k<n; k++)
X.diagonal()[i] += A[k][i] * Y.diagonal()[k] * B[k][i];
} }
template <typename FieldType, int n, class TransformationMatrix> template<class MatrixBlockA, class MatrixB, class TransformationMatrix>
static void addTransformedMatrix(Dune::ScaledIdentityMatrix<FieldType,n>& X, const TransformationMatrix& A, const Dune::ScaledIdentityMatrix<FieldType,n>& Y, const TransformationMatrix& B) static void transformMatrixPattern(Dune::BCRSMatrix<MatrixBlockA>& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{ {
for (int k=0; k<n; k++) Dune::MatrixVector::transformMatrixPattern(A,T1,B,T2);
X.scalar() += A[k][0] * Y.scalar() * B[k][0];
} }
// compute v^T*A*v for an edge vector e = e_i - e_j with i!=j ****************************** // compute v^T*A*v for an edge vector e = e_i - e_j with i!=j ******************************
template<class Matrix> template<class Matrix>
static typename Matrix::field_type simplexEdgeDiagonal(const Matrix& A, int i, int j) static typename Matrix::field_type simplexEdgeDiagonal(const Matrix& A, int i, int j)
...@@ -252,98 +93,8 @@ class StaticMatrix ...@@ -252,98 +93,8 @@ class StaticMatrix
} }
// transform matrix A = X^T * B * Y ********************************************************
template<class MatrixA, class MatrixB, class MatrixT>
static void transformMatrix(MatrixA& A, const MatrixT& X, const MatrixB& B, const MatrixT& Y)
{
A = 0.0;
for (int k=0; k<B.N(); ++k)
for (int l=0; l<B.M(); ++l)
for (int i=0; i<X.M(); ++i)
for (int j=0; j<Y.M(); ++j)
A[i][j] = X[k][i] * B[k][l] * Y[l][j];
}
template<class MatrixA, class MatrixB>
static void transformMatrix(MatrixA& A, const double& X, const MatrixB& B, const double& Y)
{
A = 0.0;
axpy(A, X*Y, B);
}
template<class MatrixA, class MatrixT>
static void transformMatrix(MatrixA& A, const MatrixT& X, const double& B, const MatrixT& Y)
{
A = X*B*Y;
}
template<class MatrixBlockA, class MatrixB, class MatrixT>
static void transformMatrix(Dune::BCRSMatrix<MatrixBlockA>& A, const MatrixT& X, const MatrixB& B, const MatrixT& Y)
{
transformMatrixPattern(A, X, B, Y);
A = 0.0;
for (int k=0; k<B.N(); ++k)
{
typename MatrixB::row_type::ConstIterator BklIt = B[k].begin();
typename MatrixB::row_type::ConstIterator BklEnd = B[k].end();
for (; BklIt!=BklEnd; ++BklIt)
{
int l = BklIt.index();
typename MatrixT::row_type::ConstIterator XkiIt = X[k].begin();
typename MatrixT::row_type::ConstIterator XkiEnd = X[k].end();
for (; XkiIt!=XkiEnd; ++XkiIt)
{
int i = XkiIt.index();
typename MatrixT::row_type::ConstIterator YljIt = Y[l].begin();
typename MatrixT::row_type::ConstIterator YljEnd = Y[l].end();
for (; YljIt!=YljEnd; ++YljIt)
{
int j = YljIt.index();
MatrixBlockA Aij;
transformMatrix(Aij,*XkiIt, *BklIt, *YljIt);
A[i][j] += Aij;
}
}
}
}
}
template<class MatrixBlockA, class MatrixB, class MatrixT>
static void transformMatrixPattern(Dune::BCRSMatrix<MatrixBlockA>& A, const MatrixT& X, const MatrixB& B, const MatrixT& Y)
{
Dune::MatrixIndexSet indices(X.M(), Y.M());
for (int k=0; k<B.N(); ++k)
{
typename MatrixB::row_type::ConstIterator BklIt = B[k].begin();
typename MatrixB::row_type::ConstIterator BklEnd = B[k].end();
for (; BklIt!=BklEnd; ++BklIt)
{
int l = BklIt.index();
typename MatrixT::row_type::ConstIterator XkiIt = X[k].begin();
typename MatrixT::row_type::ConstIterator XkiEnd = X[k].end();
for (; XkiIt!=XkiEnd; ++XkiIt)
{
int i = XkiIt.index();
typename MatrixT::row_type::ConstIterator YljIt = Y[l].begin();
typename MatrixT::row_type::ConstIterator YljEnd = Y[l].end();
for (; YljIt!=YljEnd; ++YljIt)
{
int j = YljIt.index();
indices.add(i, j);
}
}
}
}
indices.exportIdx(A);
}
// compute i-th row of A*v for an edge vector e = e_i - e_j with i!=j **********************
template <class K> template <class K>
static Dune::FieldMatrix<K,1,1>& toMatrix(K& x) static Dune::FieldMatrix<K,1,1>& toMatrix(K& x)
{ {
...@@ -380,8 +131,58 @@ class StaticMatrix ...@@ -380,8 +131,58 @@ class StaticMatrix
return *(reinterpret_cast< Dune::FieldMatrix<K,1,1>* > (&x)); return *(reinterpret_cast< Dune::FieldMatrix<K,1,1>* > (&x));
} }
/** \brief Compute an LDL^T decomposition
}; *
* The methods computes a decomposition A=LDL^T of a given dense
* symmetric matrix A such that L is lower triangular with all
* diagonal elements equal to 1 and D is diagonal. If A is positive
* definite then A=(LD^0.5)(LD^0.5)^T is the Cholesky decomposition.
* However, the LDL^T decomposition does also work for indefinite
* symmetric matrices and is more stable than the Cholesky decomposition
* since no square roots are required.
*
* The method does only change the nontrivial entries of the given matrix
* L and D, i.e., it does not set the trivial 0 and 1 entries.
* Thus one can store both in a single matrix M and use
* M as argument for L and D.
*
* The method can furthermore work in-place, i.e., it is safe to
* use A as argument for L and D. In this case the entries of A
* below and on the diagonal are changed to those to those of
* L and D, respectively.
*
* \param A Matrix to be decomposed. Only the lower triangle is used.
* \param L Matrix to store the lower triangle. Only entries below the diagonal are written.
* \param D Matrix to store the diagonal. Only diagonal entries are written.
*/
template<class SymmetricMatrix, class LowerTriangularMatrix, class DiagonalMatrix>
static void ldlt(const SymmetricMatrix& A, LowerTriangularMatrix& L, DiagonalMatrix& D)
{
Dune::MatrixVector::ldlt(A,L,D);
}
/** \brief Solve linear system using a LDL^T decomposition.
*
* The methods solves a linear system Mx=b where A is given
* by a decomposition A=LDL^T. The method does only use
* the values of L and D below and on the diagonal, respectively.
* The 1 entries on the diagonal of L are not required.
* If L and D are stored in a single matrix it is safe
* the use this matrix as argument for both.
*
* Note that the solution vector must already have the correct size.
*
* \param L Matrix containing the lower triangle of the decomposition
* \param D Matrix containing the diagonal of the decomposition
* \param b Right hand side on the linear system
* \param x Vector to store the solution of the linear system
*/
template<class LowerTriangularMatrix, class DiagonalMatrix, class RhsVector, class SolVector>
static void ldltSolve(const LowerTriangularMatrix& L, const DiagonalMatrix& D, const RhsVector& b, SolVector& x)
{
Dune::MatrixVector::ldltSolve(L,D,b,x);
}
}
#endif #endif
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_SOLVERS_COMMON_TUPLEVECTOR_HH
#define DUNE_SOLVERS_COMMON_TUPLEVECTOR_HH
#warning This file is deprecated! Use tuplevector.hh from the dune-common module instead!
#include <dune/common/tuplevector.hh>
namespace Dune
{
namespace Solvers
{
/**
* \brief A std::tuple that allows access to its element via operator[]
*/
template<class... T>
using TupleVector = Dune::TupleVector<T...>;
} // namespace Functions
} // namespace Dune
#endif
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set ts=4 sw=2 et sts=2:
#ifndef DUNE_SOLVERS_COMMON_WRAP_OWN_SHARE_HH
#define DUNE_SOLVERS_COMMON_WRAP_OWN_SHARE_HH
#include <memory>
#include <algorithm>
#include <type_traits>
namespace Dune {
namespace Solvers {
namespace Impl {
//! Wrap const reference as shared_ptr
template<class T, class S>
std::shared_ptr<const T> wrap_own_share(const S& t)
{
return std::shared_ptr<T>(&t, [](auto*){} );
}
//! Wrap reference as shared_ptr
template<class T, class S>
std::shared_ptr<T> wrap_own_share(S& t)
{
return std::shared_ptr<T>(&t, [](auto*){} );
}
//! Move r-value reference to shared_ptr
template<class T, class S>
std::shared_ptr<T> wrap_own_share(S&& t)
{
return std::make_shared<S>(std::move(t));
}
template<class T>
std::shared_ptr<T> wrap_own_share(std::nullptr_t)
{
return {};
}
//! Share ownership of shared_ptr
template<class T, class S>
std::shared_ptr<T> wrap_own_share(std::shared_ptr<S> t)
{
return t;
}
} // end namespace Impl
/** \brief Convert l-value and r-value references, and shared_ptr of an object into a shared_ptr of a convertible type.
*
* L-value references are wrapped, r-value references are moved and shared_ptr are shared.
*
* \tparam T The target type
* \tparam S The input type
* \param s An l-value or r-value reference, or shared_ptr of type S
* \returns Shared pointer of the base class type.
*/
template <class T, class S, class Enable = std::enable_if_t<not std::is_pointer<S>::value> >
auto wrap_own_share(S&& s) {
return Impl::wrap_own_share<T>(std::forward<S>(s));
}
} // end namespace Solvers
} // namespace Dune
#endif
#ifndef CONTACT_COMPUTE_ENERGY_HH // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
#define CONTACT_COMPUTE_ENERGY_HH // vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_SOLVERS_COMPUTEENERGY_HH
#define DUNE_SOLVERS_COMPUTEENERGY_HH
#include "common/arithmetic.hh" #include <dune/matrix-vector/axy.hh>
template <class M, class V> template <class M, class V>
double computeEnergy(const M& mat, const V& x, const V& rhs) double computeEnergy(const M& mat, const V& x)
{ {
return 0.5*Arithmetic::Axy(mat,x,x) - rhs*x; return 0.5*Dune::MatrixVector::Axy(mat,x,x);
} }
template <class M, class V> template <class M, class V>
double computeEnergy(const M& mat, const V& x) double computeEnergy(const M& mat, const V& x, const V& rhs)
{ {
return 0.5*Arithmetic::Axy(mat,x,x); return computeEnergy(mat,x) - rhs*x;
} }
#endif #endif
install(FILES install(FILES
amgstep.hh amgstep.hh
blockgsstep.cc blockgssteps.hh
blockgsstep.hh
cgstep.cc cgstep.cc
cgstep.hh cgstep.hh
istlseqilu0step.hh istlseqilu0step.hh
istlseqssorstep.hh
iterationstep.hh iterationstep.hh
lineariterationstep.hh lineariterationstep.hh
linegsstep.cc linegsstep.cc
...@@ -17,6 +17,8 @@ install(FILES ...@@ -17,6 +17,8 @@ install(FILES
obstacletnnmgstep.hh obstacletnnmgstep.hh
projectedblockgsstep.cc projectedblockgsstep.cc
projectedblockgsstep.hh projectedblockgsstep.hh
projectedgradientstep.hh
projectedgradientstep.cc
projectedlinegsstep.cc projectedlinegsstep.cc
projectedlinegsstep.hh projectedlinegsstep.hh
richardsonstep.hh richardsonstep.hh
......