Commit 85d9bec2 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

merge of block-vector changes

parents 6ea82cf9 7a739148
......@@ -65,8 +65,37 @@ for (auto const& e : edges(gv)) {
## Installation
We provide a *cmake*-based configuration. See also the `INSTALL.md` file in the
source directory.
We provide a *cmake*-based configuration and use the `dunecontrol` build system.
Simply run
```
dunecontrol --current all
cmake --build build-cmake --target examples
```
to build all the example problems. The `dunecontrol` script searches for the required
(and suggested) dune modules this library depends on. These include:
- [dune-common](https://gitlab.dune-project.org/core/dune-common)
- [dune-geometry](https://gitlab.dune-project.org/core/dune-geometry)
- [dune-grid](https://gitlab.dune-project.org/core/dune-grid)
- [dune-foamgrid](https://gitlab.dune-project.org/extensions/dune-foamgrid.git)
- [dune-localfunctions](https://gitlab.dune-project.org/core/dune-localfunctions)
- [dune-typetree](https://gitlab.dune-project.org/staging/dune-typetree)
- [dune-functions](https://gitlab.dune-project.org/staging/dune-functions)
(See the file `dune.module` for an up-to-date list of dependencies). The dune modules
can be obtained from https://gitlab.dune-project.org and need to be found in a
subdirectory of `DUNE_CONTROL_PATH`. See also https://dune-project.org/doc/installation
for details about the installation of dune modules.
Additionally we require the following libraries to be found:
- [Eigen3](http://eigen.tuxfamily.org)
- [SuiteSparse](http://faculty.cse.tamu.edu/davis/suitesparse.html)
And a compiler that supports the C++14 standard, e.g. g++ >= 4.9 and clang >= 3.6.
See also the Dokerfiles in [iwr/docker-images/dune-latest](https://gitlab.math.tu-dresden.de/iwr/docker-images)
for a [docker](https://www.docker.com) container based installation.
## Documentation
......@@ -74,5 +103,5 @@ Currently only a *doxygen*-based documentation of the source files is available.
Generate a html version by
```
doxygen doc/Doxyfile
cmake --build build-cmake --target documentation
```
......@@ -19,6 +19,10 @@ namespace Dec
template <class, class T>
struct InvokeType { using type = T; };
// Workaround for MSVC (problems with alias templates in pack expansion)
template <class... Args>
struct MakeVoid { using type = void; };
#ifdef HAVE_CONCEPTS
template <class T> struct ValueType;
template <class T>
......@@ -46,6 +50,31 @@ namespace Dec
{
using type = T;
};
template <class T, class = void>
struct FieldType
{
using type = T;
};
template <class T>
struct FieldType<T, typename MakeVoid<typename T::field_type>::type>
{
using type = typename T::field_type;
};
template <class T>
struct FieldType<T, typename MakeVoid<typename T::value_type>::type>
{
using type = typename T::value_type;
};
template <class T>
struct FieldType<T*>
{
using type = T;
};
}
template <class T>
......@@ -57,6 +86,12 @@ namespace Dec
template <class T>
using Value_t = typename aux::ValueType<Decay_t<T>>::type;
template <class T>
using FieldType_t = typename aux::FieldType<Decay_t<T>>::type;
template <class... Args>
using Void_t = typename aux::MakeVoid<Args...>::type;
template <class T>
using Store_t = typename aux::StoreType<Decay_t<T>>::type;
......
......@@ -2,6 +2,8 @@
#include <type_traits>
#include "Common.hpp"
#if defined(DOXYGEN)
#define REQUIRES(...)
#define CONCEPT constexpr
......@@ -12,17 +14,6 @@
namespace Dec
{
namespace aux
{
// Workaround for MSVC (problems with alias templates in pack expansion)
template <class... Args>
struct MakeVoid { using type = void; };
} // end namespace aux
template <class... Args>
using Void_t = typename aux::MakeVoid<Args...>::type;
struct valid {};
constexpr valid valid_ = {};
......
......@@ -47,9 +47,9 @@ namespace Dec
/// \brief Constructor that accepps nested brace-init-lists, e.g. { {...}, {...} }.
/**
* Example:
* \code{.cpp}
* ```
* FixMat<double, 2, 2> A{ {1.0, 1.0}, {2.0, 1.0} };
* \endcode
* ```
**/
template <class S,
REQUIRES(concepts::Convertible<S,T>) >
......@@ -72,9 +72,9 @@ namespace Dec
/// \brief Constructor that accepts brace-init-lists, e.g. {...}.
/**
* Example:
* \code{.cpp}
* ```
* FixVec<double, 2> A{1.0, 1.0};
* \endcode
* ```
**/
template <class S,
REQUIRES(concepts::Convertible<S,T>) >
......@@ -179,14 +179,14 @@ namespace Dec
/// \brief Generates a sub-matrix of indices [I0,I1)x[J0,J1)
/**
* Example:
* \code{.cpp}
* ```
* FixMat<double, 3, 3> A{ {1.0, 1.0, 1.0}, {2.0, 1.0, 0.0}, {3.0, 1.0, -1.0} };
* FixMat<double, 2, 2> B{ {1.0, 1.0}, {2.0, 1.0} };
* FixMat<double, 1, 2> C{ {1.0, 1.0} };
*
* assert( A(range_<0,2>, range_<0,2>) == B );
* assert( A(index_<0>, range_<0,2>) == C );
* \endcode
* ```
**/
template <std::size_t I0, std::size_t I1, std::size_t J0, std::size_t J1>
view::SubMat<Self, T, I0, I1, J0, J1> operator()(range_t<I0,I1>, range_t<J0,J1>)
......@@ -209,13 +209,13 @@ namespace Dec
/// \brief Generates a sub-vector of indices [I0,I1)
/**
* Example:
* \code{.cpp}
* ```
* FixVec<double, 2> a{2.0, 1.0};
* FixMat<double, 1> b{1.0};
*
* assert( a[index_<1>] == b );
* assert( a[range_<1,2>] == b );
* \endcode
* ```
**/
template <std::size_t I0, std::size_t I1>
view::SubVec<Self, T, N, M, I0, I1> operator[](range_t<I0,I1>)
......@@ -242,6 +242,11 @@ namespace Dec
/// @endcond
/** @} **/
using Super::begin;
using Super::end;
using Super::cbegin;
using Super::cend;
using Expr::num_rows;
using Expr::num_cols;
......
......@@ -24,6 +24,14 @@ namespace Dec
result += conj_op(mat1(i,j)) * mat2(i,j);
return result;
}
template <class T, class S,
REQUIRES(concepts::Arithmetic<T> && concepts::Arithmetic<S>) >
auto dot(T const& x, S const& y)
{
operation::hermitian<T> conj_op{};
return conj_op(x) * y;
}
/** \relates FixMat
* \brief Dot-product of two vectors = lhs^H * rhs
......@@ -32,14 +40,11 @@ namespace Dec
REQUIRES(concepts::Vectors<V0,V1>) >
auto dot(V0 const& lhs, V1 const& rhs)
{
using T0 = Value_t<V0>;
using T1 = Value_t<V1>;
using T = Decay_t<decltype( std::declval<T0>() * std::declval<T1>() )>;
assert (V0::size() > 0);
operation::hermitian<T1> conj_op{};
T result = 0;
for (std::size_t i = 0; i < V0::size(); ++i)
result += conj_op(lhs[i]) * rhs[i];
auto result = dot(lhs[0], rhs[0]);
for (std::size_t i = 1; i < V0::size(); ++i)
result += dot(lhs[i], rhs[i]);
return result;
}
......@@ -47,9 +52,9 @@ namespace Dec
namespace aux
{
template <class E, class T, std::size_t N, std::size_t M, class Operation>
T accumulate(MatExpr<E, T, N, M> const& mat, Operation op)
FieldType_t<T> accumulate(MatExpr<E, T, N, M> const& mat, Operation op)
{
T result = 0;
FieldType_t<T> result = 0;
for (std::size_t i = 0; i < N; ++i)
for (std::size_t j = 0; j < M; ++j)
result = op(result, mat(i,j));
......@@ -58,6 +63,10 @@ namespace Dec
} // end namespace aux
template <class T,
REQUIRES(concepts::Arithmetic<T>) >
T unary_dot(T const& x) { return sqr(std::abs(x)); }
/** \relates FixMat
* \brief Dot-product with the vector itself
......@@ -66,7 +75,7 @@ namespace Dec
REQUIRES(concepts::Vector<V>) >
auto unary_dot(V const& x)
{
auto op = [](auto const& a, auto const& b) { return a + sqr(std::abs(b)); };
auto op = [](auto const& a, auto const& b) { return a + unary_dot(b); };
return aux::accumulate(x, op);
}
......
......@@ -13,9 +13,9 @@ namespace Dec
{
template <class T, std::size_t N, std::size_t M>
class BlockMatrix
: public std::array<std::array<DOFMatrix<T>, M>, N>
: public FixMat<DOFMatrix<T>, N, M> // std::array<std::array<DOFMatrix<T>, M>, N>
{
using Super = std::array<std::array<DOFMatrix<T>, M>, N>;
using Super = FixMat<DOFMatrix<T>, N, M>; //std::array<std::array<DOFMatrix<T>, M>, N>;
public:
......@@ -30,8 +30,8 @@ namespace Dec
}
/// Access the (i,j)'th block
DOFMatrix<T>& operator()(size_type i, size_type j) { return (*this)[i][j]; }
DOFMatrix<T> const& operator()(size_type i, size_type j) const { return (*this)[i][j]; }
// DOFMatrix<T>& operator()(size_type i, size_type j) { return (*this)[i][j]; }
// DOFMatrix<T> const& operator()(size_type i, size_type j) const { return (*this)[i][j]; }
/// Return number of non-zeros
std::size_t nonZeros() const
......
......@@ -38,9 +38,9 @@ namespace Dec
template <class T, std::size_t N>
class BlockVector
: public std::array<DOFVector<T>, N>
: public FixVec<DOFVector<T>, N> //std::array<DOFVector<T>, N>
{
using Super = std::array<DOFVector<T>, N>;
using Super = FixVec<DOFVector<T>, N>; //std::array<DOFVector<T>, N>;
public:
using value_type = T;
......@@ -111,48 +111,50 @@ namespace Dec
return *this;
}
BlockVector& operator+=(BlockVector const& that)
{
for (std::size_t i = 0; i < N; ++i)
(*this)[i] += that[i];
return *this;
}
BlockVector& operator-=(BlockVector const& that)
{
for (std::size_t i = 0; i < N; ++i)
(*this)[i] -= that[i];
return *this;
}
BlockVector& operator*=(float_type factor)
{
for (std::size_t i = 0; i < N; ++i)
(*this)[i] *= factor;
return *this;
}
template <class M>
BlockVector& operator=(LinearAlgebraExpression<M> const& expr)
{
expr.assign(*this);
return *this;
}
template <class M>
BlockVector& operator+=(LinearAlgebraExpression<M> const& expr)
{
expr.add_assign(*this);
return *this;
}
template <class M>
BlockVector& operator-=(LinearAlgebraExpression<M> const& expr)
{
expr.add_assign(*this, -1.0);
return *this;
}
using Super::operator=;
// BlockVector& operator+=(BlockVector const& that)
// {
// for (std::size_t i = 0; i < N; ++i)
// (*this)[i] += that[i];
// return *this;
// }
//
// BlockVector& operator-=(BlockVector const& that)
// {
// for (std::size_t i = 0; i < N; ++i)
// (*this)[i] -= that[i];
// return *this;
// }
//
// BlockVector& operator*=(float_type factor)
// {
// for (std::size_t i = 0; i < N; ++i)
// (*this)[i] *= factor;
// return *this;
// }
//
// template <class M>
// BlockVector& operator=(LinearAlgebraExpression<M> const& expr)
// {
// expr.assign(*this);
// return *this;
// }
//
// template <class M>
// BlockVector& operator+=(LinearAlgebraExpression<M> const& expr)
// {
// expr.add_assign(*this);
// return *this;
// }
//
// template <class M>
// BlockVector& operator-=(LinearAlgebraExpression<M> const& expr)
// {
// expr.add_assign(*this, -1.0);
// return *this;
// }
/// Sets all vectors to zero
void setZero()
......
......@@ -18,30 +18,30 @@ namespace Dec
return vec;
}
template <class T, std::size_t N>
T unary_dot(BlockVector<T,N> const& a)
{
T result{};
for (std::size_t i = 0; i < N; ++i)
result += unary_dot(a[i]);
return result;
}
template <class T, std::size_t N>
inline T two_norm(BlockVector<T,N> const& vec)
{
using std::sqrt;
return sqrt(unary_dot(vec));
}
template <class T, std::size_t N>
T dot(BlockVector<T,N> const& a, BlockVector<T,N> const& b)
{
T result{};
for (std::size_t i = 0; i < N; ++i)
result += dot(a[i], b[i]);
return result;
}
// template <class T, std::size_t N>
// T unary_dot(BlockVector<T,N> const& a)
// {
// T result{};
// for (std::size_t i = 0; i < N; ++i)
// result += unary_dot(a[i]);
// return result;
// }
//
// template <class T, std::size_t N>
// inline T two_norm(BlockVector<T,N> const& vec)
// {
// using std::sqrt;
// return sqrt(unary_dot(vec));
// }
// template <class T, std::size_t N>
// T dot(BlockVector<T,N> const& a, BlockVector<T,N> const& b)
// {
// T result{};
// for (std::size_t i = 0; i < N; ++i)
// result += dot(a[i], b[i]);
// return result;
// }
template <class T, std::size_t N, std::size_t M>
......@@ -70,90 +70,90 @@ namespace Dec
}
//< lhs + rhs
template <class T, std::size_t N>
inline VecAddExpression<T,N> operator+(BlockVector<T,N> const& lhs, BlockVector<T,N> const& rhs)
{
return {lhs, rhs.self()};
}
//< lhs - rhs
template <class T, std::size_t N>
inline VecAddExpression<T,N> operator-(BlockVector<T,N> lhs, BlockVector<T,N> const& rhs)
{
return {lhs, rhs.self(), -1};
}
//< lhs + rhs
template <class T, std::size_t N, class M>
inline BlockVector<T,N> operator+(BlockVector<T,N> lhs, LinearAlgebraExpression<M> const& rhs)
{
return lhs += rhs;
}
//< lhs - rhs
template <class T, std::size_t N, class M>
inline BlockVector<T,N> operator-(BlockVector<T,N> lhs, LinearAlgebraExpression<M> const& rhs)
{
return lhs -= rhs;
}
//< lhs + rhs
template <class T, std::size_t N, class M>
inline BlockVector<T,N> operator+(VecScaleExpression<T,N> const& lhs, LinearAlgebraExpression<M> const& rhs)
{
BlockVector<T,N> result{lhs.vec_, tag::ressource{}}; result = lhs;
return result += rhs;
}
//< lhs - rhs
template <class T, std::size_t N, class M>
inline BlockVector<T,N> operator-(VecScaleExpression<T,N> const& lhs, LinearAlgebraExpression<M> const& rhs)
{
BlockVector<T,N> result{lhs.vec_, tag::ressource{}}; result = lhs;
return result -= rhs;
}
//< lhs + rhs
template <class T, std::size_t N, class M>
inline BlockVector<T,N> operator+(VecAddExpression<T,N> const& lhs, LinearAlgebraExpression<M> const& rhs)
{
BlockVector<T,N> result{lhs.x_, tag::ressource{}}; result = lhs;
return result += rhs;
}
//< lhs - rhs
template <class T, std::size_t N, class M>
inline BlockVector<T,N> operator-(VecAddExpression<T,N> const& lhs, LinearAlgebraExpression<M> const& rhs)
{
BlockVector<T,N> result{lhs.x_, tag::ressource{}}; result = lhs;
return result -= rhs;
}
//< lhs * factor
template <class T, std::size_t N>
inline VecScaleExpression<T,N> operator*(BlockVector<T,N> const& v, T const& factor)
{
return {v, factor};
}
//< factor * lhs
template <class T, std::size_t N>
inline VecScaleExpression<T,N> operator*(T const& factor, BlockVector<T,N> const& v)
{
return {v, factor};
}
//< lhs / factor
template <class T, std::size_t N>
inline VecScaleExpression<T,N> operator/(BlockVector<T,N> const& v, T const& factor)
{
return {v, 1.0/factor};
}
// //< lhs + rhs
// template <class T, std::size_t N>
// inline VecAddExpression<T,N> operator+(BlockVector<T,N> const& lhs, BlockVector<T,N> const& rhs)
// {
// return {lhs, rhs.self()};
// }
//
// //< lhs - rhs
// template <class T, std::size_t N>
// inline VecAddExpression<T,N> operator-(BlockVector<T,N> lhs, BlockVector<T,N> const& rhs)
// {
// return {lhs, rhs.self(), -1};
// }
//
//
// //< lhs + rhs
// template <class T, std::size_t N, class M>
// inline BlockVector<T,N> operator+(BlockVector<T,N> lhs, LinearAlgebraExpression<M> const& rhs)
// {
// return lhs += rhs;
// }
//
// //< lhs - rhs
// template <class T, std::size_t N, class M>
// inline BlockVector<T,N> operator-(BlockVector<T,N> lhs, LinearAlgebraExpression<M> const& rhs)
// {
// return lhs -= rhs;
// }
//
//
// //< lhs + rhs
// template <class T, std::size_t N, class M>
// inline BlockVector<T,N> operator+(VecScaleExpression<T,N> const& lhs, LinearAlgebraExpression<M> const& rhs)
// {
// BlockVector<T,N> result{lhs.vec_, tag::ressource{}}; result = lhs;
// return result += rhs;
// }
//
// //< lhs - rhs
// template <class T, std::size_t N, class M>
// inline BlockVector<T,N> operator-(VecScaleExpression<T,N> const& lhs, LinearAlgebraExpression<M> const& rhs)
// {
// BlockVector<T,N> result{lhs.vec_, tag::ressource{}}; result = lhs;
// return result -= rhs;
// }
//
//
//
// //< lhs + rhs
// template <class T, std::size_t N, class M>
// inline BlockVector<T,N> operator+(VecAddExpression<T,N> const& lhs, LinearAlgebraExpression<M> const& rhs)
// {
// BlockVector<T,N> result{lhs.x_, tag::ressource{}}; result = lhs;
// return result += rhs;
// }
//
// //< lhs - rhs
// template <class T, std::size_t N, class M>
// inline BlockVector<T,N> operator-(VecAddExpression<T,N> const& lhs, LinearAlgebraExpression<M> const& rhs)
// {
// BlockVector<T,N> result{lhs.x_, tag::ressource{}}; result = lhs;
// return result -= rhs;
// }
//
//
// //< lhs * factor
// template <class T, std::size_t N>
// inline VecScaleExpression<T,N> operator*(BlockVector<T,N> const& v, T const& factor)
// {
// return {v, factor};
// }
//
// //< factor * lhs
// template <class T, std::size_t N>