PLaSK library
Loading...
Searching...
No Matches
tensor3.hpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#ifndef PLASK__TENSOR3_H
15#define PLASK__TENSOR3_H
16
21#include <iostream>
22
23#include "../math.hpp"
24#include "2d.hpp"
25#include "3d.hpp"
26#include "tensor2.hpp"
27
28namespace plask {
29
36template <typename T> struct Tensor3 {
37 // clang-format off
38 T c00,
39 c01,
40 c02,
41 c10,
42 c11,
43 c12,
44 c20,
45 c21,
46 c22;
47 // clang-format on
48
49 T& tran() { return c00; }
50 const T& lon() const { return c00; }
51
52 T& lon() { return c11; }
53 const T& tran() const { return c11; }
54
55 T& vert() { return c22; }
56 const T& vert() const { return c22; }
57
59 Tensor3() {}
60
65 template <typename OtherT>
66 Tensor3(const Tensor3<OtherT>& p)
67 : c00(p.c00), c01(p.c01), c02(p.c02), c10(p.c10), c11(p.c11), c12(p.c12), c20(p.c20), c21(p.c21), c22(p.c22) {}
68
73 Tensor3(const T& val) : c00(val), c01(0.), c02(0.), c10(0.), c11(val), c12(0.), c20(0.), c21(0.), c22(val) {}
74
79 Tensor3(const T& c00, const T& c22) : c00(c00), c01(0.), c02(0.), c10(0.), c11(c00), c12(0.), c20(0.), c21(0.), c22(c22) {}
80
85 Tensor3(const T& c00, const T& c11, const T& c22)
86 : c00(c00), c01(0.), c02(0.), c10(0.), c11(c11), c12(0.), c20(0.), c21(0.), c22(c22) {}
87
92 Tensor3(const T& c00, const T& c11, const T& c22, const T& c01)
93 : c00(c00), c01(c01), c02(0.), c10(conj(c01)), c11(c11), c12(0.), c20(0.), c21(0.), c22(c22) {}
94
99 Tensor3(const T& c00, const T& c11, const T& c22, const T& c01, const T& c10)
100 : c00(c00), c01(c01), c02(0.), c10(c10), c11(c11), c12(0.), c20(0.), c21(0.), c22(c22) {}
101
106 Tensor3(const T& c00, const T& c11, const T& c22, const T& c01, const T& c02, const T& c12)
107 : c00(c00), c01(c01), c02(c02), c10(conj(c01)), c11(c11), c12(c12), c20(conj(c02)), c21(conj(c12)), c22(c22) {}
108
113 Tensor3(const T& c00,
114 const T& c11,
115 const T& c22,
116 const T& c01,
117 const T& c10,
118 const T& c02,
119 const T& c20,
120 const T& c12,
121 const T& c21)
122 : c00(c00), c01(c01), c02(c02), c10(c10), c11(c11), c12(c12), c20(c20), c21(c21), c22(c22) {}
123
128 template <typename T0, typename T1>
129 Tensor3(const std::pair<T0, T1>& comp)
130 : c00(comp.first), c01(0.), c02(0.), c10(0.), c11(comp.first), c12(0.), c20(0.), c21(0.), c22(comp.second) {}
131
136 Tensor3(const Tensor2<T>& tens)
137 : c00(tens.c00), c01(0.), c02(0.), c10(0.), c11(tens.c00), c12(0.), c20(0.), c21(0.), c22(tens.c11) {}
138
143 Tensor3(const Vec<2, T>& vec) : c00(vec.c0), c01(0.), c02(0.), c10(0.), c11(vec.c0), c12(0.), c20(0.), c21(0.), c22(vec.c1) {}
144
149 Tensor3(const Vec<3, T>& vec) : c00(vec.c0), c01(0.), c02(0.), c10(0.), c11(vec.c1), c12(0.), c20(0.), c21(0.), c22(vec.c2) {}
150
155 Tensor3(const T* data)
156 : c00(data[0]),
157 c01(data[1]),
158 c02(data[2]),
159 c10(data[3]),
160 c11(data[4]),
161 c12(data[5]),
162 c20(data[6]),
163 c21(data[7]),
164 c22(data[8]) {}
165
172 inline T& operator[](size_t i) {
173 assert(i < 9);
174 return *(&c00 + i);
175 }
176
183 inline const T& operator[](size_t i) const {
184 assert(i < 9);
185 return *(&c00 + i);
186 }
187
194 inline T& operator()(size_t i, size_t j) {
195 assert(i < 3 && j < 3);
196 return *(&c00 + 3 * i + j);
197 }
198
205 inline const T& operator()(size_t i, size_t j) const {
206 assert(i < 3 && j < 3);
207 return *(&c00 + 3 * i + j);
208 }
209
211 operator std::tuple<T, T, T, T, T, T, T, T, T>() const { return std::make_tuple(c00, c11, c22, c01, c10, c02, c20, c12, c21); }
212
218 template <typename OtherT> bool operator==(const Tensor3<OtherT>& p) const {
219 return p.c00 == c00 && p.c01 == c01 && p.c02 == c02 && p.c10 == c10 && p.c11 == c11 && p.c12 == c12 && p.c20 == c20 &&
220 p.c21 == c21 && p.c22 == c22;
221 }
222
228 template <typename OtherT> constexpr bool equals(const Tensor3<OtherT>& p) const {
229 return is_zero(p.c00 - c00) && is_zero(p.c01 - c01) && is_zero(p.c02 - c02) && is_zero(p.c10 - c10) &&
230 is_zero(p.c11 - c11) && is_zero(p.c12 - c12) && is_zero(p.c20 - c20) && is_zero(p.c21 - c21) && is_zero(p.c22 - c22);
231 }
232
238 template <typename OtherT> bool operator!=(const Tensor3<OtherT>& p) const {
239 return p.c00 != c00 || p.c01 != c01 || p.c02 != c02 || p.c10 != c10 || p.c11 != c11 || p.c12 != c12 || p.c20 != c20 ||
240 p.c21 != c21 || p.c22 != c22;
241 }
242
248 template <typename OtherT> auto operator+(const Tensor3<OtherT>& other) const -> Tensor3<decltype(c00 + other.c00)> {
249 return Tensor3<decltype(this->c00 + other.c00)>(c00 + other.c00, c11 + other.c11, c22 + other.c22, c01 + other.c01,
250 c10 + other.c10, c02 + other.c02, c20 + other.c20, c12 + other.c12,
251 c21 + other.c21);
252 }
253
259 Tensor3<T>& operator+=(const Tensor3<T>& other) {
260 c00 += other.c00;
261 c01 += other.c01;
262 c02 += other.c02;
263 c10 += other.c10;
264 c11 += other.c11;
265 c12 += other.c12;
266 c20 += other.c20;
267 c21 += other.c21;
268 c22 += other.c22;
269 return *this;
270 }
271
278 template <typename OtherT> auto operator-(const Tensor3<OtherT>& other) const -> Tensor3<decltype(c00 - other.c00)> {
279 return Tensor3<decltype(this->c00 - other.c00)>(c00 - other.c00, c11 - other.c11, c22 - other.c22, c01 - other.c01,
280 c10 - other.c10, c02 - other.c02, c20 - other.c20, c12 - other.c12,
281 c21 - other.c21);
282 }
283
289 Tensor3<T>& operator-=(const Tensor3<T>& other) {
290 c00 -= other.c00;
291 c01 -= other.c01;
292 c02 -= other.c02;
293 c10 -= other.c10;
294 c11 -= other.c11;
295 c12 -= other.c12;
296 c20 -= other.c20;
297 c21 -= other.c21;
298 c22 -= other.c22;
299 return *this;
300 }
301
307 template <typename OtherT> auto operator*(const OtherT scale) const -> Tensor3<decltype(c00 * scale)> {
308 return Tensor3<decltype(c00 * scale)>(c00 * scale, c11 * scale, c22 * scale, c01 * scale, c10 * scale, c02 * scale,
309 c20 * scale, c12 * scale, c21 * scale);
310 }
311
317 Tensor3<T>& operator*=(const T scalar) {
318 c00 *= scalar;
319 c01 *= scalar;
320 c02 *= scalar;
321 c10 *= scalar;
322 c11 *= scalar;
323 c12 *= scalar;
324 c20 *= scalar;
325 c21 *= scalar;
326 c22 *= scalar;
327 return *this;
328 }
329
335 Tensor3<T> operator/(const T scale) const {
336 return Tensor3<decltype(c00 / scale)>(c00 / scale, c11 / scale, c22 / scale, c01 / scale, c10 / scale, c02 / scale,
337 c20 / scale, c12 / scale, c21 / scale);
338 }
339
345 Tensor3<T>& operator/=(const T scalar) {
346 c00 /= scalar;
347 c01 /= scalar;
348 c02 /= scalar;
349 c10 /= scalar;
350 c11 /= scalar;
351 c12 /= scalar;
352 c20 /= scalar;
353 c21 /= scalar;
354 c22 /= scalar;
355 return *this;
356 }
357
362 Tensor3<T> operator-() const { return Tensor3<T>(-c00, -c11, -c22, -c01, -c10, -c02, -c20, -c12, -c21); }
363
368 Tensor3<T> sqr() const {
369 return Tensor3<T>(c00 * c00 + c01 * c10 + c02 * c20, // c00
370 c10 * c01 + c11 * c11 + c12 * c21, // c11
371 c20 * c02 + c21 * c12 + c22 * c22, // c22
372 c00 * c01 + c01 * c11 + c02 * c21, // c01
373 c10 * c00 + c11 * c10 + c12 * c20, // c10
374 c00 * c02 + c01 * c12 + c02 * c22, // c02
375 c20 * c00 + c21 * c10 + c22 * c20, // c20
376 c10 * c02 + c11 * c12 + c12 * c22, // c12
377 c20 * c01 + c21 * c11 + c22 * c21 // c21
378 );
379 }
380
385 Tensor3<T> pow(int n) const {
386 if (n < 0) {
387 return inv().pow(-n);
388 } else if (n == 0) {
389 return Tensor3<T>(1.);
390 } else if (n == 1) {
391 return *this;
392 } else if (n == 2) {
393 return sqr();
394 } else if (n % 2 == 0) {
395 return sqr().pow(n / 2);
396 } else {
397 Tensor3<T> a = sqr().pow(n / 2), b(*this);
398 return Tensor3<T>(a.c00 * b.c00 + a.c01 * b.c10 + a.c02 * b.c20, // c00
399 a.c10 * b.c01 + a.c11 * b.c11 + a.c12 * b.c21, // c11
400 a.c20 * b.c02 + a.c21 * b.c12 + a.c22 * b.c22, // c22
401 a.c00 * b.c01 + a.c01 * b.c11 + a.c02 * b.c21, // c01
402 a.c10 * b.c00 + a.c11 * b.c10 + a.c12 * b.c20, // c10
403 a.c00 * b.c02 + a.c01 * b.c12 + a.c02 * b.c22, // c02
404 a.c20 * b.c00 + a.c21 * b.c10 + a.c22 * b.c20, // c20
405 a.c10 * b.c02 + a.c11 * b.c12 + a.c12 * b.c22, // c12
406 a.c20 * b.c01 + a.c21 * b.c11 + a.c22 * b.c21 // c21
407 );
408 }
409 }
410
411 // /**
412 // * Square root of each component of tensor
413 // * \return squared tensor
414 // */
415 // Tensor3<T> sqrt() const {
416 // TODO
417 // }
418
424 Tensor3<T> inv() const {
425 // clang-format off
426 T a00 = c11*c22 - c12*c21, a01 = c02*c21 - c01*c22, a02 = c01*c12 - c02*c11,
427 a10 = c12*c20 - c10*c22, a11 = c00*c22 - c02*c20, a12 = c02*c10 - c00*c12,
428 a20 = c10*c21 - c11*c20, a21 = c01*c20 - c00*c21, a22 = c00*c11 - c01*c10;
429 // clang-format on
430
431 T det = c00 * c11 * c22 + c01 * c12 * c20 + c02 * c10 * c21 - c02 * c11 * c20 - c00 * c12 * c21 - c10 * c01 * c22;
432 return Tensor3<T>(a00 / det, a11 / det, a22 / det, a01 / det, a10 / det, a02 / det, a20 / det, a12 / det, a21 / det);
433 }
434
441 friend inline std::ostream& operator<<(std::ostream& out, const Tensor3<T>& to_print) {
442 return out << '[' << str(to_print.c00) << ", " << str(to_print.c01) << ", " << str(to_print.c02) << "; "
443 << str(to_print.c10) << ", " << str(to_print.c11) << ", " << str(to_print.c12) << "; " << str(to_print.c20)
444 << ", " << str(to_print.c21) << ", " << str(to_print.c22) << "]";
445 }
446};
447
454template <typename T, typename OtherT> auto operator*(const OtherT scale, const Tensor3<T>& tensor) -> decltype(tensor * scale) {
455 return tensor * scale;
456}
457
463template <typename T> inline Tensor3<T> conj(const Tensor3<T>& v) {
464 return Tensor3<T>(conj(v.c00), conj(v.c11), conj(v.c22), conj(v.c01), conj(v.c10), conj(v.c02), conj(v.c20), conj(v.c12),
465 conj(v.c21));
466}
467
469template <typename T> struct NaNImpl<Tensor3<T>> {
470 static constexpr Tensor3<T> get() { return Tensor3<T>(NaN<T>()); }
471};
472
474template <typename T> struct ZeroImpl<Tensor3<T>> {
475 static constexpr Tensor3<T> get() { return Tensor3<T>(0.); }
476};
477
480template <typename T> inline bool is_zero(const Tensor3<T>& v) {
481 return is_zero(v.c00) && is_zero(v.c01) && is_zero(v.c02) && is_zero(v.c10) && is_zero(v.c11) && is_zero(v.c12) &&
482 is_zero(v.c20) && is_zero(v.c21) && is_zero(v.c22);
483}
484
485/*
486PLASK_API_EXTERN_TEMPLATE_STRUCT(Tensor3<double>)
487PLASK_API_EXTERN_TEMPLATE_STRUCT(Tensor3< std::complex<double> >)
488*/
489
490template <typename T>
491inline bool isnan(plask::Tensor3<T> tens) {
492 return isnan(tens.c00) || isnan(tens.c01) || isnan(tens.c02) || isnan(tens.c10) || isnan(tens.c11) || isnan(tens.c12) ||
493 isnan(tens.c20) || isnan(tens.c21) || isnan(tens.c22);
494}
495
496} // namespace plask
497
498namespace std {
499
500template <typename T> plask::Tensor3<T> pow(plask::Tensor3<T> tens, int n) { return tens.pow(n); }
501
502} // namespace std
503
504#endif // PLASK__TENSOR3_H