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
#ifndef DUNE_TENSOR_SSD_HH
#define DUNE_TENSOR_SSD_HH
/** \file
\brief A third-rank tensor with two static (SS) and one dynamic (D) dimension
*/
#include <dune/common/array.hh>
#include <dune/common/fmatrix.hh>
/** \brief A third-rank tensor with two static (SS) and one dynamic (D) dimension
*
* \tparam T Type of the entries
* \tparam N1 Size of the first dimension
* \tparam N2 Size of the second dimension
*/
template <class T, int N1, int N2>
class TensorSSD
{
public:
/** \brief Constructor with the third dimension */
explicit TensorSSD(size_t N3)
: N3_(N3)
{
for (int i=0; i<N1; i++)
for (int j=0; j<N2; j++)
data_[i][j].resize(N3_);
}
size_t dim(int index) const
{
switch (index) {
case 0:
return N1;
case 1:
return N2;
case 2:
return N3_;
default:
assert(false);
}
}
/** \brief Direct access to individual entries */
T& operator()(size_t i, size_t j, size_t k)
{
return data_[i][j][k];
}
/** \brief Direct const access to individual entries */
const T& operator()(size_t i, size_t j, size_t k) const
{
return data_[i][j][k];
}
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
/** \brief Assignment from scalar */
TensorSSD<T,N1,N2>& operator=(const T& scalar)
{
for (int i=0; i<N1; i++)
for (int j=0; j<N2; j++)
for (size_t k=0; k<dim(2); k++)
data_[i][j][k] = scalar;
return *this;
}
friend TensorSSD<T,N1,N2> operator*(const TensorSSD<T,N1,N2>& a, const Dune::Matrix<Dune::FieldMatrix<T,1,1> >& b)
{
TensorSSD<T,N1,N2> result(b.M());
assert(a.dim(2)==b.N());
size_t N4 = a.dim(2); // third dimension of a
for (int i=0; i<N1; i++)
for (int j=0; j<N2; j++)
for (size_t k=0; k<b.M(); k++) {
result.data_[i][j][k] = 0;
for (size_t l=0; l<N4; l++)
result.data_[i][j][k] += a.data_[i][j][l]*b[l][k][0][0];
}
return result;
}
friend TensorSSD<T,N1,N2> operator+(const TensorSSD<T,N1,N2>& a, const TensorSSD<T,N1,N2>& b)
{
assert(a.dim(2)==b.dim(2));
size_t N3 = a.dim(2);
TensorSSD<T,N1,N2> result(N3);
for (int i=0; i<N1; i++)
for (int j=0; j<N2; j++)
for (size_t k=0; k<N3; k++)
result.data_[i][j][k] = a.data_[i][j][k] + b.data_[i][j][k];
return result;
}
// having the dynamic data type on the inside is kind of a stupid data layout
Dune::array<Dune::array<std::vector<T>, N2>, N1> data_;
// size of the third dimension
size_t N3_;
};
#endif