Commit fe4e5d71 authored by Carsten Gräser's avatar Carsten Gräser

Cleanup symmetric matrix

* Move into Dune::Fufem
* Unify indentation
* Add typedef for stored data type
* Return const Data&
parent 4bad0670
Pipeline #16236 canceled with stage
......@@ -4,6 +4,7 @@
#include <dune/common/fvector.hh>
namespace Dune {
namespace Fufem {
/** \brief A class implementing a symmetric matrix with compile-time size
*
......@@ -19,91 +20,94 @@ public:
/** \brief The type used for scalars
*/
typedef T field_type;
/** \brief Default constructor, creates uninitialized matrix
*/
SymmetricMatrix()
{}
/** \brief Construct from raw memory */
template <class Iterator>
SymmetricMatrix(Iterator it)
{
for (size_t i=0; i<data_.size(); ++i, ++it)
data_[i] = *it;
}
SymmetricMatrix<T,N>& operator=(const T& s)
{
data_ = s;
return *this;
}
SymmetricMatrix<T,N>& 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<T,N>& v1, const FieldVector<T,N>& v2) const
{
T result = 0;
for (size_t i=0; i<N; i++)
for (size_t j=0; j<=i; j++)
result += (1-0.5*(i==j)) * operator()(i,j) * (v1[i] * v2[j] + v1[j] * v2[i]);
return result;
}
void axpy(const T& a, const SymmetricMatrix<T,N>& other)
{
data_.axpy(a,other.data_);
}
Dune::FieldVector<T,N*(N+1)/2> data() const
{
return data_;
}
Dune::FieldVector<T,N*(N+1)/2>& data()
{
return data_;
}
using field_type = T;
using Data = Dune::FieldVector<T,N*(N+1)/2>;
/** \brief Default constructor, creates uninitialized matrix
*/
SymmetricMatrix()
{}
/** \brief Construct from raw memory */
template <class Iterator>
SymmetricMatrix(Iterator it)
{
for (size_t i=0; i<data_.size(); ++i, ++it)
data_[i] = *it;
}
SymmetricMatrix<T,N>& operator=(const T& s)
{
data_ = s;
return *this;
}
SymmetricMatrix<T,N>& 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<T,N>& v1, const FieldVector<T,N>& v2) const
{
T result = 0;
for (size_t i=0; i<N; i++)
for (size_t j=0; j<=i; j++)
result += (1-0.5*(i==j)) * operator()(i,j) * (v1[i] * v2[j] + v1[j] * v2[i]);
return result;
}
void axpy(const T& a, const SymmetricMatrix<T,N>& other)
{
data_.axpy(a,other.data_);
}
const Data& data() const
{
return data_;
}
Data& data()
{
return data_;
}
private:
Dune::FieldVector<T,N*(N+1)/2> data_;
Data data_;
};
}
} // namespace Fufem
} // namespace Dune
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment