Skip to content
Snippets Groups Projects
Commit 88b250be authored by Max Kahnt's avatar Max Kahnt
Browse files

Add MultiType types that extend the dune core functionality.

parent c0d7b6e0
Branches
No related tags found
No related merge requests found
Pipeline #
add_subdirectory(test) add_subdirectory(test)
add_subdirectory(traits) add_subdirectory(traits)
add_subdirectory(types)
#install headers #install headers
install(FILES install(FILES
......
#install headers
install(FILES
multitypematrix.hh
multitypevector.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/matrix-vector/traits)
#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
#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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment