v0.13.2
Loading...
Searching...
No Matches
Dg_value.hpp
Go to the documentation of this file.
1/* A general version, not for pointers. */
2
3#pragma once
4
5#include "../Tensor3.hpp"
6
7namespace FTensor
8{
9 template <class T, int Tensor_Dim01, int Tensor_Dim2> class Dg
10 {
11 T data[(Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2][Tensor_Dim2];
12
13 public:
14 template <class... U> Dg(U... d) : data{d...}
15 {
16 static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
17 "Incorrect number of Arguments. Constructor should "
18 "initialize the entire Tensor");
19 }
20
21 Dg() {}
22
23 /* There are two operator(int,int,int)'s, one for non-consts that lets you
24 change the value, and one for consts that doesn't. */
25
26 T &operator()(const int N1, const int N2, const int N3)
27 {
28#ifdef FTENSOR_DEBUG
29 if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
30 || N3 >= Tensor_Dim2 || N3 < 0)
31 {
32 std::stringstream s;
33 s << "Bad index in Dg<T," << Tensor_Dim01 << "," << Tensor_Dim2
34 << ">.operator(" << N1 << "," << N2 << "," << N3 << ")"
35 << std::endl;
36 throw std::out_of_range(s.str());
37 }
38#endif
39 return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2][N3]
40 : data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2][N3];
41 }
42
43 T operator()(const int N1, const int N2, const int N3) const
44 {
45#ifdef FTENSOR_DEBUG
46 if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
47 || N3 >= Tensor_Dim2 || N3 < 0)
48 {
49 std::stringstream s;
50 s << "Bad index in Dg<T," << Tensor_Dim01 << "," << Tensor_Dim2
51 << ">.operator(" << N1 << "," << N2 << "," << N3 << ") const"
52 << std::endl;
53 throw std::out_of_range(s.str());
54 }
55#endif
56 return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2][N3]
57 : data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2][N3];
58 }
59
60 /* These operator()'s are the first part in constructing template
61 expressions. */
62
63 template <char i, char j, char k, int Dim01, int Dim2>
64 typename std::enable_if<
65 (Tensor_Dim01 >= Dim01 && Tensor_Dim2 >= Dim2),
66 Dg_Expr<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim01, Dim2, i, j, k>>::type
67 operator()(const Index<i, Dim01>, const Index<j, Dim01>,
68 const Index<k, Dim2>)
69 {
70 return Dg_Expr<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim01, Dim2, i, j, k>(
71 *this);
72 }
73 // TODO Add support for multiple dimensions ^^^^^
74 template <char i, char j, char k, int Dim01, int Dim2>
75 typename std::enable_if<(Tensor_Dim01 >= Dim01 && Tensor_Dim2 >= Dim2),
76 Dg_Expr<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T,
77 Dim01, Dim2, i, j, k>>::type
78 operator()(const Index<i, Dim01>, const Index<j, Dim01>,
79 const Index<k, Dim2>) const
80 {
81 return Dg_Expr<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim01, Dim2, i,
82 j, k>(*this);
83 }
84
85 /* These operators are for internal contractions. We can not do
86 something like A(i,j,j) where i and j have different dimensions,
87 because it becomes ambiguous. */
88
89 /* A(i,j,j) */
90
91 template <char i, char j, int Dim>
92 typename std::enable_if<
93 (Tensor_Dim01 >= Dim && Tensor_Dim2 >= Dim),
94 Tensor1_Expr<
95 Tensor3_contracted_12<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim>, T,
96 Dim, i>>::type
97 operator()(const Index<i, Dim>, const Index<j, Dim>,
98 const Index<j, Dim>) const
99 {
100 using TensorExpr
101 = Tensor3_contracted_12<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim>;
102 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
103 }
104
105 /* A(j,i,j) */
106
107 template <char i, char j, int Dim>
108 typename std::enable_if<
109 (Tensor_Dim01 >= Dim && Tensor_Dim2 >= Dim),
110 Tensor1_Expr<
111 Tensor3_contracted_02<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim>, T,
112 Dim, i>>::type
113 operator()(const Index<j, Dim>, const Index<i, Dim>,
114 const Index<j, Dim>) const
115 {
116 using TensorExpr
117 = Tensor3_contracted_02<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim>;
118 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
119 }
120
121 /* A(j,j,i) */
122
123 template <char i, char j, int Dim01, int Dim2>
124 typename std::enable_if<
125 (Tensor_Dim01 >= Dim01 && Tensor_Dim2 >= Dim2),
126 Tensor1_Expr<
127 Tensor3_contracted_01<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim01>, T,
128 Dim2, i>>::type
129 operator()(const Index<j, Dim01>, const Index<j, Dim01>,
130 const Index<i, Dim2>) const
131 {
132 using TensorExpr
133 = Tensor3_contracted_01<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim01>;
134 return Tensor1_Expr<TensorExpr, T, Dim2, i>(TensorExpr(*this));
135 }
136
137 /* This is for expressions where a number is used for one slot, and
138 indices for the others, yielding a Tensor2_Expr or
139 Tensor2_symmetric_Expr. The non-const versions don't actually
140 create a Dg_number_rhs_* object, while the const versions
141 do create a Dg_number_*. */
142
143 /* First slot. */
144
145 template <char i, char j, int N, int Dim1, int Dim2>
146 typename std::enable_if<
147 (Tensor_Dim01 > N && Tensor_Dim01 >= Dim1 && Tensor_Dim2 >= Dim2),
148 Tensor2_Expr<Dg_number_rhs_0<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>, T,
149 Dim1, Dim2, i, j>>::type
150 operator()(const Number<N>, const Index<i, Dim1>, const Index<j, Dim2>)
151 {
152 using TensorExpr
153 = Dg_number_rhs_0<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>;
154 return Tensor2_Expr<TensorExpr, T, Dim1, Dim2, i, j>(*this);
155 }
156
157 template <char i, char j, int N, int Dim1, int Dim2>
158 typename std::enable_if<
159 (Tensor_Dim01 > N && Tensor_Dim01 >= Dim1 && Tensor_Dim2 >= Dim2),
160 Tensor2_Expr<Dg_number_0<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>,
161 T, Dim1, Dim2, i, j>>::type
162 operator()(const Number<N>, const Index<i, Dim1>,
163 const Index<j, Dim2>) const
164 {
165 using TensorExpr
166 = Dg_number_0<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>;
167 return Tensor2_Expr<TensorExpr, T, Dim1, Dim2, i, j>(TensorExpr(*this));
168 }
169
170 /* Second slot. */
171
172 template <char i, char j, int N, int Dim0, int Dim2>
173 typename std::enable_if<
174 (Tensor_Dim01 >= Dim0 && Tensor_Dim01 > N && Tensor_Dim2 >= Dim2),
175 Tensor2_Expr<Dg_number_rhs_0<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>, T,
176 Dim0, Dim2, i, j>>::type
177 operator()(const Index<i, Dim0>, const Number<N>, const Index<j, Dim2>)
178 {
179 using TensorExpr
180 = Dg_number_rhs_0<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>;
181 return Tensor2_Expr<TensorExpr, T, Dim0, Dim2, i, j>(*this);
182 }
183
184 template <char i, char j, int N, int Dim0, int Dim2>
185 typename std::enable_if<
186 (Tensor_Dim01 >= Dim0 && Tensor_Dim01 > N && Tensor_Dim2 >= Dim2),
187 Tensor2_Expr<Dg_number_0<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>,
188 T, Dim0, Dim2, i, j>>::type
189 operator()(const Index<i, Dim0>, const Number<N>,
190 const Index<j, Dim2>) const
191 {
192 using TensorExpr
193 = Dg_number_0<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>;
194 return Tensor2_Expr<TensorExpr, T, Dim0, Dim2, i, j>(TensorExpr(*this));
195 }
196
197 /* Third slot. */
198
199 template <char i, char j, int N, int Dim>
200 typename std::enable_if<
201 (Tensor_Dim01 >= Dim && Tensor_Dim2 > N),
202 Tensor2_symmetric_Expr<
203 Dg_number_rhs_2<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>, T, Dim, i,
204 j>>::type
205 operator()(const Index<i, Dim>, const Index<j, Dim>, const Number<N>)
206 {
207 using TensorExpr
208 = Dg_number_rhs_2<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>;
209 return Tensor2_symmetric_Expr<TensorExpr, T, Dim, i, j>(*this);
210 }
211
212 template <char i, char j, int N, int Dim>
213 typename std::enable_if<
214 (Tensor_Dim01 >= Dim && Tensor_Dim2 > N),
215 Tensor2_symmetric_Expr<
216 Dg_number_2<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>, T, Dim, i,
217 j>>::type
218 operator()(const Index<i, Dim>, const Index<j, Dim>, const Number<N>) const
219 {
220 using TensorExpr
221 = Dg_number_2<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N>;
222 return Tensor2_symmetric_Expr<TensorExpr, T, Dim, i, j>(
223 TensorExpr(*this));
224 }
225 // TODO add multiple dimensions up here ^^^^^
226 /* This is for expressions where a number is used for two slots, and
227 an Index for the other, yielding a Tensor1_Expr. The non-const
228 versions don't actually create a Dg_number_rhs_* object,
229 while the const versions do create a Dg_number_*. */
230
231 /* Index in first slot. */
232
233 template <char i, int N1, int N2, int Dim>
234 typename std::enable_if<
235 (Tensor_Dim01 >= Dim && Tensor_Dim01 > N1 && Tensor_Dim2 > N2),
236 Tensor1_Expr<Dg_number_rhs_12<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N1, N2>,
237 T, Dim, i>>::type
238 operator()(const Index<i, Dim> index, const Number<N1>, const Number<N2>)
239 {
240 using TensorExpr
241 = Dg_number_rhs_12<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N1, N2>;
242 return Tensor1_Expr<TensorExpr, T, Dim, i>(*this);
243 }
244
245 template <char i, int N1, int N2, int Dim>
246 typename std::enable_if<
247 (Tensor_Dim01 >= Dim && Tensor_Dim01 > N1 && Tensor_Dim2 > N2),
248 Tensor1_Expr<
249 Dg_number_12<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N1, N2>, T,
250 Dim, i>>::type
251 operator()(const Index<i, Dim> index, const Number<N1>,
252 const Number<N2>) const
253 {
254 using TensorExpr
255 = Dg_number_12<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N1, N2>;
256 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
257 }
258
259 /* Index in second slot. I use the same structures as for the Index
260 in the first slot since the tensor is symmetric on the first two
261 indices. */
262
263 template <char i, int N0, int N2, int Dim>
264 typename std::enable_if<
265 (Tensor_Dim01 > N0 && Tensor_Dim01 >= Dim && Tensor_Dim2 > N2),
266 Tensor1_Expr<Dg_number_rhs_12<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N0, N2>,
267 T, Dim, i>>::type
268 operator()(const Number<N0>, const Index<i, Dim>, const Number<N2>)
269 {
270 using TensorExpr
271 = Dg_number_rhs_12<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N0, N2>;
272 return Tensor1_Expr<TensorExpr, T, Dim, i>(*this);
273 }
274
275 template <char i, int N0, int N2, int Dim>
276 typename std::enable_if<
277 (Tensor_Dim01 > N0 && Tensor_Dim01 >= Dim && Tensor_Dim2 > N2),
278 Tensor1_Expr<
279 Dg_number_12<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N0, N2>, T,
280 Dim, i>>::type
281 operator()(const Number<N0>, const Index<i, Dim> index,
282 const Number<N2>) const
283 {
284 using TensorExpr
285 = Dg_number_12<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N0, N2>;
286 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
287 }
288
289 /* Index in third slot. */
290
291 template <char i, int N0, int N1, int Dim>
292 typename std::enable_if<
293 (Tensor_Dim01 > N0 && Tensor_Dim01 > N1 && Tensor_Dim2 >= Dim),
294 Tensor1_Expr<Dg_number_rhs_01<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N0, N1>,
295 T, Dim, i>>::type
296 operator()(const Number<N0>, const Number<N1>, const Index<i, Dim> index)
297 {
298 using TensorExpr
299 = Dg_number_rhs_01<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N0, N1>;
300 return Tensor1_Expr<TensorExpr, T, Dim, i>(*this);
301 }
302
303 template <char i, int N0, int N1, int Dim>
304 typename std::enable_if<
305 (Tensor_Dim01 > N0 && Tensor_Dim01 > N1 && Tensor_Dim2 >= Dim),
306 Tensor1_Expr<
307 Dg_number_01<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N0, N1>, T,
308 Dim, i>>::type
309 operator()(const Number<N0>, const Number<N1>, const Index<i, Dim>) const
310 {
311 using TensorExpr
312 = Dg_number_01<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T, N0, N1>;
313 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
314 }
315
316 /* Specializations for using actual numbers instead of Number<>.
317 This is for expressions where an actual number is used for one
318 slot, and indices for the others, yielding a Tensor2_Expr or
319 Tensor2_symmetric_Expr. I only define the const versions. */
320
321 /* First slot. */
322
323 template <char i, char j, int Dim1, int Dim2>
324 typename std::enable_if<
325 (Tensor_Dim01 >= Dim1 && Tensor_Dim2 >= Dim2),
326 Tensor2_Expr<Dg_numeral_0<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>, T,
327 Dim1, Dim2, i, j>>::type
328 operator()(const int N, const Index<i, Dim1>, const Index<j, Dim2>) const
329 {
330 using TensorExpr
331 = Dg_numeral_0<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>;
332 return Tensor2_Expr<TensorExpr, T, Dim1, Dim2, i, j>(
333 TensorExpr(*this, N));
334 }
335
336 /* Second slot. */
337
338 template <char i, char j, int Dim0, int Dim2>
339 typename std::enable_if<
340 (Tensor_Dim01 >= Dim0 && Tensor_Dim2 >= Dim2),
341 Tensor2_Expr<Dg_numeral_0<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>, T,
342 Dim0, Dim2, i, j>>::type
343 operator()(const Index<i, Dim0>, const int N, const Index<j, Dim2>) const
344 {
345 using TensorExpr
346 = Dg_numeral_0<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>;
347 return Tensor2_Expr<TensorExpr, T, Dim0, Dim2, i, j>(
348 TensorExpr(*this, N));
349 }
350
351 /* Third slot. */
352
353 template <char i, char j, int Dim>
354 typename std::enable_if<
355 (Tensor_Dim01 >= Dim),
356 Tensor2_symmetric_Expr<
357 Dg_numeral_2<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>, T, Dim, i,
358 j>>::type
359 operator()(const Index<i, Dim>, const Index<j, Dim>, const int N) const
360 {
361 using TensorExpr
362 = Dg_numeral_2<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>;
363 return Tensor2_symmetric_Expr<TensorExpr, T, Dim, i, j>(
364 TensorExpr(*this, N));
365 }
366 // TODO Allow two different dimensions in 1 and 2
367 /* This is for expressions where a numeral is used for two slots, and
368 an Index for the other, yielding a Tensor1_Expr. */
369
370 /* Index in first slot. */
371
372 template <char i, int Dim>
373 typename std::enable_if<
374 (Tensor_Dim01 >= Dim),
375 Tensor1_Expr<Dg_numeral_12<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>, T,
376 Dim, i>>::type
377 operator()(const Index<i, Dim> index, const int N1, const int N2) const
378 {
379 using TensorExpr
380 = Dg_numeral_12<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>;
381 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
382 }
383
384 /* Index in second slot. I use the same structures as for the Index
385 in the first slot since the tensor is symmetric on the first two
386 indices. */
387
388 template <char i, int Dim>
389 typename std::enable_if<
390 (Tensor_Dim01 >= Dim),
391 Tensor1_Expr<Dg_numeral_12<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>, T,
392 Dim, i>>::type
393 operator()(const int N1, const Index<i, Dim> index, const int N2) const
394 {
395 using TensorExpr
396 = Dg_numeral_12<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>;
397 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
398 }
399
400 /* Index in third slot. */
401
402 template <char i, int Dim>
403 typename std::enable_if<
404 (Tensor_Dim2 >= Dim),
405 Tensor1_Expr<Dg_numeral_01<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>, T,
406 Dim, i>>::type
407 operator()(const int N1, const int N2, const Index<i, Dim>) const
408 {
409 using TensorExpr
410 = Dg_numeral_01<const Dg<T, Tensor_Dim01, Tensor_Dim2>, T>;
411 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
412 }
413 };
414}
N2
static Number< 2 > N2
Definition: BasicFeTools.hpp:101
N1
static Number< 1 > N1
Definition: BasicFeTools.hpp:100
N0
static Number< 0 > N0
Definition: BasicFeTools.hpp:99
FTensor::Dg_Expr
Definition: Dg_Expr.hpp:26
FTensor::Dg_number_01
Definition: Dg_number.hpp:54
FTensor::Dg_number_0
Definition: Dg_number.hpp:12
FTensor::Dg_number_12
Definition: Dg_number.hpp:40
FTensor::Dg_number_2
Definition: Dg_number.hpp:26
FTensor::Dg_number_rhs_01
Definition: Dg_number.hpp:63
FTensor::Dg_number_rhs_0
Definition: Dg_number.hpp:21
FTensor::Dg_number_rhs_12
Definition: Dg_number.hpp:49
FTensor::Dg_number_rhs_2
Definition: Dg_number.hpp:35
FTensor::Dg_numeral_01
Definition: Dg_numeral.hpp:50
FTensor::Dg_numeral_0
Definition: Dg_numeral.hpp:12
FTensor::Dg_numeral_12
Definition: Dg_numeral.hpp:36
FTensor::Dg_numeral_2
Definition: Dg_numeral.hpp:24
FTensor::Dg::operator()
T & operator()(const int N1, const int N2, const int N3)
Definition: Dg_value.hpp:26
FTensor::Dg::data
T data[(Tensor_Dim01 *(Tensor_Dim01+1))/2][Tensor_Dim2]
Definition: Dg_value.hpp:11
FTensor::Dg::operator()
T operator()(const int N1, const int N2, const int N3) const
Definition: Dg_value.hpp:43
FTensor::Dg::Dg
Dg(U... d)
Definition: Dg_value.hpp:14
FTensor::Tensor1_Expr
Definition: Tensor1_Expr.hpp:28
FTensor::Tensor2_Expr
Definition: Tensor2_Expr.hpp:27
FTensor::Tensor2_symmetric_Expr
Definition: Tensor2_symmetric_Expr.hpp:37
FTensor::Tensor3_contracted_01
Definition: Tensor3_contracted.hpp:49
FTensor::Tensor3_contracted_02
Definition: Tensor3_contracted.hpp:30
FTensor::Tensor3_contracted_12
Definition: Tensor3_contracted.hpp:11
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
T
const double T
Definition: mixed_reac_diff.cpp:158
FTensor
Tensors class implemented by Walter Landry.
Definition: FTensor.hpp:51
FTensor::d
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
N
const int N
Definition: speed_test.cpp:3
std::enable_if
Definition: enable_if.hpp:8
Generated by Doxygen 1.9.5 and hosted at University of Glasgow