diff --git a/.gitsuper b/.gitsuper
index a287144c915b2b3b9933f93cc49aef47bb9e40d3..d6c9f35cccdd9befdfd61eb52c5f3b6e595098e9 100644
--- a/.gitsuper
+++ b/.gitsuper
@@ -1,26 +1,26 @@
[supermodule]
remote = git@github.com:dune-community/dune-gdt-super.git
status = 1a3bcab04b011a5d6e44f9983cae6ff89fa695e8 bin (heads/master)
- 5af7224dda0addd6cde1a0073ff93fc9be1f5095 config.opts (heads/master-7-g5af7224)
+ 20a673b9dad7e2e25bd97defa8849debb59d247c config.opts (heads/master)
8f2c5aba441417bf2c42f22272f538c68a89cc4a dune-alugrid (remotes/origin/releases/2.5)
707acf201d5a754c80f87cc4d71aa36bf29a6e3f dune-common (v2.5.1-9-g707acf20)
- +a57cba2fd8f8434cd9690566d5af6090cbd5a38f dune-gdt (heads/dailywork_tleibner)
+ +fa850d2b3e957e3083f8ef49c5f1b1d097abe2aa dune-gdt (heads/entropy_flux_even_newer)
390a2c503783bbed778a8ff610f8c5ca09c238d0 dune-geometry (v2.5.1-5-g390a2c5)
d7b20bbc5f6fdcfc312beb0ea5d16d39ea26904e dune-grid (v2.5.1-2-gd7b20bbc5)
- e9d9a3336735090648637e044e279866bbea3597 dune-grid-glue (v2.4.0-60-ge9d9a33)
+ +e9d9a3336735090648637e044e279866bbea3597 dune-grid-glue (v2.4.0-60-ge9d9a33)
63df56a54f81eda308233a683eb329e77e69f0a9 dune-istl (v2.5.1rc1)
0d757d65e5d57134a7ecf304e35d063f4ccc7116 dune-localfunctions (v2.5.1rc1)
8a69fc68165780921bbba77da338b6932daf983c dune-pybindxi (v2.2.1-16-g8a69fc6)
741e4f8e53bdd3e1b6e19d84eb22b6e3dc48526c dune-python (remotes/origin/releases/2.5)
26cc8cb4161a3a51002ab2a81b8c81d2c951ee79 dune-testtools (remotes/origin/p/renemilk/testname_listing_hack_no-skiptest)
- 8fe883e99c58c9f0c2f92457d546a0ac9f5a9bf9 dune-uggrid (v2.5.2-1-g8fe883e9)
- +0caf50617b246fca4d8e350ddd9d16dbd4f7216f dune-xt-common (heads/dailywork_tleibner)
- +b8cf8efc32e46511990b691afee55da3b7af15c8 dune-xt-data (heads/master)
- +d02b6e380acce0b86bac7cb1de0e7e9dbb7401c1 dune-xt-functions (heads/master)
- +dd30fcd7d4485eb2a8158d5ddf01333f58502c40 dune-xt-grid (heads/master)
- +8bd7ad451015182fc009914da3527846e78d08ca dune-xt-la (heads/remove_cow)
- 60cd896bf3f1eb99066563aa1a07113e5a791d47 scripts (heads/master)
-commit = c15b7428201be08e84b3fb9e766bbdf51d9ea2b4
+ 0a74e7dd0b2115778a5d490dab08a2ed07fcaa1e dune-uggrid (v2.5.2)
+ 7d248b5ae4b59e5a9b73034fe9943e22a10cfde0 dune-xt-common (heads/dailywork_tleibner)
+ b8cf8efc32e46511990b691afee55da3b7af15c8 dune-xt-data (heads/master)
+ 0c9df39934e23b950f357c912b4c800b86de683f dune-xt-functions (heads/dailywork_tleibner)
+ dd30fcd7d4485eb2a8158d5ddf01333f58502c40 dune-xt-grid (heads/master)
+ fab7ebec48ebafb5cf541ea3ff696955b480a768 dune-xt-la (heads/master)
+ 09d0378f616b94d68bcdd9fc6114813181849ec0 scripts (heads/master)
+commit = bdd5c2f0a34cf97a088447ea76c3866dac9a3576
[submodule.bin]
remote = git@github.com:dune-community/local-bin.git
@@ -30,7 +30,7 @@ commit = 1a3bcab04b011a5d6e44f9983cae6ff89fa695e8
[submodule.config.opts]
remote = git@github.com:dune-community/config.opts.git
status =
-commit = 5af7224dda0addd6cde1a0073ff93fc9be1f5095
+commit = 20a673b9dad7e2e25bd97defa8849debb59d247c
[submodule.dune-alugrid]
remote = https://github.com/dune-mirrors/dune-alugrid.git
@@ -45,7 +45,7 @@ commit = 707acf201d5a754c80f87cc4d71aa36bf29a6e3f
[submodule.dune-gdt]
remote = git@github.com:dune-community/dune-gdt.git
status = 2424627f0ad5de7e4aaa5e7f48bc2a02414d95a1 .vcsetup (heads/master)
-commit = a57cba2fd8f8434cd9690566d5af6090cbd5a38f
+commit = fa850d2b3e957e3083f8ef49c5f1b1d097abe2aa
[submodule.dune-geometry]
remote = git@github.com:dune-community/dune-geometry.git
@@ -90,12 +90,12 @@ commit = 26cc8cb4161a3a51002ab2a81b8c81d2c951ee79
[submodule.dune-uggrid]
remote = https://github.com/dune-mirrors/dune-uggrid.git
status =
-commit = 8fe883e99c58c9f0c2f92457d546a0ac9f5a9bf9
+commit = 0a74e7dd0b2115778a5d490dab08a2ed07fcaa1e
[submodule.dune-xt-common]
remote = git@github.com:dune-community/dune-xt-common.git
status = 2424627f0ad5de7e4aaa5e7f48bc2a02414d95a1 .vcsetup (heads/master)
-commit = 0caf50617b246fca4d8e350ddd9d16dbd4f7216f
+commit = 7d248b5ae4b59e5a9b73034fe9943e22a10cfde0
[submodule.dune-xt-data]
remote = https://github.com/dune-community/dune-xt-data
@@ -104,8 +104,8 @@ commit = b8cf8efc32e46511990b691afee55da3b7af15c8
[submodule.dune-xt-functions]
remote = git@github.com:dune-community/dune-xt-functions.git
-status = 2424627f0ad5de7e4aaa5e7f48bc2a02414d95a1 .vcsetup (heads/master)
-commit = d02b6e380acce0b86bac7cb1de0e7e9dbb7401c1
+status = 2424627f0ad5de7e4aaa5e7f48bc2a02414d95a1 .vcsetup ((null))
+commit = 0c9df39934e23b950f357c912b4c800b86de683f
[submodule.dune-xt-grid]
remote = git@github.com:dune-community/dune-xt-grid.git
@@ -115,10 +115,10 @@ commit = dd30fcd7d4485eb2a8158d5ddf01333f58502c40
[submodule.dune-xt-la]
remote = git@github.com:dune-community/dune-xt-la.git
status = 2424627f0ad5de7e4aaa5e7f48bc2a02414d95a1 .vcsetup (heads/master)
-commit = 8bd7ad451015182fc009914da3527846e78d08ca
+commit = fab7ebec48ebafb5cf541ea3ff696955b480a768
[submodule.scripts]
remote = https://github.com/wwu-numerik/scripts.git
status = fb5ebc10e647d637c69497af2ec2560847eb2112 python/pylicense (v0.2.0~10)
-commit = 60cd896bf3f1eb99066563aa1a07113e5a791d47
+commit = 09d0378f616b94d68bcdd9fc6114813181849ec0
diff --git a/dune/xt/common/fmatrix.hh b/dune/xt/common/fmatrix.hh
index abed7171b848b2117c31fb0439fe72260cbe82d5..3686656a7bdd732928de4830faac6db1cadf5448 100644
--- a/dune/xt/common/fmatrix.hh
+++ b/dune/xt/common/fmatrix.hh
@@ -201,7 +201,8 @@ inline void FieldMatrix<K, ROWS, COLS>::luDecomposition(FieldMatrix<K, ROWS, COL
if (imax != i) {
for (size_type j = 0; j < ROWS; j++)
std::swap(A[i][j], A[imax][j]);
- func.swap(i, imax); // swap the pivot or rhs
+ assert(imax < std::numeric_limits<int>::max() && i < std::numeric_limits<int>::max());
+ func.swap(static_cast<int>(i), static_cast<int>(imax)); // swap the pivot or rhs
}
// singular ?
@@ -214,7 +215,8 @@ inline void FieldMatrix<K, ROWS, COLS>::luDecomposition(FieldMatrix<K, ROWS, COL
A[k][i] = factor;
for (size_type j = i + 1; j < ROWS; j++)
A[k][j] -= factor * A[i][j];
- func(factor, k, i);
+ assert(k < std::numeric_limits<int>::max() && i < std::numeric_limits<int>::max());
+ func(factor, static_cast<int>(k), static_cast<int>(i));
}
}
}
@@ -370,6 +372,143 @@ public:
}; // class FieldMatrix
+template <class K, size_t num_blocks, size_t block_rows, size_t block_cols = block_rows>
+class BlockedFieldMatrix
+{
+ using ThisType = BlockedFieldMatrix;
+
+public:
+ static constexpr size_t num_rows = num_blocks * block_rows;
+ static constexpr size_t num_cols = num_blocks * block_cols;
+ using MatrixType = Dune::FieldMatrix<K, num_rows, num_cols>;
+ using BlockType = FieldMatrix<K, block_rows, block_cols>;
+
+ BlockedFieldMatrix(const K& val = K(0.))
+ : backend_(BlockType(val))
+ {
+ }
+
+ BlockedFieldMatrix(const MatrixType& other)
+ {
+ *this = other;
+ }
+
+ BlockedFieldMatrix(const BlockType& block)
+ : backend_(block)
+ {
+ }
+
+ ThisType& operator=(const MatrixType& other)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t ll = 0; ll < block_rows; ++ll)
+ for (size_t mm = 0; mm < block_cols; ++mm)
+ backend_[jj][ll][mm] = other[jj * block_rows + ll][jj * block_cols + mm];
+ return *this;
+ }
+
+ bool operator==(const ThisType& other) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ if (block(jj) != other.block(jj))
+ return false;
+ return true;
+ }
+
+ K get_entry(const size_t ii, const size_t jj) const
+ {
+ assert(ii < num_rows && jj < num_cols);
+ if (ii / block_rows != jj / block_cols)
+ return K(0.);
+ return backend_[ii / block_rows][ii % block_rows][jj % block_cols];
+ }
+
+ K& get_entry(const size_t jj, const size_t ll, const size_t mm)
+ {
+ assert(jj < num_blocks && ll < block_rows && mm < block_cols);
+ return backend_[jj][ll][mm];
+ }
+
+ const K& get_entry(const size_t jj, const size_t ll, const size_t mm) const
+ {
+ assert(jj < num_blocks && ll < block_rows && mm < block_cols);
+ return backend_[jj][ll][mm];
+ }
+
+ BlockType& block(const size_t jj)
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ const BlockType& block(const size_t jj) const
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ void mv(const Dune::FieldVector<K, num_cols>& x, Dune::FieldVector<K, num_rows>& ret) const
+ {
+ std::fill(ret.begin(), ret.end(), 0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto row_offset = block_rows * jj;
+ const auto col_offset = block_cols * jj;
+ for (size_t ll = 0; ll < block_rows; ++ll)
+ for (size_t mm = 0; mm < block_cols; ++mm)
+ ret[row_offset + ll] += backend_[jj][ll][mm] * x[col_offset + mm];
+ } // jj
+ } // void mv(...)
+
+ void mv(const BlockedFieldVector<K, num_blocks, block_cols>& x,
+ BlockedFieldVector<K, num_blocks, block_rows>& ret) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj].mv(x.block(jj), ret.block(jj));
+ } // void mv(...)
+
+ void mtv(const Dune::FieldVector<K, num_rows>& x, Dune::FieldVector<K, num_cols>& ret) const
+ {
+ std::fill(ret.begin(), ret.end(), 0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto row_offset = block_rows * jj;
+ const auto col_offset = block_cols * jj;
+ for (size_t mm = 0; mm < block_cols; ++mm)
+ for (size_t ll = 0; ll < block_rows; ++ll)
+ ret[col_offset + mm] += backend_[jj][ll][mm] * x[row_offset + ll];
+ } // jj
+ } // void mtv(...)
+
+ void mtv(const BlockedFieldVector<K, num_blocks, block_rows>& x,
+ BlockedFieldVector<K, num_blocks, block_cols>& ret) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj].mtv(x.block(jj), ret.block(jj));
+ } // void mv(...)
+
+ template <size_t br, size_t bc>
+ ThisType& rightmultiply(const BlockedFieldMatrix<K, num_blocks, br, bc>& other)
+ {
+ assert((this != &other) && "Multiplying a matrix by itself gives wrong results, please copy before!");
+ static_assert(br == bc, "Cannot rightmultiply with non-square matrix");
+ static_assert(br == block_cols, "Size mismatch");
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj].rightmultiply(other.backend_[jj]);
+ return *this;
+ }
+
+ BlockedFieldMatrix<K, num_blocks, block_cols, block_rows> transpose()
+ {
+ BlockedFieldMatrix<K, num_blocks, block_cols, block_rows> ret;
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ ret.block(jj) = block(jj).transpose();
+ return ret;
+ }
+
+private:
+ FieldVector<BlockType, num_blocks> backend_;
+};
+
+
template <class K, int N, int M>
struct MatrixAbstraction<Dune::XT::Common::FieldMatrix<K, N, M>>
{
diff --git a/dune/xt/common/fvector.hh b/dune/xt/common/fvector.hh
index 3aa946ad34892b948d8221ecdcf94702bd3e6764..b4e47092f778d3ae840b45635b89d1de21c5e47f 100644
--- a/dune/xt/common/fvector.hh
+++ b/dune/xt/common/fvector.hh
@@ -13,10 +13,11 @@
#ifndef DUNE_XT_COMMON_FVECTOR_HH
#define DUNE_XT_COMMON_FVECTOR_HH
+#include <functional>
#include <initializer_list>
+#include <numeric>
#include <type_traits>
#include <vector>
-#include <functional>
#include <boost/functional/hash.hpp>
@@ -116,7 +117,7 @@ public:
return ret;
}
- template <int S>
+ template <int S = SIZE>
operator typename std::enable_if<(S == SIZE) && (SIZE != 1), Dune::FieldMatrix<K, S, 1>>::type() const
{
Dune::FieldMatrix<K, SIZE, 1> ret;
@@ -168,6 +169,196 @@ public:
}; // class FieldVector
+template <class K, size_t block_num, size_t size_block>
+class BlockedFieldVector
+{
+ using ThisType = BlockedFieldVector;
+
+public:
+ static constexpr size_t num_blocks = block_num;
+ static constexpr size_t block_size = size_block;
+ static constexpr size_t static_size = num_blocks * block_size;
+ using VectorType = Dune::FieldVector<K, static_size>;
+ using BlockType = FieldVector<K, block_size>;
+
+ BlockedFieldVector(const K& val = K(0.))
+ : backend_(BlockType(val))
+ {
+ }
+
+ BlockedFieldVector(const VectorType& other)
+ {
+ *this = other;
+ }
+
+ BlockedFieldVector(const BlockType& block)
+ : backend_(block)
+ {
+ }
+
+ ThisType& operator=(const VectorType& other)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t ii = 0; ii < block_size; ++ii)
+ backend_[jj][ii] = other[jj * block_size + ii];
+ return *this;
+ }
+
+ size_t size() const
+ {
+ return static_size;
+ }
+
+ K& operator[](const size_t ii)
+ {
+ assert(ii < static_size);
+ return backend_[ii / block_size][ii % block_size];
+ }
+
+ const K& operator[](const size_t ii) const
+ {
+ assert(ii < static_size);
+ return backend_[ii / block_size][ii % block_size];
+ }
+
+ K& get_entry(const size_t jj, const size_t ii)
+ {
+ assert(jj < num_blocks);
+ assert(ii < block_size);
+ return backend_[jj][ii];
+ }
+
+ const K& get_entry(const size_t jj, const size_t ii) const
+ {
+ assert(jj < num_blocks);
+ assert(ii < block_size);
+ return backend_[jj][ii];
+ }
+
+ BlockType& block(const size_t jj)
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ const BlockType& block(const size_t jj) const
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ K one_norm() const
+ {
+ K ret(0);
+ for (const auto& block_vec : backend_)
+ for (const auto& entry : block_vec)
+ ret += std::abs(entry);
+ return ret;
+ }
+
+ K two_norm() const
+ {
+ return std::sqrt(two_norm2());
+ }
+
+ K two_norm2() const
+ {
+ K ret(0);
+ for (const auto& block_vec : backend_)
+ for (const auto& entry : block_vec)
+ ret += std::pow(entry, 2);
+ return ret;
+ }
+
+ operator VectorType() const
+ {
+ VectorType ret(0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t ii = 0; ii < block_size; ++ii)
+ ret[jj * block_size + ii] = backend_[jj][ii];
+ return ret;
+ }
+
+ K* data()
+ {
+ return &(backend_[0][0]);
+ }
+
+ const K* data() const
+ {
+ return &(backend_[0][0]);
+ }
+
+ K* begin()
+ {
+ return data();
+ }
+
+ const K* begin() const
+ {
+ return data();
+ }
+
+ K* end()
+ {
+ return &(backend_[num_blocks - 1][block_size]);
+ }
+
+ const K* end() const
+ {
+ return &(backend_[num_blocks - 1][block_size]);
+ }
+
+ ThisType& operator*=(const K& val)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj] *= val;
+ return *this;
+ }
+
+ ThisType operator*(const K& val) const
+ {
+ auto ret = *this;
+ ret *= val;
+ return ret;
+ }
+
+ K operator*(const ThisType& other) const
+ {
+ return std::inner_product(begin(), end(), other.begin(), 0.);
+ }
+
+ ThisType& operator+=(const ThisType& other)
+ {
+ backend_ += other.backend_;
+ return *this;
+ }
+
+ ThisType operator+(const ThisType& other) const
+ {
+ auto ret = *this;
+ ret += other;
+ return ret;
+ }
+
+ ThisType& operator-=(const ThisType& other)
+ {
+ backend_ -= other.backend_;
+ return *this;
+ }
+
+ ThisType operator-(const ThisType& other) const
+ {
+ auto ret = *this;
+ ret -= other;
+ return ret;
+ }
+
+private:
+ FieldVector<BlockType, num_blocks> backend_;
+};
+
+
//! this allows to set the init value of the FieldVector at compile time
template <class K, int SIZE, K value>
class ValueInitFieldVector : public Dune::XT::Common::FieldVector<K, SIZE>
@@ -228,6 +419,34 @@ struct VectorAbstraction<Dune::XT::Common::FieldVector<K, SIZE>>
};
+//! Specialization of VectorAbstraction for Dune::XT::Common::BlockedFieldVector
+template <class K, size_t num_blocks, size_t block_size>
+struct VectorAbstraction<Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>>
+ : public internal::VectorAbstractionBase<Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>, K>,
+ public internal::
+ HasSubscriptOperatorForVectorAbstraction<Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>,
+ typename Dune::FieldTraits<K>::field_type>
+{
+ using VectorType = Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>;
+ static constexpr bool has_static_size = true;
+ static constexpr size_t static_size = VectorType::static_size;
+ static constexpr bool is_contiguous = true;
+
+ template <size_t SZ = static_size>
+ static inline VectorType create(const size_t sz, const K& val = suitable_default<K>::value())
+ {
+ static_assert(SZ == static_size, "Creation of Vector with different size not implemented!");
+ if (sz != SZ)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldVector< ..., " << SZ << " > (of "
+ << "static size) with "
+ << sz
+ << " elements!");
+ return VectorType(val);
+ }
+};
+
+
template <class V>
typename std::
enable_if<is_vector<V>::value && VectorAbstraction<V>::has_static_size,
diff --git a/dune/xt/common/test/fmatrix.cc b/dune/xt/common/test/fmatrix.cc
new file mode 100644
index 0000000000000000000000000000000000000000..924276bcdb33b26a61798b85dbef31aa57659986
--- /dev/null
+++ b/dune/xt/common/test/fmatrix.cc
@@ -0,0 +1,283 @@
+// This file is part of the dune-xt-common project:
+// https://github.com/dune-community/dune-xt-common
+// Copyright 2009-2018 dune-xt-common developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2017)
+// Rene Milk (2018)
+
+#include <dune/xt/common/test/main.hxx> // <- Needs to come first, include the config.h.
+
+#include <dune/xt/common/fmatrix.hh>
+#include <dune/xt/common/float_cmp.hh>
+
+GTEST_TEST(dune_xt_common_field_matrix, creation_and_calculations)
+{
+ static constexpr size_t rows = 2;
+ static constexpr size_t cols = 3;
+ using MatrixType = Dune::XT::Common::FieldMatrix<double, rows, cols>;
+ using SquareMatrixType = Dune::XT::Common::FieldMatrix<double, rows, rows>;
+ using BaseMatrixType = Dune::FieldMatrix<double, rows, cols>;
+ using ScalarMatrixType = Dune::XT::Common::FieldMatrix<double, 1, 1>;
+ using ScalarBaseMatrixType = Dune::FieldMatrix<double, 1, 1>;
+ using RowVectorType = Dune::XT::Common::FieldVector<double, cols>;
+ using ColVectorType = Dune::XT::Common::FieldVector<double, rows>;
+ MatrixType zeros;
+ ScalarMatrixType zero;
+ MatrixType ones(1.);
+ MatrixType ones2(rows, cols, 1.);
+ ScalarMatrixType one(1.);
+ ScalarMatrixType one2{1.};
+ MatrixType mat({{0., 1., 2.}, {3., 4., 5.}});
+ EXPECT_ANY_THROW(MatrixType({{0., 2.}, {3., 4., 5.}}));
+ EXPECT_ANY_THROW(MatrixType({{0., 1., 2.}, {4., 5.}}));
+ EXPECT_ANY_THROW(MatrixType({{0., 2.}, {4., 5.}}));
+ EXPECT_ANY_THROW(ScalarMatrixType({{0., 2.}, {4., 5.}}));
+#ifdef NDEBUG
+ EXPECT_ANY_THROW(MatrixType(rows + 1, cols - 1, 1.));
+ EXPECT_ANY_THROW(ScalarMatrixType(2, 3, 1.));
+#endif
+ EXPECT_EQ(zero[0][0], 0.);
+ EXPECT_EQ(one[0][0], 1.);
+ EXPECT_EQ(one2[0][0], 1.);
+ EXPECT_EQ(one, 1.);
+ EXPECT_EQ(ones, ones2);
+ for (size_t ii = 0; ii < rows; ++ii) {
+ for (size_t jj = 0; jj < cols; ++jj) {
+ EXPECT_DOUBLE_EQ(mat[ii][jj], ii * cols + jj);
+ EXPECT_EQ(zeros[ii][jj], 0.);
+ EXPECT_EQ(ones[ii][jj], 1.);
+ }
+ }
+ MatrixType copy_assigned = mat;
+ MatrixType copy_constructed(mat);
+ ScalarMatrixType scalar_copy_assigned = one;
+ ScalarMatrixType scalar_copy_constructed(one);
+ BaseMatrixType base_mat = mat;
+ MatrixType base_constructed(base_mat);
+ ScalarBaseMatrixType scalar_base_mat = one;
+ ScalarMatrixType scalar_base_constructed(scalar_base_mat);
+ EXPECT_EQ(copy_assigned, mat);
+ EXPECT_EQ(copy_constructed, mat);
+ EXPECT_EQ(base_constructed, mat);
+ EXPECT_EQ(scalar_copy_assigned, one);
+ EXPECT_EQ(scalar_copy_constructed, one);
+ EXPECT_EQ(scalar_base_constructed, one);
+ auto transposed_mat = mat.transpose();
+ auto scalar_transposed_mat = one.transpose();
+ EXPECT_EQ(scalar_transposed_mat, one);
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj = 0; jj < cols; ++jj)
+ EXPECT_EQ(transposed_mat[jj][ii], mat[ii][jj]);
+ RowVectorType ones_vec(1.);
+ ColVectorType mat_times_vec2;
+ auto mat_times_vec = mat * ones_vec;
+ mat.mv(ones_vec, mat_times_vec2);
+ for (size_t ii = 0; ii < rows; ++ii) {
+ EXPECT_DOUBLE_EQ(mat_times_vec[ii], std::accumulate(mat[ii].begin(), mat[ii].end(), 0.));
+ EXPECT_DOUBLE_EQ(mat_times_vec[ii], mat_times_vec2[ii]);
+ }
+ auto mat_scaled = mat;
+ mat_scaled *= 2.;
+ auto mat_scaled2 = mat * 2.;
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj = 0; jj < cols; ++jj)
+ EXPECT_DOUBLE_EQ(mat_scaled[ii][jj], mat_scaled2[ii][jj]);
+ SquareMatrixType square_mat{{0, 1}, {2, 3}};
+ SquareMatrixType expected_inv{{-1.5, 0.5}, {1, 0}};
+ auto inv_square_mat = square_mat;
+ inv_square_mat.invert();
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj = 0; jj < rows; ++jj)
+ EXPECT_DOUBLE_EQ(expected_inv[ii][jj], inv_square_mat[ii][jj]);
+ ScalarMatrixType fvec_constructed(Dune::XT::Common::FieldVector<double, 1>{1.});
+ ScalarMatrixType fvec_constructed2(Dune::FieldVector<double, 1>{1.});
+ ScalarMatrixType fvec_assigned = Dune::FieldVector<double, 1>{1.};
+ EXPECT_EQ(fvec_constructed, one);
+ EXPECT_EQ(fvec_constructed2, one);
+ EXPECT_EQ(fvec_assigned, one);
+ auto two = one * 2.;
+ EXPECT_EQ(two[0][0], 2.);
+
+
+#if 0
+template <class L, int ROWS, int COLS, class R>
+Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, ROWS, COLS>
+operator-(const Dune::FieldMatrix<L, ROWS, COLS>& left, const Dune::FieldMatrix<R, ROWS, COLS>& right)
+{
+ Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, ROWS, COLS> ret = left;
+ ret -= right;
+ return ret;
+}
+
+
+template <class L, int ROWS, int COLS, class R>
+Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, ROWS, COLS>
+operator+(const Dune::FieldMatrix<L, ROWS, COLS>& left, const Dune::FieldMatrix<R, ROWS, COLS>& right)
+{
+ Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, ROWS, COLS> ret = left;
+ ret += right;
+ return ret;
+}
+
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS> operator*(const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ return left.rightmultiplyany(right);
+}
+
+// we need this explicit overload to fix an ambiguous operator* error due to the automatic conversion from
+// FieldMatrix<K, 1, 1> to K
+template <class K, int L_ROWS>
+Dune::XT::Common::FieldMatrix<K, L_ROWS, 1> operator*(const Dune::XT::Common::FieldMatrix<K, L_ROWS, 1>& left,
+ const Dune::FieldMatrix<K, 1, 1>& right)
+{
+ return left.rightmultiplyany(right);
+}
+
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+void rightmultiply(Dune::FieldMatrix<K, L_ROWS, R_COLS>& ret,
+ const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ for (size_t ii = 0; ii < L_ROWS; ++ii) {
+ for (size_t jj = 0; jj < R_COLS; ++jj) {
+ ret[ii][jj] = 0.;
+ for (size_t kk = 0; kk < L_COLS; ++kk)
+ ret[ii][jj] += left[ii][kk] * right[kk][jj];
+ }
+ }
+}
+
+template <class L, int L_ROWS, int L_COLS, class R, int R_COLS>
+typename std::enable_if<!std::is_same<L, R>::value,
+ Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, L_ROWS, R_COLS>>::
+ type
+ operator*(const Dune::FieldMatrix<L, L_ROWS, L_COLS>& left, const Dune::FieldMatrix<R, L_COLS, R_COLS>& right)
+{
+ using Promoted = Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, L_ROWS, R_COLS>;
+ using Dune::XT::Common::convert_to;
+ return convert_to<Promoted>(left).rightmultiplyany(convert_to<Promoted>(right));
+}
+
+// versions that do not allocate matrices on the stack (for large matrices)
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
+ rightmultiply(*ret, *left, right);
+ return ret;
+}
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
+{
+ auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
+ rightmultiply(*ret, left, *right);
+ return ret;
+}
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
+ const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
+{
+ return left * *right;
+}
+#endif
+}
+
+GTEST_TEST(blockedfieldmatrix, creation_and_calculations)
+{
+ static constexpr size_t num_blocks = 2;
+ static constexpr size_t block_rows = 2;
+ static constexpr size_t block_cols = 3;
+ static constexpr size_t rows = num_blocks * block_rows;
+ static constexpr size_t cols = num_blocks * block_cols;
+ using BlockedMatrixType = Dune::XT::Common::BlockedFieldMatrix<double, num_blocks, block_rows, block_cols>;
+ using MatrixType = Dune::FieldMatrix<double, rows, cols>;
+ using XtMatrixType = Dune::XT::Common::FieldMatrix<double, rows, cols>;
+ using RowVectorType = Dune::XT::Common::FieldVector<double, cols>;
+ using ColVectorType = Dune::XT::Common::FieldVector<double, rows>;
+ using BlockedRowVectorType = Dune::XT::Common::BlockedFieldVector<double, num_blocks, block_cols>;
+ using BlockedColVectorType = Dune::XT::Common::BlockedFieldVector<double, num_blocks, block_rows>;
+ using BlockType = BlockedMatrixType::BlockType;
+ BlockedMatrixType zeros;
+ BlockedMatrixType ones(1.);
+ MatrixType mat({{0, 1, 2, 0, 0, 0}, {3, 4, 5, 0, 0, 0}, {0, 0, 0, 0, 1, 2}, {0, 0, 0, 3, 4, 5}});
+ BlockedMatrixType block_mat(mat);
+ BlockedMatrixType block_mat2 = mat;
+ BlockType block{{0, 1, 2}, {3, 4, 5}};
+ BlockedMatrixType block_mat3(block);
+ const BlockedMatrixType const_block_mat = block_mat;
+ EXPECT_EQ(block_mat, block_mat2);
+ EXPECT_EQ(block_mat, block_mat3);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ EXPECT_EQ(block_mat.block(jj), block);
+ EXPECT_EQ(const_block_mat.block(jj), block);
+ }
+ for (size_t ii = 0; ii < rows; ++ii) {
+ for (size_t jj = 0; jj < cols; ++jj) {
+ EXPECT_EQ(zeros.get_entry(ii, jj), 0.);
+ EXPECT_EQ(block_mat.get_entry(ii, jj), mat[ii][jj]);
+ EXPECT_EQ(const_block_mat.get_entry(ii, jj), mat[ii][jj]);
+ }
+ }
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ for (size_t ll = 0; ll < block_rows; ++ll) {
+ for (size_t mm = 0; mm < block_cols; ++mm) {
+ EXPECT_EQ(ones.block(jj)[ll][mm], 1.);
+ EXPECT_EQ(block_mat.get_entry(jj, ll, mm), mat[jj * block_rows + ll][jj * block_cols + mm]);
+ EXPECT_EQ(const_block_mat.get_entry(jj, ll, mm), mat[jj * block_rows + ll][jj * block_cols + mm]);
+ EXPECT_EQ(block_mat.block(jj)[ll][mm], mat[jj * block_rows + ll][jj * block_cols + mm]);
+ EXPECT_EQ(const_block_mat.block(jj)[ll][mm], mat[jj * block_rows + ll][jj * block_cols + mm]);
+ }
+ }
+ }
+ RowVectorType ones_vec(1.);
+ ColVectorType mat_times_vec;
+ BlockedRowVectorType block_ones_vec(1.);
+ BlockedColVectorType block_mat_times_vec;
+ block_mat.mv(ones_vec, mat_times_vec);
+ block_mat.mv(block_ones_vec, block_mat_times_vec);
+ for (size_t ii = 0; ii < rows; ++ii) {
+ EXPECT_DOUBLE_EQ(mat_times_vec[ii], std::accumulate(mat[ii].begin(), mat[ii].end(), 0.));
+ EXPECT_DOUBLE_EQ(mat_times_vec[ii], block_mat_times_vec[ii]);
+ }
+ ColVectorType ones_vec2(1.);
+ RowVectorType mtv_vec;
+ BlockedColVectorType block_ones_vec2(1.);
+ BlockedRowVectorType block_mtv_vec;
+ block_mat.mtv(ones_vec2, mtv_vec);
+ block_mat.mtv(block_ones_vec2, block_mtv_vec);
+ auto transposed_mat = XtMatrixType(mat).transpose();
+ auto transposed_block_mat = block_mat.transpose();
+ for (size_t ii = 0; ii < cols; ++ii) {
+ EXPECT_DOUBLE_EQ(mtv_vec[ii], std::accumulate(transposed_mat[ii].begin(), transposed_mat[ii].end(), 0.));
+ EXPECT_DOUBLE_EQ(mtv_vec[ii], block_mtv_vec[ii]);
+ }
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj = 0; jj < cols; ++jj)
+ EXPECT_EQ(transposed_block_mat.get_entry(jj, ii), transposed_mat[jj][ii]);
+ Dune::XT::Common::BlockedFieldMatrix<double, 2, 2, 2> block_square_mat;
+ size_t val = 0;
+ for (size_t jj = 0; jj < 2; ++jj)
+ for (size_t ll = 0; ll < 2; ++ll)
+ for (size_t mm = 0; mm < 2; ++mm)
+ block_square_mat.block(jj)[ll][mm] = val++;
+ auto block_square_mat_copy = block_square_mat;
+ block_square_mat.rightmultiply(block_square_mat_copy);
+ Dune::FieldMatrix<double, 4, 4> expected_mat{{2, 3, 0, 0}, {6, 11, 0, 0}, {0, 0, 46, 55}, {0, 0, 66, 79}};
+ for (size_t ii = 0; ii < 4; ++ii)
+ for (size_t jj = 0; jj < 4; ++jj)
+ EXPECT_EQ(block_square_mat.get_entry(ii, jj), expected_mat[ii][jj]);
+}
diff --git a/dune/xt/common/test/fvector.cc b/dune/xt/common/test/fvector.cc
index 78d9df81f5fba0caa0f8e3b9b392fc471b8723fd..73bdfff784e7c20fe2a8efd894f6458e56c3bc2c 100644
--- a/dune/xt/common/test/fvector.cc
+++ b/dune/xt/common/test/fvector.cc
@@ -11,7 +11,7 @@
#include <dune/xt/common/test/main.hxx> // <- Needs to come first, include the config.h.
#include <dune/xt/common/fvector.hh>
-
+#include <dune/xt/common/float_cmp.hh>
// using this function to allow for {...} syntax below (would not work in the EXPECT_EQ macro directly).
template <class A, int a, class B, int b>
@@ -20,6 +20,113 @@ void check(const Dune::XT::Common::FieldVector<A, a>& expected_result, const Dun
EXPECT_EQ(expected_result, actual_result);
}
+GTEST_TEST(dune_xt_common_field_vector, creation_and_calculations)
+{
+ using VectorType = Dune::XT::Common::FieldVector<double, 4>;
+ using BaseVectorType = Dune::FieldVector<double, 4>;
+ using RowMatrixType = Dune::FieldMatrix<double, 1, 4>;
+ using ColMatrixType = Dune::FieldMatrix<double, 4, 1>;
+ using ArrayType = std::array<double, 4>;
+ VectorType zeros;
+ for (const auto& val : zeros)
+ EXPECT_EQ(val, 0.);
+ VectorType test_vec{0, 1, 2, 3};
+ for (size_t ii = 0; ii < 4; ++ii)
+ EXPECT_EQ(test_vec[ii], ii);
+ EXPECT_ANY_THROW(VectorType({1, 2, 3}));
+ VectorType copy_assigned = test_vec;
+ check(copy_assigned, test_vec);
+ VectorType copy_constructed(test_vec);
+ check(copy_constructed, test_vec);
+ std::vector<double> std_vec{0, 1, 2, 3};
+ VectorType constructed_from_std_vec(std_vec);
+ check(constructed_from_std_vec, test_vec);
+ std::vector<double> std_vec2 = test_vec;
+ for (size_t ii = 0; ii < 4; ++ii)
+ EXPECT_EQ(std_vec2[ii], ii);
+ std_vec.resize(2);
+ EXPECT_ANY_THROW(VectorType{std_vec});
+ BaseVectorType base_vec{0, 1, 2, 3};
+ VectorType constructed_from_base(base_vec);
+ check(constructed_from_base, test_vec);
+ RowMatrixType row_mat = test_vec;
+ for (size_t ii = 0; ii < 4; ++ii)
+ EXPECT_EQ(row_mat[0][ii], ii);
+ // TODO: reenable this test
+ // ColMatrixType col_mat = test_vec;
+ // for (size_t ii = 0; ii < 4; ++ii)
+ // EXPECT_EQ(col_mat[ii][0], ii);
+ ArrayType array = test_vec;
+ for (size_t ii = 0; ii < 4; ++ii)
+ EXPECT_EQ(array[ii], ii);
+ // TODO: remove this when the construction above is reenabled
+ ColMatrixType col_mat;
+ for (size_t ii = 0; ii < 4; ++ii)
+ col_mat[ii][0] = ii;
+ EXPECT_DOUBLE_EQ(test_vec * test_vec, 14.);
+ EXPECT_DOUBLE_EQ(test_vec * col_mat, 14.);
+ for (size_t ii = 0; ii < 4; ++ii)
+ EXPECT_DOUBLE_EQ((test_vec * 2.)[ii], 2. * ii);
+ for (size_t ii = 0; ii < 4; ++ii)
+ EXPECT_DOUBLE_EQ((test_vec / 2.)[ii], ii / 2.);
+}
+
+
+GTEST_TEST(blockedfieldvector, creation_and_calculations)
+{
+ using BlockedVectorType = Dune::XT::Common::BlockedFieldVector<double, 2, 2>;
+ using VectorType = Dune::FieldVector<double, 4>;
+ using BlockType = typename BlockedVectorType::BlockType;
+ // using RowMatrixType = Dune::FieldMatrix<double, 1, 4>;
+ // using ColMatrixType = Dune::FieldMatrix<double, 4, 1>;
+ // using ArrayType = std::array<double, 4>;
+ BlockedVectorType zeros;
+ const BlockedVectorType const_zeros;
+ for (const auto& val : zeros)
+ EXPECT_EQ(val, 0.);
+ for (const auto& val : const_zeros)
+ EXPECT_EQ(val, 0.);
+ VectorType vec{0, 1, 2, 3};
+ BlockedVectorType vec_blocked(vec);
+ BlockedVectorType vec_blocked2 = vec;
+ const BlockedVectorType const_vec_blocked(vec);
+ BlockedVectorType vec_blocked3(BlockType{0., 1.});
+ for (size_t ii = 0; ii < 4; ++ii) {
+ EXPECT_EQ(vec_blocked[ii], ii);
+ EXPECT_EQ(vec_blocked2[ii], ii);
+ EXPECT_EQ(const_vec_blocked[ii], ii);
+ }
+ for (size_t jj = 0; jj < 2; ++jj) {
+ for (size_t ii = 0; ii < 2; ++ii) {
+ EXPECT_EQ(vec_blocked.block(jj)[ii], ii + 2 * jj);
+ EXPECT_EQ(vec_blocked.get_entry(jj, ii), ii + 2 * jj);
+ EXPECT_EQ(vec_blocked3.get_entry(jj, ii), ii);
+ EXPECT_EQ(const_vec_blocked.block(jj)[ii], ii + 2 * jj);
+ EXPECT_EQ(const_vec_blocked.get_entry(jj, ii), ii + 2 * jj);
+ }
+ }
+ EXPECT_DOUBLE_EQ(vec_blocked.one_norm(), 6);
+ EXPECT_DOUBLE_EQ(vec_blocked.two_norm(), std::sqrt(14));
+ EXPECT_DOUBLE_EQ(vec_blocked.two_norm2(), 14);
+ EXPECT_DOUBLE_EQ(vec_blocked * vec_blocked, 14);
+ VectorType assigned_vec = vec_blocked;
+ for (size_t ii = 0; ii < 4; ++ii)
+ EXPECT_EQ(assigned_vec[ii], vec_blocked[ii]);
+ EXPECT_EQ(*vec_blocked.data(), 0.);
+ EXPECT_TRUE(Dune::XT::Common::FloatCmp::eq(vec_blocked - vec_blocked, zeros));
+ EXPECT_TRUE(Dune::XT::Common::FloatCmp::eq(vec_blocked + (vec_blocked * (-1)), zeros));
+ EXPECT_TRUE(Dune::XT::Common::FloatCmp::eq(vec_blocked * 0., zeros));
+ auto vec_blocked_copy = vec_blocked;
+ vec_blocked_copy *= 0.;
+ EXPECT_TRUE(Dune::XT::Common::FloatCmp::eq(vec_blocked_copy, zeros));
+ vec_blocked_copy = vec_blocked;
+ vec_blocked_copy += vec_blocked;
+ for (size_t ii = 0; ii < 4; ++ii)
+ EXPECT_DOUBLE_EQ(vec_blocked_copy[ii], 2. * ii);
+ vec_blocked_copy -= vec_blocked;
+ EXPECT_TRUE(Dune::XT::Common::FloatCmp::eq(vec_blocked_copy, vec_blocked));
+}
+
GTEST_TEST(hstack, dune_field_vector)
{