diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt index cbfd6a59eb21069e86cff070fecde0685935e3af..50ae8466fe99ed16c3d1fa305d43c5e3a78ebd21 100644 --- a/dune/matrix-vector/CMakeLists.txt +++ b/dune/matrix-vector/CMakeLists.txt @@ -11,6 +11,7 @@ install(FILES axy.hh blockmatrixview.hh componentwisematrixmap.hh + componentwisevectormaps.hh concepts.hh crossproduct.hh genericvectortools.hh diff --git a/dune/matrix-vector/componentwisevectormaps.hh b/dune/matrix-vector/componentwisevectormaps.hh new file mode 100644 index 0000000000000000000000000000000000000000..fe80a4bf831e11a2facaf762dd639099a5f77e4a --- /dev/null +++ b/dune/matrix-vector/componentwisevectormaps.hh @@ -0,0 +1,44 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2 +#ifndef DUNE_MATRIX_VECTOR_COMPONENTWISE_VECTOR_MAPS +#define DUNE_MATRIX_VECTOR_COMPONENTWISE_VECTOR_MAPS +#include <dune/common/hybridutilities.hh> +#include <dune/common/typetraits.hh> +#include <dune/common/iteratorfacades.hh> + +namespace Dune { +namespace MatrixVector { + + /** Apply a function f componentwise to a vector v and store it in u, i.e. + * u_i = f(v_i) + */ + template<typename F, typename V> + void applyComponentwise(F&& f, V& u, const V& v) { + namespace H = Dune::Hybrid; + H::ifElse(Dune::IsNumber<std::decay_t<V>>(), [&](auto&& id) { + u= f(id(v));}, + [&](auto id) { + H::forEach(H::integralRange(H::size(id(v))), [&](auto i) { + applyComponentwise(f,H::elementAt(u,i), H::elementAt(v,i)); + }); + }); + } + + /** Multiply the i-th component of u by the i-th component of v, i.e. + * u_i <- u_i*v_i + */ + template<typename V> + void multiplyComponentwise(V& u, const V& v) { + namespace H = Dune::Hybrid; + H::ifElse(Dune::IsNumber<std::decay_t<V>>(), [&](auto&& id) { + u*=id(v);}, + [&](auto id) { + H::forEach(H::integralRange(H::size(id(v))), [&](auto i) { + multiplyComponentwise(H::elementAt(u,i), H::elementAt(v,i)); + }); + }); + } + +} +} +#endif diff --git a/dune/matrix-vector/test/CMakeLists.txt b/dune/matrix-vector/test/CMakeLists.txt index 4d2e82f6d77daf9294af6431241609d32433c1e2..698e3457eabb62c8e23098c5aafe6e7c6123eddd 100644 --- a/dune/matrix-vector/test/CMakeLists.txt +++ b/dune/matrix-vector/test/CMakeLists.txt @@ -1,4 +1,5 @@ dune_add_test(SOURCES arithmetictest.cc) +dune_add_test(SOURCES componentwisevectortest.cc) dune_add_test(SOURCES resizetest.cc) dune_add_test(SOURCES staticmatrixtoolstest.cc) dune_add_test(SOURCES triangularsolvetest.cc) diff --git a/dune/matrix-vector/test/componentwisevectortest.cc b/dune/matrix-vector/test/componentwisevectortest.cc new file mode 100644 index 0000000000000000000000000000000000000000..ebfec34027e1597b9ceb25fab9810f11add2d365 --- /dev/null +++ b/dune/matrix-vector/test/componentwisevectortest.cc @@ -0,0 +1,108 @@ +#include <config.h> +#include <dune/common/parallel/mpihelper.hh> +#include <dune/common/test/testsuite.hh> +#include <dune/common/indices.hh> + +#include <dune/istl/bvector.hh> +#include <dune/istl/vbvector.hh> +#include <dune/istl/multitypeblockvector.hh> + +#include <dune/matrix-vector/componentwisevectormaps.hh> + +#include <cmath> +#include <string> + +template<typename V> +auto test_vector(V v, std::string type) { + using namespace Dune; + auto suite = TestSuite("test componentwise on vector"); + + auto u = v; + + // test function application + auto f = [](const auto& in) { + return in*2.; + }; + + MatrixVector::applyComponentwise(f, u, v); + + // now, u_i should be f(v_i) = 2*v_i. + // we test this by computing the sum + // over both vectors. + // It should hold sum(u) == 2.*sum(v) (modulo floating point arithmetic) + auto ones =v; + ones = 1.; + + auto v_sum = ones*v; + + suite.check(std::abs(ones*u - 2.*v_sum) < 1e-15, "Test applyComponentwise") << "Failed for type " << type; + + // test componentwise multiply + u = 2.; // reset u's state + + MatrixVector::multiplyComponentwise(u, v); // u <- u.*v + + // test again sums + suite.check(std::abs(ones*u - 2.*v_sum) < 1e-15, "Test multiplyComponentwise") << "Failed for type " << type; + + return suite; + +} + +int main(int argc, char** argv) { + + auto& helper = Dune::MPIHelper::instance(argc, argv); + + auto suite = Dune::TestSuite("Test componentwise vector maps"); + + const auto n = 10; + + // Test BlockVector<FieldVector<double, 1>> + { + auto vector = Dune::BlockVector<Dune::FieldVector<double, 1>>(n); + vector = 0.5; + suite.subTest(test_vector(std::move(vector), "BlockVector<FieldVector<double, 1>>")); + } + + // Test BlockVector<FieldVector<double, 2>> + { + auto vector = Dune::BlockVector<Dune::FieldVector<double, 2>>(n); + vector = 0.5; + suite.subTest(test_vector(std::move(vector), "BlockVector<FieldVector<double, 2>>")); + } + + // Test BlockVector<double>> + { + auto vector = Dune::BlockVector<double>(n); + vector = 0.5; + suite.subTest(test_vector(std::move(vector), "BlockVector<double>")); + } + + // Test VariableBlockVector<double> + { + auto vector = Dune::VariableBlockVector<double>(n); + int i = 1; + for(auto it = vector.createbegin(); it != vector.createend(); ++it) { + *it = i; + i = i%2 +1; // non homogenous block structure + } + vector = 0.5; + suite.subTest(test_vector(std::move(vector), "VariableBlockVector<double>")); + } + + // Test MultiTypeBlockVector + { + using V1 = Dune::BlockVector<Dune::FieldVector<double, 1>>; + using V2 = Dune::BlockVector<Dune::FieldVector<double, 2>>; + + using MT = Dune::MultiTypeBlockVector<V1, V2>; + MT mt{}; + using namespace Dune::Indices; + mt[_0].resize(n); + mt[_1].resize(n); + mt = 0.5; + suite.subTest(test_vector(std::move(mt), "MultiTypeBlockVector")); + } + + return suite.exit(); +}