From 5b4326e37bbfce6b75b7fd9fc93dfb7c18f46339 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Thu, 28 Mar 2019 11:58:06 +0100 Subject: [PATCH] Implementation of symmetric matrices --- dune/fufem/symmetricmatrix.hh | 109 ++++++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 dune/fufem/symmetricmatrix.hh diff --git a/dune/fufem/symmetricmatrix.hh b/dune/fufem/symmetricmatrix.hh new file mode 100644 index 00000000..4fc92ca0 --- /dev/null +++ b/dune/fufem/symmetricmatrix.hh @@ -0,0 +1,109 @@ +#ifndef DUNE_FUFEM_SYMMETRICMATRIX_HH +#define DUNE_FUFEM_SYMMETRICMATRIX_HH + +#include + +namespace Dune { + +/** \brief A class implementing a symmetric matrix with compile-time size + * + * A \f$dim\times dim \f$ matrix is stored internally as a Dune::FieldVector + * The components are assumed to be \f$[ E(1,1),\ E(2,1),\ E(2,2),\ E(3,1),\ E(3,2),\ E(3,3) ]\f$ + * and analogous for other dimensions + * \tparam dim number of lines/columns of the matrix + */ +template +class SymmetricMatrix +{ +public: + + /** \brief The type used for scalars + */ + typedef T field_type; + + /** \brief Default constructor, creates uninitialized matrix + */ + SymmetricMatrix() + {} + + + /** \brief Construct from raw memory */ + template + SymmetricMatrix(Iterator it) + { + for (size_t i=0; i& operator=(const T& s) + { + data_ = s; + return *this; + } + + SymmetricMatrix& operator*=(const T& s) + { + data_ *= s; + return *this; + } + + /** \brief Matrix style random read/write access to components + * \param i line index + * \param j column index + * \note You need to know what you are doing: You can only access the lower + * left triangular submatrix using this methods. It requires i>=j. + */ + T& operator() (int i, int j) + { + if( i >= j ) + return data_[i*(i+1)/2+j]; + else + return data_[j*(j+1)/2+i]; + } + + /** \brief Matrix style random read access to components + * \param i line index + * \param j column index + * \note You need to know what you are doing: You can only access the lower + * left triangular submatrix using this methods. It requires i>=j. + */ + const T& operator() (int i, int j) const + { + if( i >= j ) + return data_[i*(i+1)/2+j]; + else + return data_[j*(j+1)/2+i]; + } + + T energyScalarProduct (const FieldVector& v1, const FieldVector& v2) const + { + T result = 0; + for (size_t i=0; i& other) + { + data_.axpy(a,other.data_); + } + + Dune::FieldVector data() const + { + return data_; + } + + Dune::FieldVector& data() + { + return data_; + } + +private: + Dune::FieldVector data_; +}; + +} +#endif + -- GitLab