Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#ifndef DUNE_GFE_SYMMETRICMATRIX_HH
#define DUNE_GFE_SYMMETRICMATRIX_HH
#include <dune/common/fvector.hh>
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 <tt>Dune::FieldVector<double, dim*(dim+1)/2></tt>
* 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 T, int N>
class SymmetricMatrix
{
public:
/** \brief Default constructor
*
* Tensor is initialized containing zeros if no argument is given.
* \param eye if true tensor is initialized as identity
*/
SymmetricMatrix()
{}
/** \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)
{
assert(i>=j);
return data_[i*(i+1)/2+j];
}
/** \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
{
assert(i>=j);
return data_[i*(i+1)/2+j];
}
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+(i!=j)) * operator()(i,j) * v1[i] * v2[j];
return result;
}
private:
Dune::FieldVector<T,N*(N+1)/2> data_;
};
}
#endif