diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt index c1f8cb1029b10f3d5669f6606609576a39aee1ac..d25563c693555f5879255ec97904dfaf8a8b0e21 100644 --- a/dune/matrix-vector/CMakeLists.txt +++ b/dune/matrix-vector/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(test) add_subdirectory(traits) +add_subdirectory(types) #install headers install(FILES diff --git a/dune/matrix-vector/types/CMakeLists.txt b/dune/matrix-vector/types/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..1feb2f3fc7c2de67efcef88d4c72c45fecfa3f0c --- /dev/null +++ b/dune/matrix-vector/types/CMakeLists.txt @@ -0,0 +1,5 @@ +#install headers +install(FILES + multitypematrix.hh + multitypevector.hh + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/matrix-vector/traits) diff --git a/dune/matrix-vector/types/multitypematrix.hh b/dune/matrix-vector/types/multitypematrix.hh new file mode 100644 index 0000000000000000000000000000000000000000..5e86decd44f341b390b736aa9ce4778bdd01b84d --- /dev/null +++ b/dune/matrix-vector/types/multitypematrix.hh @@ -0,0 +1,111 @@ +#ifndef DUNE_MATRIX_VECTOR_TYPES_MULTITYPEMATRIX_HH +#define DUNE_MATRIX_VECTOR_TYPES_MULTITYPEMATRIX_HH + +#include <dune/common/hybridutilities.hh> +#include <dune/common/typetraits.hh> +#include <dune/istl/multitypeblockmatrix.hh> +#include <dune/matrix-vector/traits/matrixtraits.hh> +#include <dune/matrix-vector/traits/utilities.hh> + +namespace Dune { +namespace MatrixVector { + +template <class...> +struct MultiTypeMatrix; + +template <class... Cols> +using MultiTypeMatrixRow = MultiTypeMatrix<Cols...>; + +/** + * \brief Is a dune-istl MultiTypeBlockMatrix with additional functionality: + * - is nestable in other dune matrix types, e.g. BCRSMatrix + * - construction from scalar + * - (recursive) infinity_norm() + * - scaling through multiplication and division by a scalar + */ +template <class... Rows> +struct MultiTypeMatrix : public MultiTypeBlockMatrix<Rows...> { + using This = MultiTypeMatrix<Rows...>; + using Base = MultiTypeBlockMatrix<Rows...>; + using Base::operator=; + using typename Base::field_type; + + // Fake value to silence matrices that extract this parameter. + constexpr static int blocklevel = -1000; + + //! \brief Default constructor + MultiTypeMatrix() {} + + //! \brief Constructor initializing all blocks with given scalar + template <class K, typename = EnableNumber<K>> + MultiTypeMatrix(K value) { + *this = value; + } + + template <class K, typename = EnableNumber<K>> + void operator*=(const K& w) { + Dune::Hybrid::forEach(*this, [&](auto&& entry) { entry *= w; }); + } + + template <class K, typename = EnableNumber<K>> + void operator/=(const K& w) { + Dune::Hybrid::forEach(*this, [&](auto&& entry) { entry /= w; }); + } + + // Note: This code is mostly copied from MultiTypeBlockVector. + // Note: We enforce lazy template instantiation to avoid compiler + // erros for nested types that do not provide infinity_norm() if + // this method is never called. + template <class = This> + auto infinity_norm() const { + using namespace Hybrid; + using std::max; + using real_type = real_t<This>; + + real_type result = 0.0; + // Compute max norm tracking appearing nan values + // if the field type supports nan. + ifElse(has_nan<field_type>(), + [&](auto) { + // This variable will preserve any nan value + real_type nanTracker = result; + forEach(*this, [&](auto&& entry) { + real_type entryNorm = entry.infinity_norm(); + result = max(entryNorm, result); + nanTracker += entryNorm; + }); + // Incorporate possible nan value into result + result *= (nanTracker / nanTracker); + }, + [&](auto) { + forEach(*this, [&](auto&& entry) { + result = max(entry.infinity_norm(), result); + }); + }); + + return result; + } +}; + +// inject matrix traits +namespace Traits { + +template <class... Rows> +struct MatrixTraits<MultiTypeMatrix<Rows...>> { + constexpr static bool isMatrix = true; +}; + +} // end namespace Traits + +} // end namespace MatrixVector + +// TODO promote type from all entries. +template <class Row, class... Rows> +struct FieldTraits<MatrixVector::MultiTypeMatrix<Row, Rows...>> { + using field_type = field_t<Row>; + using real_type = real_t<Row>; +}; + +} // end namespace Dune + +#endif // DUNE_MATRIX_VECTOR_TYPES_MULTITYPEMATRIX_HH diff --git a/dune/matrix-vector/types/multitypevector.hh b/dune/matrix-vector/types/multitypevector.hh new file mode 100644 index 0000000000000000000000000000000000000000..cdfdaa348d69e59b37c37f2d5a574deccc5db4f9 --- /dev/null +++ b/dune/matrix-vector/types/multitypevector.hh @@ -0,0 +1,63 @@ +#ifndef DUNE_MATRIX_VECTOR_TYPES_MULTITYPEVECTOR_HH +#define DUNE_MATRIX_VECTOR_TYPES_MULTITYPEVECTOR_HH + +#include <dune/common/typetraits.hh> +#include <dune/istl/multitypeblockvector.hh> +#include <dune/matrix-vector/traits/utilities.hh> +#include <dune/matrix-vector/traits/vectortraits.hh> + +namespace Dune { +namespace MatrixVector { + +/** + * \brief Is a dune-istl MultiTypeBlockVector with additional functionality: + * - is nestable in other dune matrix types, e.g. BlockVector + * - construction from scalar + * - scaling through division by scalar (multiplication is available) + */ +template <class... Args> +struct MultiTypeVector : public MultiTypeBlockVector<Args...> { + using Base = MultiTypeBlockVector<Args...>; + using Base::operator=; + + // Fake value to silence vectors that extract this parameter. + constexpr static int blocklevel = -1000; + + /** \brief Default constructor */ + MultiTypeVector() {} + + /** \brief Constructor initializing all blocks with given scalar */ + template <class K, typename = EnableNumber<K>> + MultiTypeVector(K value) { + *this = value; + } + + template <class K, typename = EnableNumber<K>> + void operator/=(const K& w) { + Dune::Hybrid::forEach(*this, [&](auto&& entry) { entry /= w; }); + } +}; + +// inject vector traits +namespace Traits { + +template <class... Args> +struct VectorTraits<MultiTypeVector<Args...>> { + constexpr static bool isVector = true; +}; + +} // end namespace Traits + +} // end namespace MatrixVector + +// inject field traits from vector +template <class Arg, class... Args> +struct FieldTraits<MatrixVector::MultiTypeVector<Arg, Args...>> { + // TODO generally one should promote the field type from all args. + using field_type = field_t<Arg>; + using real_type = real_t<Arg>; +}; + +} // end namespace Dune + +#endif // DUNE_MATRIX_VECTOR_TYPES_MULTITYPEVECTOR_HH