v0.13.1
Tensor3_value.hpp
Go to the documentation of this file.
1/* A general version, not for pointers. */
2
3#pragma once
4
6
7#include <ostream>
8
9namespace FTensor
10{
11 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
12 class Tensor3
13 {
14 T data[Tensor_Dim0][Tensor_Dim1][Tensor_Dim2];
15
16 public:
17 template <class... U> Tensor3(U... d) : data{d...}
18 {
19 static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
20 "Incorrect number of Arguments. Constructor should "
21 "initialize the entire Tensor");
22 }
23
25
26 /* There are two operator(int,int,int)'s, one for non-consts that lets you
27 change the value, and one for consts that doesn't. */
28
29 T &operator()(const int N1, const int N2, const int N3)
30 {
31#ifdef FTENSOR_DEBUG
32 if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim1 || N2 < 0
33 || N3 >= Tensor_Dim2 || N3 < 0)
34 {
35 std::stringstream s;
36 s << "Bad index in Tensor3<T," << Tensor_Dim0 << "," << Tensor_Dim1
37 << "," << Tensor_Dim2 << ">.operator(" << N1 << "," << N2 << ","
38 << N3 << ")" << std::endl;
39 throw std::out_of_range(s.str());
40 }
41#endif
42 return data[N1][N2][N3];
43 }
44
45 T operator()(const int N1, const int N2, const int N3) const
46 {
47#ifdef FTENSOR_DEBUG
48 if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim1 || N2 < 0
49 || N3 >= Tensor_Dim2 || N3 < 0)
50 {
51 std::stringstream s;
52 s << "Bad index in Tensor3<T," << Tensor_Dim0 << "," << Tensor_Dim1
53 << "," << Tensor_Dim2 << ">.operator(" << N1 << "," << N2 << ","
54 << N3 << ") const" << std::endl;
55 throw std::out_of_range(s.str());
56 }
57#endif
58 return data[N1][N2][N3];
59 }
60
61 /* These operator()'s are the first part in constructing template
62 expressions. */
63
64 template <char i, char j, char k, int Dim0, int Dim1, int Dim2>
65 typename std::enable_if<
66 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
68 Dim1, Dim2, i, j, k>>::type
69 operator()(const Index<i, Dim0>, const Index<j, Dim1>,
70 const Index<k, Dim2>)
71 {
73 Dim0, Dim1, Dim2, i, j, k>(*this);
74 }
75
76 template <char i, char j, char k, int Dim0, int Dim1, int Dim2>
77 typename std::enable_if<
78 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
80 Dim0, Dim1, Dim2, i, j, k>>::type
81 operator()(const Index<i, Dim0>, const Index<j, Dim1>,
82 const Index<k, Dim2>) const
83 {
84 return Tensor3_Expr<
86 Dim2, i, j, k>(*this);
87 }
88
89 /* These operators are for internal contractions. We can not do
90 something like A(i,j,j) where i and j have different dimensions,
91 because it becomes ambiguous. */
92 /*TODO I dont know how i and j having different dimensions can be
93 * ambiguous?? Since you do Use that in A(j,j,i)*/
94 /* A(i,j,j) */
95
96 template <char i, char j, int Dim>
97 typename std::enable_if<
98 (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 >= Dim),
102 T, Dim, i>>::type
103 operator()(const Index<i, Dim>, const Index<j, Dim>,
104 const Index<j, Dim>) const
105 {
106 using TensorExpr = Tensor3_contracted_12<
108 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
109 }
110
111 /* A(j,i,j) */
112
113 template <char i, char j, int Dim>
114 typename std::enable_if<
115 (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 >= Dim),
119 T, Dim, i>>::type
120 operator()(const Index<j, Dim>, const Index<i, Dim>,
121 const Index<j, Dim>) const
122 {
123 using TensorExpr = Tensor3_contracted_02<
125 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
126 }
127
128 /* A(j,j,i) */
129
130 template <char i, char j, int Dim, int Dim2>
131 typename std::enable_if<
132 (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 >= Dim2),
136 T, Dim2, i>>::type
137 operator()(const Index<j, Dim>, const Index<j, Dim>,
138 const Index<i, Dim2>) const
139 {
140 using TensorExpr = Tensor3_contracted_01<
142 return Tensor1_Expr<TensorExpr, T, Dim2, i>(TensorExpr(*this));
143 }
144
145 /* This is for expressions where a number is used for one slot, and
146 indices for the others, yielding a Tensor2_Expr or
147 Tensor2_symmetric_Expr. The non-const versions don't actually
148 create a Tensor3_number_rhs_* object, while the const versions
149 do create a Tensor3_number_*. */
150
151 /* First slot. */
152
153 template <char i, char j, int N, int Dim1, int Dim2>
154 typename std::enable_if<
155 (Tensor_Dim0 > N && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
158 T, Dim1, Dim2, i, j>>::type
159 operator()(const Number<N>, const Index<i, Dim1>, const Index<j, Dim2>)
160 {
161 using TensorExpr = Tensor3_number_rhs_0<
164 }
165
166 template <char i, char j, int N, int Dim1, int Dim2>
167 typename std::enable_if<
168 (Tensor_Dim0 > N && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
172 T, Dim1, Dim2, i, j>>::type
173 operator()(const Number<N>, const Index<i, Dim1>,
174 const Index<j, Dim2>) const
175 {
176 using TensorExpr = Tensor3_number_0<
178 return Tensor2_Expr<TensorExpr, T, Dim1, Dim2, i, j>(TensorExpr(*this));
179 }
180
181 /* Second slot. */
182
183 template <char i, char j, int N, int Dim0, int Dim2>
184 typename std::enable_if<
185 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 > N && Tensor_Dim2 >= Dim2),
188 T, Dim0, Dim2, i, j>>::type
189 operator()(const Index<i, Dim0>, const Number<N>, const Index<j, Dim2>)
190 {
191 using TensorExpr = Tensor3_number_rhs_0<
194 }
195
196 template <char i, char j, int N, int Dim0, int Dim2>
197 typename std::enable_if<
198 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 > N && Tensor_Dim2 >= Dim2),
202 T, Dim0, Dim2, i, j>>::type
203 operator()(const Index<i, Dim0>, const Number<N>,
204 const Index<j, Dim2>) const
205 {
206 using TensorExpr = Tensor3_number_0<
208 return Tensor2_Expr<TensorExpr, T, Dim0, Dim2, i, j>(TensorExpr(*this));
209 }
210
211 /* Third slot. */
212 // TODO allow two different dimensions here
213 template <char i, char j, int N, int Dim>
214 typename std::enable_if<
215 (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 > N),
218 T, N>,
219 T, Dim, i, j>>::type
220 operator()(const Index<i, Dim>, const Index<j, Dim>, const Number<N>)
221 {
222 using TensorExpr = Tensor3_number_rhs_2<
225 }
226 // TODO allow two different dimensions here
227 template <char i, char j, int N, int Dim>
228 typename std::enable_if<
229 (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim && Tensor_Dim2 > N),
233 T, Dim, i, j>>::type
234 operator()(const Index<i, Dim>, const Index<j, Dim>, const Number<N>) const
235 {
236 using TensorExpr = Tensor3_number_2<
239 TensorExpr(*this));
240 }
241
242 /* This is for expressions where a number is used for two slots, and
243 an Index for the other, yielding a Tensor1_Expr. The non-const
244 versions don't actually create a Tensor3_number_rhs_* object,
245 while the const versions do create a Tensor3_number_*. */
246
247 /* Index in first slot. */
248
249 template <char i, int N1, int N2, int Dim>
250 typename std::enable_if<
251 (Tensor_Dim0 >= Dim && Tensor_Dim1 > N1 && Tensor_Dim2 > N2),
255 T, Dim, i>>::type
256 operator()(const Index<i, Dim> index, const Number<N1>, const Number<N2>)
257 {
258 using TensorExpr = Tensor3_number_rhs_12<
261 }
262
263 template <char i, int N1, int N2, int Dim>
264 typename std::enable_if<
265 (Tensor_Dim0 >= Dim && Tensor_Dim1 > N1 && Tensor_Dim2 > N2),
269 T, Dim, i>>::type
270 operator()(const Index<i, Dim> index, const Number<N1>,
271 const Number<N2>) const
272 {
273 using TensorExpr = Tensor3_number_12<
275 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
276 }
277
278 /* Index in second slot. I use the same structures as for the Index
279 in the first slot since the tensor is symmetric on the first two
280 indices. */
281
282 template <char i, int N0, int N2, int Dim>
283 typename std::enable_if<
284 (Tensor_Dim0 > N0 && Tensor_Dim1 >= Dim && Tensor_Dim2 > N2),
288 T, Dim, i>>::type
289 operator()(const Number<N0>, const Index<i, Dim> index, const Number<N2>)
290 {
291 using TensorExpr = Tensor3_number_rhs_12<
294 }
295
296 template <char i, int N0, int N2, int Dim>
297 typename std::enable_if<
298 (Tensor_Dim0 > N0 && Tensor_Dim1 >= Dim && Tensor_Dim2 > N2),
302 T, Dim, i>>::type
303 operator()(const Number<N0>, const Index<i, Dim> index,
304 const Number<N2>) const
305 {
306 using TensorExpr = Tensor3_number_12<
308 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
309 }
310
311 /* Index in third slot. */
312
313 template <char i, int N0, int N1, int Dim>
314 typename std::enable_if<
315 (Tensor_Dim0 > N0 && Tensor_Dim1 > N1 && Tensor_Dim2 >= Dim),
319 T, Dim, i>>::type
320 operator()(const Number<N0>, const Number<N1>, const Index<i, Dim> index)
321 {
322 using TensorExpr = Tensor3_number_rhs_01<
325 }
326
327 template <char i, int N0, int N1, int Dim>
328 typename std::enable_if<
329 (Tensor_Dim0 > N0 && Tensor_Dim1 > N1 && Tensor_Dim2 >= Dim),
333 T, Dim, i>>::type
334 operator()(const Number<N0>, const Number<N1>,
335 const Index<i, Dim> index) const
336 {
337 using TensorExpr = Tensor3_number_01<
339 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
340 }
341
342 /* Specializations for using actual numbers instead of Number<>.
343 This is for expressions where an actual number is used for one
344 slot, and indices for the others, yielding a Tensor2_Expr or
345 Tensor2_symmetric_Expr. I only define the const versions. */
346
347 /* First slot. */
348
349 template <char i, char j, int Dim1, int Dim2>
350 typename std::enable_if<
351 (Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2),
354 T, Dim1, Dim2, i, j>>::type
355 operator()(const int N, const Index<i, Dim1>, const Index<j, Dim2>) const
356 {
357 using TensorExpr = Tensor3_numeral_0<
360 TensorExpr(*this, N));
361 }
362
363 /* Second slot. */
364
365 template <char i, char j, int Dim0, int Dim2>
366 typename std::enable_if<
367 (Tensor_Dim0 >= Dim0 && Tensor_Dim2 >= Dim2),
370 T, Dim0, Dim2, i, j>>::type
371 operator()(const Index<i, Dim0>, const int N, const Index<j, Dim2>) const
372 {
373 using TensorExpr = Tensor3_numeral_0<
376 TensorExpr(*this, N));
377 }
378
379 /* Third slot. */
380 // TODO Allow multiple dimensions here
381 template <char i, char j, int Dim>
382 typename std::enable_if<
383 (Tensor_Dim0 >= Dim && Tensor_Dim1 >= Dim),
387 T, Dim, i, j>>::type
388 operator()(const Index<i, Dim>, const Index<j, Dim>, const int N) const
389 {
390 using TensorExpr = Tensor3_numeral_2<
393 TensorExpr(*this, N));
394 }
395
396 /* This is for expressions where a numeral is used for two slots, and
397 an Index for the other, yielding a Tensor1_Expr. */
398
399 /* Index in first slot. */
400
401 template <char i, int Dim>
402 typename std::enable_if<
403 (Tensor_Dim0 >= Dim),
406 T, Dim, i>>::type
407 operator()(const Index<i, Dim> index, const int N1, const int N2) const
408 {
409 using TensorExpr = Tensor3_numeral_12<
411 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
412 }
413
414 /* Index in second slot. I use the same structures as for the Index
415 in the first slot since the tensor is symmetric on the first two
416 indices. */
417
418 template <char i, int Dim>
419 typename std::enable_if<
420 (Tensor_Dim1 >= Dim),
423 T, Dim, i>>::type
424 operator()(const int N1, const Index<i, Dim> index, const int N2) const
425 {
426 using TensorExpr = Tensor3_numeral_12<
428 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
429 }
430
431 /* Index in third slot. */
432
433 template <char i, int Dim>
434 typename std::enable_if<
435 (Tensor_Dim2 >= Dim),
438 T, Dim, i>>::type
439 operator()(const int N1, const int N2, const Index<i, Dim> index) const
440 {
441 using TensorExpr = Tensor3_numeral_01<
443 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
444 }
445 };
446}
447
448/// JSON compatible output for numbers. If T is something funky
449/// (strings, arrays, etc.) then you are going to have a bad time.
450
451namespace FTensor
452{
453 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
454 std::ostream &Tensor3_ostream_row(
455 std::ostream &os,
457 const int &i, const int &j)
458 {
459 os << '[';
460 for(int k = 0; k + 1 < Tensor_Dim2; ++k)
461 {
462 os << t(i, j, k) << ',';
463 }
464 if(Tensor_Dim2 > 0)
465 {
466 os << t(i, j, Tensor_Dim2 - 1);
467 }
468 os << ']';
469 return os;
470 }
471
472 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
474 std::ostream &os,
476 const int &i)
477 {
478 os << '[';
479 for(int j = 0; j + 1 < Tensor_Dim1; ++j)
480 {
482 os << ',';
483 }
484 if(Tensor_Dim1 > 0)
485 {
486 FTensor::Tensor3_ostream_row(os, t, i, Tensor_Dim1 - 1);
487 }
488 os << ']';
489 return os;
490 }
491}
492
493template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
494std::ostream &
495operator<<(std::ostream &os,
497{
498 os << '[';
499 for(int i = 0; i + 1 < Tensor_Dim0; ++i)
500 {
502 os << ',';
503 }
504 if(Tensor_Dim0 > 0)
505 {
506 FTensor::Tensor3_ostream_block(os, t, Tensor_Dim0 - 1);
507 }
508 os << ']';
509 return os;
510}
511
512/// JSON compatible input. It does not validate that separators are
513/// actually braces '[' or commas ','. It also ignores errors from
514/// missing trailing characters. So you could do something like
515///
516/// Tensor3<double,2,2,2> t3_1;
517/// std::stringstream ss(":::3:4:::7:8:::::11:12:::13:14");
518/// ss >> t3_1;
519
520namespace FTensor
521{
522 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
523 std::istream &Tensor3_istream_row(
524 std::istream &is,
526 const int &i, const int &j)
527 {
528 char c;
529 is >> c;
530 for(int k = 0; k + 1 < Tensor_Dim2; ++k)
531 {
532 is >> t(i, j, k) >> c;
533 }
534 if(Tensor_Dim2 > 0)
535 {
536 is >> t(i, j, Tensor_Dim2 - 1);
537 }
538 is >> c;
539 return is;
540 }
541
542 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
544 std::istream &is,
546 const int &i)
547 {
548 char c;
549 is >> c;
550 for(int j = 0; j + 1 < Tensor_Dim1; ++j)
551 {
552 Tensor3_istream_row(is, t, i, j);
553 is >> c;
554 }
555 if(Tensor_Dim1 > 0)
556 {
557 Tensor3_istream_row(is, t, i, Tensor_Dim1 - 1);
558 }
559 is >> c;
560 return is;
561 }
562}
563
564template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
565std::istream &
566operator>>(std::istream &is,
568{
569 char c;
570 is >> c;
571 for(int i = 0; i + 1 < Tensor_Dim0; ++i)
572 {
574 is >> c;
575 }
576 if(Tensor_Dim0 > 0)
577 {
578 FTensor::Tensor3_istream_block(is, t, Tensor_Dim0 - 1);
579 }
580 is >> c;
581 return is;
582}
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0
std::istream & operator>>(std::istream &is, FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t)
std::ostream & operator<<(std::ostream &os, const FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t)
T & operator()(const int N1, const int N2, const int N3)
T operator()(const int N1, const int N2, const int N3) const
T data[Tensor_Dim0][Tensor_Dim1][Tensor_Dim2]
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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 & Tensor3_ostream_block(std::ostream &os, const FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t, const int &i)
std::istream & Tensor3_istream_row(std::istream &is, FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t, const int &i, const int &j)
std::ostream & Tensor3_ostream_row(std::ostream &os, const FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t, const int &i, const int &j)
std::istream & Tensor3_istream_block(std::istream &is, FTensor::Tensor3< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > &t, const int &i)
constexpr double t
plate stiffness
Definition: plate.cpp:76
const int N
Definition: speed_test.cpp:3