v0.14.0
Loading...
Searching...
No Matches
Tensor2_symmetric_value.hpp
Go to the documentation of this file.
1#pragma once
2
3/* A general version, not for pointers. */
4
5#include <ostream>
6#ifdef FTENSOR_DEBUG
7#include <sstream>
8#include <stdexcept>
9#endif
10
11namespace FTensor
12{
13 template <class T, int Tensor_Dim> class Tensor2_symmetric
14 {
15 T data[(Tensor_Dim * (Tensor_Dim + 1)) / 2];
16
17 public:
18 template <class... U> Tensor2_symmetric(U... d) : data{d...}
19 {
20 static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
21 "Incorrect number of Arguments. Constructor should "
22 "initialize the entire Tensor");
23 }
24
26
27 /* There are two operator(int,int)'s, one for non-consts that lets you
28 change the value, and one for consts that doesn't. */
29
30 T &operator()(const int N1, const int N2)
31 {
32#ifdef FTENSOR_DEBUG
33 if(N1 >= Tensor_Dim || N1 < 0 || N2 >= Tensor_Dim || N2 < 0)
34 {
35 std::stringstream s;
36 s << "Bad index in Tensor2_symmetric<T," << Tensor_Dim
37 << ">.operator(" << N1 << "," << N2 << ")" << std::endl;
38 throw std::out_of_range(s.str());
39 }
40#endif
41 return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim - N2 - 1)) / 2]
42 : data[N2 + (N1 * (2 * Tensor_Dim - N1 - 1)) / 2];
43 }
44
45 T operator()(const int N1, const int N2) const
46 {
47#ifdef FTENSOR_DEBUG
48 if(N1 >= Tensor_Dim || N1 < 0 || N2 >= Tensor_Dim || N2 < 0)
49 {
50 std::stringstream s;
51 s << "Bad index in Tensor2_symmetric<T," << Tensor_Dim
52 << ">.operator(" << N1 << "," << N2 << ") const" << std::endl;
53 throw std::out_of_range(s.str());
54 }
55#endif
56 return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim - N2 - 1)) / 2]
57 : data[N2 + (N1 * (2 * Tensor_Dim - N1 - 1)) / 2];
58 }
59
60 template <int N1, int N2>
61 inline T operator()(const Number<N1> &, const Number<N2>()) const {
62
63 static_assert(N1 < Tensor_Dim, "Bad index");
64 static_assert(N2 < Tensor_Dim, "Bad index");
65
66 if constexpr (N1 > N2)
67 return data[N1 + (N2 * (2 * Tensor_Dim - N2 - 1)) / 2];
68 else
69 return data[N2 + (N1 * (2 * Tensor_Dim - N1 - 1)) / 2];
70 }
71
72 template <int N1, int N2>
73 inline T &operator()(const Number<N1> &, const Number<N2>()) {
74
75 static_assert(N1 < Tensor_Dim, "Bad index");
76 static_assert(N2 < Tensor_Dim, "Bad index");
77
78 if constexpr (N1 > N2)
79 return data[N1 + (N2 * (2 * Tensor_Dim - N2 - 1)) / 2];
80 else
81 return data[N2 + (N1 * (2 * Tensor_Dim - N1 - 1)) / 2];
82 }
83
84 /* These operator()'s are the first part in constructing template
85 expressions. They can be used to slice off lower dimensional
86 parts. They are not entirely safe, since you can accidentally use a
87 higher dimension than what is really allowed (like Dim=5). */
88
89 /* This returns a Tensor2_Expr, since the indices are not really
90 symmetric anymore since they cover different dimensions. */
91
92 template <char i, char j, int Dim0, int Dim1>
93 typename std::enable_if<
94 (Tensor_Dim >= Dim0 && Tensor_Dim >= Dim1),
95 Tensor2_Expr<Tensor2_symmetric<T, Tensor_Dim>, T, Dim0, Dim1, i, j>>::type
96 operator()(const Index<i, Dim0>, const Index<j, Dim1>)
97 {
99 j>(*this);
100 }
101
102 template <char i, char j, int Dim0, int Dim1>
103 typename std::enable_if<(Tensor_Dim >= Dim0 && Tensor_Dim >= Dim1),
105 T, Dim0, Dim1, i, j>>::type
106 operator()(const Index<i, Dim0>, const Index<j, Dim1>) const
107 {
109 Dim1, i, j>(*this);
110 }
111
112 /* This returns a Tensor2_symmetric_Expr, since the indices are still
113 symmetric on the lower dimensions. */
114
115 template <char i, char j, int Dim>
116 typename std::enable_if<
117 (Tensor_Dim >= Dim),
119 operator()(const Index<i, Dim>, const Index<j, Dim>)
120 {
122 i, j>(*this);
123 }
124
125 template <char i, char j, int Dim>
126 typename std::enable_if<
127 (Tensor_Dim >= Dim),
129 j>>::type
130 operator()(const Index<i, Dim>, const Index<j, Dim>) const
131 {
133 Dim, i, j>(*this);
134 }
135
136 /* This is for expressions where a number is used for one slot, and
137 an index for another, yielding a Tensor1_Expr. The non-const
138 versions don't actually create a Tensor2_number_rhs_[01] object.
139 They create a Tensor1_Expr directly, which provides the
140 appropriate indexing operators. The const versions do create a
141 Tensor2_number_[01]. */
142
143 template <char i, int N, int Dim>
144 typename std::enable_if<
145 (Tensor_Dim >= Dim && Tensor_Dim > N),
147 T, Dim, i>>::type
148 operator()(const Index<i, Dim>, const Number<N> &)
149 {
150 using TensorExpr
153 }
154
155 template <char i, int N, int Dim>
156 typename std::enable_if<
157 (Tensor_Dim >= Dim && Tensor_Dim > N),
159 T, Dim, i>>::type
160 operator()(const Index<i, Dim>, const Number<N> &) const
161 {
162 using TensorExpr
164 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
165 }
166
167 template <char i, int N, int Dim>
168 typename std::enable_if<
169 (Tensor_Dim > N && Tensor_Dim >= Dim),
171 T, Dim, i>>::type
172 operator()(const Number<N> &, const Index<i, Dim>)
173 {
174 using TensorExpr
177 }
178
179 template <char i, int N, int Dim>
180 typename std::enable_if<
181 (Tensor_Dim > N && Tensor_Dim >= Dim),
183 T, Dim, i>>::type
184 operator()(const Number<N> &, const Index<i, Dim>) const
185 {
186 using TensorExpr
188 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
189 }
190
191 /* Specializations for using actual numbers instead of Number<> */
192
193 template <char i, int Dim>
194 typename std::enable_if<
195 (Tensor_Dim >= Dim),
197 T, Dim, i>>::type
198 operator()(const Index<i, Dim>, const int N) const
199 {
200 using TensorExpr
202 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N));
203 }
204
205 template <char i, int Dim>
206 typename std::enable_if<
207 (Tensor_Dim >= Dim),
209 T, Dim, i>>::type
210 operator()(const int N, const Index<i, Dim>) const
211 {
212 using TensorExpr
214 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N));
215 }
216
217 /* These two operator()'s return the Tensor2 with internal
218 contractions, yielding a T. I have to specify one for both
219 const and non-const because otherwise the compiler will use the
220 operator() which gives a Tensor2_Expr<>. */
221
222 template <char i, int Dim>
223 typename std::enable_if<(Tensor_Dim >= Dim), T>::type
224 operator()(const Index<i, Dim>, const Index<i, Dim>)
225 {
227 }
228
229 template <char i, int Dim>
230 typename std::enable_if<(Tensor_Dim >= Dim), T>::type
231 operator()(const Index<i, Dim>, const Index<i, Dim>) const
232 {
234 }
235
236 private:
237 template <int N> T internal_contract(const Number<N> &) const
238 {
239 return data[N - 1 + ((N - 1) * (2 * Tensor_Dim - N)) / 2]
241 }
242
243 T internal_contract(const Number<1> &) const { return data[0]; }
244 };
245}
246
247/// JSON compatible output. It only outputs unique elements, so a 3x3
248/// symmetric matrix only outputs 6 elements.
249
250namespace FTensor
251{
252 template <class T, int Tensor_Dim>
254 std::ostream &os, const FTensor::Tensor2_symmetric<T, Tensor_Dim> &t,
255 const int &i)
256 {
257 os << '[';
258 for(int j = i; j + 1 < Tensor_Dim; ++j)
259 {
260 os << t(i, j) << ',';
261 }
262 if(Tensor_Dim > 0)
263 {
264 os << t(i, Tensor_Dim - 1);
265 }
266 os << ']';
267 return os;
268 }
269
270 template <class T, int Tensor_Dim>
271 std::ostream &operator<<(std::ostream &os,
273 os << '[';
274 for (int i = 0; i + 1 < Tensor_Dim; ++i) {
276 os << ',';
277 }
278 if (Tensor_Dim > 0) {
279 FTensor::Tensor2_symmetric_ostream_row(os, t, Tensor_Dim - 1);
280 }
281 os << ']';
282 return os;
283 }
284}
285
286namespace FTensor
287{
288 template <class T, int Tensor_Dim>
289 std::istream &
292 const int &i)
293 {
294 char c;
295 is >> c;
296 for(int j = i; j + 1 < Tensor_Dim; ++j)
297 {
298 is >> t(i, j) >> c;
299 }
300 if(Tensor_Dim > 0)
301 {
302 is >> t(i, Tensor_Dim - 1);
303 }
304 is >> c;
305 return is;
306 }
307
308 template <class T, int Tensor_Dim>
309 std::istream &operator>>(std::istream &is,
311 char c;
312 is >> c;
313 for (int i = 0; i + 1 < Tensor_Dim; ++i) {
315 is >> c;
316 }
317 if (Tensor_Dim > 0) {
318 FTensor::Tensor2_symmetric_istream_row(is, t, Tensor_Dim - 1);
319 }
320 is >> c;
321 return is;
322 }
323}
static Number< 2 > N2
static Number< 1 > N1
T operator()(const int N1, const int N2) const
T & operator()(const int N1, const int N2)
T data[(Tensor_Dim *(Tensor_Dim+1))/2]
T internal_contract(const Number< 1 > &) const
T operator()(const Number< N1 > &, const Number< N2 >()) const
T & operator()(const Number< N1 > &, const Number< N2 >())
T internal_contract(const Number< N > &) const
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'j', 3 > j
const double T
Tensors class implemented by Walter Landry.
Definition FTensor.hpp:51
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition dTensor0.hpp:27
std::ostream & Tensor2_symmetric_ostream_row(std::ostream &os, const FTensor::Tensor2_symmetric< T, Tensor_Dim > &t, const int &i)
std::istream & operator>>(std::istream &is, FTensor::Tensor1< T, Tensor_Dim > &t)
std::ostream & operator<<(std::ostream &os, const FTensor::Tensor1< T, Tensor_Dim > &t)
std::istream & Tensor2_symmetric_istream_row(std::istream &is, FTensor::Tensor2_symmetric< T, Tensor_Dim > &t, const int &i)
constexpr double t
plate stiffness
Definition plate.cpp:59
const int N
Definition speed_test.cpp:3