MiniTensor Version of the Day
Loading...
Searching...
No Matches
MiniTensor_Tensor.h
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MiniTensor Package
4//
5// Copyright 2016 NTESS and the MiniTensor contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#if !defined(MiniTensor_Tensor_h)
11#define MiniTensor_Tensor_h
12
13#include <algorithm>
14#include <cassert>
15#include <iomanip>
16#include <iostream>
17#include <vector>
18
19#include <Kokkos_Core.hpp>
20#include "MiniTensor_Vector.h"
21
22namespace minitensor {
23
26
30template<typename T, Index N>
32
36template<typename T, Index N = DYNAMIC>
37class Tensor: public TensorBase<T, tensor_store<T, N>>
38{
39public:
40
44 static constexpr
45 Index
46 ORDER = 2;
47
51 static constexpr
52 bool
54
59
64 static constexpr
65 Index
67 {
68 return ORDER;
69 }
70
74 explicit
77
82 explicit
84 Tensor(Index const dimension);
85
90 explicit
92 Tensor(Filler const value);
93
99 explicit
101 Tensor(Index const dimension, Filler const value);
102
109 template<class ArrayT>
110 explicit
113 Source const source,
114 ArrayT & data,
115 Index index1);
116
124 template<class ArrayT>
125 explicit
128 Source const source,
129 ArrayT & data,
130 Index index1,
131 Index index2);
132
141 template<class ArrayT>
142 explicit
145 Source const source,
146 ArrayT & data,
147 Index index1,
148 Index index2,
149 Index index3);
150
160 template<class ArrayT>
161 explicit
164 Source const source,
165 ArrayT & data,
166 Index index1,
167 Index index2,
168 Index index3,
169 Index index4);
170
181 template<class ArrayT>
182 explicit
185 Source const source,
186 ArrayT & data,
187 Index index1,
188 Index index2,
189 Index index3,
190 Index index4,
191 Index index5);
192
204 template<class ArrayT>
205 explicit
208 Source const source,
209 ArrayT & data,
210 Index index1,
211 Index index2,
212 Index index3,
213 Index index4,
214 Index index5,
215 Index index6);
216
224 template<class ArrayT>
225 explicit
228 Source const source,
229 Index const dimension,
230 ArrayT & data,
231 Index index1);
232
241 template<class ArrayT>
242 explicit
245 Source const source,
246 Index const dimension,
247 ArrayT & data,
248 Index index1,
249 Index index2);
250
260 template<class ArrayT>
261 explicit
264 Source const source,
265 Index const dimension,
266 ArrayT & data,
267 Index index1,
268 Index index2,
269 Index index3);
270
281 template<class ArrayT>
282 explicit
285 Source const source,
286 Index const dimension,
287 ArrayT & data,
288 Index index1,
289 Index index2,
290 Index index3,
291 Index index4);
292
304 template<class ArrayT>
305 explicit
308 Source const source,
309 Index const dimension,
310 ArrayT & data,
311 Index index1,
312 Index index2,
313 Index index3,
314 Index index4,
315 Index index5);
316
329 template<class ArrayT>
330 explicit
333 Source const source,
334 Index const dimension,
335 ArrayT & data,
336 Index index1,
337 Index index2,
338 Index index3,
339 Index index4,
340 Index index5,
341 Index index6);
342
347 explicit
349 Tensor(T const * data_ptr);
350
356 explicit
358 Tensor(Index const dimension, T const * data_ptr);
359
365
371
379 //
380 explicit
382 Tensor(T const & s00, T const & s01, T const & s10, T const & s11);
383
396 explicit
399 T const & s00, T const & s01, T const & s02,
400 T const & s10, T const & s11, T const & s12,
401 T const & s20, T const & s21, T const & s22);
402
408 explicit
410 Tensor(T const * data_ptr, ComponentOrder const component_order);
411
418 explicit
421 Index const dimension,
422 T const * data_ptr,
423 ComponentOrder const component_order);
424
428 virtual
431
438 T const &
439 operator()(Index const i, Index const j) const;
440
447 T &
448 operator()(Index const i, Index const j);
449
454 Index
456
461 Index
463
468 Index
470
475 void
476 set_dimension(Index const dimension);
477};
478
483template<typename S, typename T, Index N>
486operator+(Tensor<S, N> const & A, Tensor<T, N> const & B);
487
492template<typename S, typename T, Index N>
495operator-(Tensor<S, N> const & A, Tensor<T, N> const & B);
496
501template<typename T, Index N>
504operator-(Tensor<T, N> const & A);
505
511template<typename T, Index N>
513bool
514operator==(Tensor<T, N> const & A, Tensor<T, N> const & B);
515
521template<typename T, Index N>
523bool
524operator!=(Tensor<T, N> const & A, Tensor<T, N> const & B);
525
532template<typename S, typename T, Index N>
535operator*(Tensor<T, N> const & A, Vector<S, N> const & u);
536
543template<typename S, typename T, Index N>
546operator*(Vector<S, N> const & u, Tensor<T, N> const & A);
547
552template<typename S, typename T, Index N>
555operator*(Tensor<S, N> const & A, Tensor<T, N> const & B);
556
563template<typename S, typename T, Index N>
565typename lazy_disable_if<order_1234<S>, apply_tensor<Promote<S, T>, N>>::type
566operator*(S const & s, Tensor<T, N> const & A);
567
574template<typename S, typename T, Index N>
576typename lazy_disable_if<order_1234<S>, apply_tensor<Promote<S, T>, N>>::type
577operator*(Tensor<T, N> const & A, S const & s);
578
585template<typename S, typename T, Index N>
588operator/(Tensor<T, N> const & A, S const & s);
589
596template<typename S, typename T, Index N>
599operator/(S const & s, Tensor<T, N> const & A);
600
607template<typename T, Index N>
608std::istream &
609operator>>(std::istream & is, Tensor<T, N> & A);
610
617template<typename T, Index N>
618std::ostream &
619operator<<(std::ostream & os, Tensor<T, N> const & A);
620
627template<typename T, Index N>
630row(Tensor<T, N> const & A, Index const i);
631
638template<typename T, Index N>
641col(Tensor<T, N> const & A, Index const j);
642
649template<typename S, typename T, Index N>
652dot(Tensor<T, N> const & A, Vector<S, N> const & u);
653
660template<typename S, typename T, Index N>
663dot(Vector<S, N> const & u, Tensor<T, N> const & A);
664
671template<typename S, typename T, Index N>
674dot(Tensor<S, N> const & A, Tensor<T, N> const & B);
675
682template<typename S, typename T, Index N>
685t_dot(Tensor<S, N> const & A, Tensor<T, N> const & B);
686
693template<typename S, typename T, Index N, Index P>
696t_dot(Tensor<S, N> const & A, Matrix<T, N, P> const & B);
697
704template<typename S, typename T, Index N>
707t_dot(Tensor<S, N> const & A, Vector<T, N> const & B);
708
715template<typename S, typename T, Index N>
718dot_t(Tensor<S, N> const & A, Tensor<T, N> const & B);
719
726template<typename S, typename T, Index N>
729t_dot_t(Tensor<S, N> const & A, Tensor<T, N> const & B);
730
737template<typename S, typename T, Index N>
739typename Promote<S, T>::type
740dotdot(Tensor<S, N> const & A, Tensor<T, N> const & B);
741
748template<typename S, typename T, Index N>
751dyad(Vector<S, N> const & u, Vector<T, N> const & v);
752
759template<typename S, typename T, Index N>
762bun(Vector<S, N> const & u, Vector<T, N> const & v);
763
770template<typename S, typename T, Index N>
773tensor(Vector<S, N> const & u, Vector<T, N> const & v);
774
780template<typename T, Index N>
783diag(Vector<T, N> const & v);
784
790template<typename T, Index N>
793diag(Tensor<T, N> const & A);
794
799template<typename T, Index N>
801Tensor<T, N> const
802zero();
803
804template<typename T>
807zero(Index const dimension);
808
809template<typename T, Index N>
811Tensor<T, N> const
812zero(Index const dimension);
813
817template<typename T, Index N>
819Tensor<T, N> const
820identity();
821
822template<typename T>
825identity(Index const dimension);
826
827template<typename T, Index N>
829Tensor<T, N> const
830identity(Index const dimension);
831
835template<typename T, Index N>
837Tensor<T, N> const
838eye();
839
840template<typename T>
843eye(Index const dimension);
844
845template<typename T, Index N>
847Tensor<T, N> const
848eye(Index const dimension);
849
853template<typename T, Index N>
855Tensor<T, N> const
857
858template<typename T>
861levi_civita_2(Index const dimension);
862
863template<typename T, Index N>
865Tensor<T, N> const
866levi_civita_2(Index const dimension);
867
871template<typename T, Index N>
873Tensor<T, N> const
875
876template<typename T>
879permutation_2(Index const dimension);
880
881template<typename T, Index N>
883Tensor<T, N> const
884permutation_2(Index const dimension);
885
889template<typename T, Index N>
891Tensor<T, N> const
893
894template<typename T>
897alternator_2(Index const dimension);
898
899template<typename T, Index N>
901Tensor<T, N> const
902alternator_2(Index const dimension);
903
907template<typename T, Index N>
910transpose(Tensor<T, N> const & A);
911
915template<typename T, Index N>
919
924template<typename T, Index N>
927sym(Tensor<T, N> const & A);
928
933template<typename T, Index N>
936skew(Tensor<T, N> const & A);
937
944template<typename T, Index N>
947skew(Vector<T, N> const & u);
948
949} // namespace minitensor
950
951namespace minitensor {
952
953//
954// Constructor that initializes to NaNs
955//
956template<typename T, Index N>
960{
961 set_dimension(N);
962 return;
963}
964
965template<typename T, Index N>
967Tensor<T, N>::Tensor(Index const dimension) :
968 TensorBase<T, Store>::TensorBase(dimension, ORDER)
969{
970 return;
971}
972
973//
974// Create tensor from a specified value
975//
976template<typename T, Index N>
979 TensorBase<T, Store>::TensorBase(N, ORDER, value)
980{
981 return;
982}
983
984template<typename T, Index N>
986Tensor<T, N>::Tensor(Index const dimension, Filler const value) :
987 TensorBase<T, Store>::TensorBase(dimension, ORDER, value)
988{
989 return;
990}
991
992//
993// Create tensor from array
994//
995//
996template<typename T, Index N>
997template<class ArrayT>
1000 Source const,
1001 ArrayT & data,
1002 Index index1) :
1003 TensorBase<T, Store>::TensorBase(N, ORDER, data, index1)
1004{
1005 return;
1006}
1007
1008template<typename T, Index N>
1009template<class ArrayT>
1012 Source const,
1013 ArrayT & data,
1014 Index index1,
1015 Index index2) :
1016 TensorBase<T, Store>::TensorBase(N, ORDER, data, index1, index2)
1017{
1018 return;
1019}
1020
1021template<typename T, Index N>
1022template<class ArrayT>
1025 Source const,
1026 ArrayT & data,
1027 Index index1,
1028 Index index2,
1029 Index index3) :
1030 TensorBase<T, Store>::TensorBase(N, ORDER, data, index1, index2, index3)
1031{
1032 return;
1033}
1034
1035template<typename T, Index N>
1036template<class ArrayT>
1039 Source const,
1040 ArrayT & data,
1041 Index index1,
1042 Index index2,
1043 Index index3,
1044 Index index4) :
1046 N,
1047 ORDER,
1048 data,
1049 index1,
1050 index2,
1051 index3,
1052 index4)
1053{
1054 return;
1055}
1056
1057template<typename T, Index N>
1058template<class ArrayT>
1061 Source const,
1062 ArrayT & data,
1063 Index index1,
1064 Index index2,
1065 Index index3,
1066 Index index4,
1067 Index index5) :
1069 N,
1070 ORDER,
1071 data,
1072 index1,
1073 index2,
1074 index3,
1075 index4,
1076 index5)
1077{
1078 return;
1079}
1080
1081template<typename T, Index N>
1082template<class ArrayT>
1085 Source const,
1086 ArrayT & data,
1087 Index index1,
1088 Index index2,
1089 Index index3,
1090 Index index4,
1091 Index index5,
1092 Index index6) :
1094 N,
1095 ORDER,
1096 data,
1097 index1,
1098 index2,
1099 index3,
1100 index4,
1101 index5,
1102 index6)
1103{
1104 return;
1105}
1106
1107template<typename T, Index N>
1108template<class ArrayT>
1111 Source const,
1112 Index const dimension,
1113 ArrayT & data,
1114 Index index1) :
1115 TensorBase<T, Store>::TensorBase(dimension, ORDER, data, index1)
1116{
1117 return;
1118}
1119
1120template<typename T, Index N>
1121template<class ArrayT>
1124 Source const,
1125 Index const dimension,
1126 ArrayT & data,
1127 Index index1,
1128 Index index2) :
1129 TensorBase<T, Store>::TensorBase(dimension, ORDER, data, index1, index2)
1130{
1131 return;
1132}
1133
1134template<typename T, Index N>
1135template<class ArrayT>
1138 Source const,
1139 Index const dimension,
1140 ArrayT & data,
1141 Index index1,
1142 Index index2,
1143 Index index3) :
1145 dimension,
1146 ORDER,
1147 data,
1148 index1,
1149 index2,
1150 index3)
1151{
1152 return;
1153}
1154
1155template<typename T, Index N>
1156template<class ArrayT>
1159 Source const,
1160 Index const dimension,
1161 ArrayT & data,
1162 Index index1,
1163 Index index2,
1164 Index index3,
1165 Index index4) :
1167 dimension,
1168 ORDER,
1169 data,
1170 index1,
1171 index2,
1172 index3,
1173 index4)
1174{
1175 return;
1176}
1177
1178template<typename T, Index N>
1179template<class ArrayT>
1182 Source const,
1183 Index const dimension,
1184 ArrayT & data,
1185 Index index1,
1186 Index index2,
1187 Index index3,
1188 Index index4,
1189 Index index5) :
1191 dimension,
1192 ORDER,
1193 data,
1194 index1,
1195 index2,
1196 index3,
1197 index4,
1198 index5)
1199{
1200 return;
1201}
1202
1203template<typename T, Index N>
1204template<class ArrayT>
1207 Source const,
1208 Index const dimension,
1209 ArrayT & data,
1210 Index index1,
1211 Index index2,
1212 Index index3,
1213 Index index4,
1214 Index index5,
1215 Index index6) :
1217 dimension,
1218 ORDER,
1219 data,
1220 index1,
1221 index2,
1222 index3,
1223 index4,
1224 index5,
1225 index6)
1226{
1227 return;
1228}
1229
1230template<typename T, Index N>
1232Tensor<T, N>::Tensor(T const * data_ptr) :
1233 TensorBase<T, Store>::TensorBase(N, ORDER, data_ptr)
1234{
1235 return;
1236}
1237
1238template<typename T, Index N>
1240Tensor<T, N>::Tensor(Index const dimension, T const * data_ptr) :
1241 TensorBase<T, Store>::TensorBase(dimension, ORDER, data_ptr)
1242{
1243 return;
1244}
1245//
1246// Copy constructor
1247//
1248template<typename T, Index N>
1252{
1253 return;
1254}
1255
1256//
1257// Create tensor specifying components
1258// \param s00, s01, ... components in the R^2 canonical basis
1259//
1260template<typename T, Index N>
1263 T const & s00, T const & s01,
1264 T const & s10, T const & s11)
1265{
1266 Tensor<T, N> &
1267 self = (*this);
1268
1269 self.set_dimension(2);
1270
1271 self[0] = s00;
1272 self[1] = s01;
1273
1274 self[2] = s10;
1275 self[3] = s11;
1276
1277 return;
1278}
1279
1280//
1281// Create tensor specifying components
1282// \param s00, s01, ... components in the R^3 canonical basis
1283//
1284template<typename T, Index N>
1287 T const & s00, T const & s01, T const & s02,
1288 T const & s10, T const & s11, T const & s12,
1289 T const & s20, T const & s21, T const & s22)
1290{
1291 Tensor<T, N> &
1292 self = (*this);
1293
1294 self.set_dimension(3);
1295
1296 self[0] = s00;
1297 self[1] = s01;
1298 self[2] = s02;
1299
1300 self[3] = s10;
1301 self[4] = s11;
1302 self[5] = s12;
1303
1304 self[6] = s20;
1305 self[7] = s21;
1306 self[8] = s22;
1307
1308 return;
1309}
1310
1311//
1312// Create tensor from array with component order
1313//
1314template<typename T, Index N>
1317 T const * data_ptr,
1318 ComponentOrder const component_order)
1319{
1320 assert(data_ptr != NULL);
1321
1322 fill(data_ptr, component_order);
1323
1324 return;
1325}
1326
1327template<typename T, Index N>
1330 Index const dimension,
1331 T const * data_ptr,
1332 ComponentOrder const component_order)
1333{
1334 assert(data_ptr != NULL);
1335
1336 Tensor<T, N> &
1337 self = (*this);
1338
1339 self.set_dimension(dimension);
1340
1341 fill(data_ptr, component_order);
1342
1343 return;
1344}
1345
1346//
1347// 2nd-order tensor from 4th-order tensor
1348//
1349template<typename T, Index N>
1352{
1353 Index const
1354 dimension_4th = A.get_dimension();
1355
1356 Index const
1357 dimension_2nd = dimension_4th * dimension_4th;
1358
1359 Tensor<T, N> &
1360 self = (*this);
1361
1362 self.set_dimension(dimension_2nd);
1363
1364 Index const
1365 number_components = dimension_2nd * dimension_2nd;
1366
1367 for (Index i = 0; i < number_components; ++i) {
1368 self[i] = A[i];
1369 }
1370
1371 return;
1372}
1373
1374//
1375// Simple destructor
1376//
1377template<typename T, Index N>
1380{
1381 return;
1382}
1383
1384//
1385// Get dimension
1386//
1387template<typename T, Index N>
1389Index
1394
1395//
1396// Get number rows
1397//
1398template<typename T, Index N>
1400Index
1402{
1403 return get_dimension();
1404}
1405
1406//
1407// Get number columns
1408//
1409template<typename T, Index N>
1411Index
1413{
1414 return get_dimension();
1415}
1416
1417//
1418// Set dimension
1419//
1420template<typename T, Index N>
1422void
1424{
1425 TensorBase<T, Store>::set_dimension(dimension, ORDER);
1426 return;
1427}
1428
1429//
1430// Indexing for constant tensor
1431//
1432template<typename T, Index N>
1434T const &
1435Tensor<T, N>::operator()(Index const i, Index const j) const
1436 {
1437 Tensor<T, N> const &
1438 self = (*this);
1439
1440 Index const
1441 dimension = self.get_dimension();
1442
1443 return self[i * dimension + j];
1444}
1445
1446//
1447//Tensor indexing
1448//
1449template<typename T, Index N>
1451T &
1453{
1454 Tensor<T, N> &
1455 self = (*this);
1456
1457 Index const
1458 dimension = self.get_dimension();
1459
1460 return self[i * dimension + j];
1461}
1462
1463namespace {
1464
1465template<typename S>
1467bool
1468greater_than(S const & a, S const & b)
1469{
1470 return a.first > b.first;
1471}
1472
1473} // anonymous namespace
1474
1475//
1476// Extract a row as a vector
1477//
1478template<typename T, Index N>
1480Vector<T, N>
1481row(Tensor<T, N> const & A, Index const i)
1482{
1483 Index const
1484 dimension = A.get_dimension();
1485
1487 v(dimension);
1488
1489 switch (dimension) {
1490 default:
1491 for (Index j = 0; j < dimension; ++j) {
1492 v(j) = A(i, j);
1493 }
1494 break;
1495
1496 case 2:
1497 v(0) = A(i, 0);
1498 v(1) = A(i, 1);
1499 break;
1500
1501 case 3:
1502 v(0) = A(i, 0);
1503 v(1) = A(i, 1);
1504 v(2) = A(i, 2);
1505 break;
1506 }
1507
1508 return v;
1509}
1510
1511//
1512// Extract a column as a vector
1513//
1514template<typename T, Index N>
1516Vector<T, N>
1517col(Tensor<T, N> const & A, Index const j)
1518{
1519 Index const
1520 dimension = A.get_dimension();
1521
1523 v(dimension);
1524
1525 switch (dimension) {
1526 default:
1527 for (Index i = 0; i < dimension; ++i) {
1528 v(i) = A(i, j);
1529 }
1530 break;
1531
1532 case 2:
1533 v(0) = A(0, j);
1534 v(1) = A(1, j);
1535 break;
1536
1537 case 3:
1538 v(0) = A(0, j);
1539 v(1) = A(1, j);
1540 v(2) = A(2, j);
1541 break;
1542 }
1543
1544 return v;
1545}
1546
1547//
1548// Tensor addition
1549//
1550template<typename S, typename T, Index N>
1552Tensor<typename Promote<S, T>::type, N>
1554{
1556 C(A.get_dimension());
1557
1558 add(A, B, C);
1559
1560 return C;
1561}
1562
1563//
1564// Tensor subtraction
1565//
1566template<typename S, typename T, Index N>
1568Tensor<typename Promote<S, T>::type, N>
1570{
1572 C(A.get_dimension());
1573
1574 subtract(A, B, C);
1575
1576 return C;
1577}
1578
1579//
1580// Tensor minus
1581//
1582template<typename T, Index N>
1584Tensor<T, N>
1586{
1588 B(A.get_dimension());
1589
1590 minus(A, B);
1591
1592 return B;
1593}
1594
1595//
1596// Tensor equality
1597//
1598template<typename T, Index N>
1600bool
1602{
1603 return equal(A, B);
1604}
1605
1606//
1607// Tensor inequality
1608//
1609template<typename T, Index N>
1611bool
1613{
1614 return not_equal(A, B);
1615}
1616
1617//
1618// Scalar tensor product
1619//
1620template<typename S, typename T, Index N>
1622typename
1623lazy_disable_if<order_1234<S>, apply_tensor<Promote<S, T>, N>>::type
1624operator*(S const & s, Tensor<T, N> const & A)
1625{
1627 B(A.get_dimension());
1628
1629 scale(A, s, B);
1630
1631 return B;
1632}
1633
1634//
1635// Tensor scalar product
1636//
1637template<typename S, typename T, Index N>
1639typename
1640lazy_disable_if<order_1234<S>, apply_tensor<Promote<S, T>, N>>::type
1641operator*(Tensor<T, N> const & A, S const & s)
1642{
1644 B(A.get_dimension());
1645
1646 scale(A, s, B);
1647
1648 return B;
1649}
1650
1651//
1652// Tensor scalar division
1653//
1654template<typename S, typename T, Index N>
1656Tensor<typename Promote<S, T>::type, N>
1657operator/(Tensor<T, N> const & A, S const & s)
1658{
1660 B(A.get_dimension());
1661
1662 divide(A, s, B);
1663
1664 return B;
1665}
1666
1667//
1668// Scalar tensor division
1669//
1670template<typename S, typename T, Index N>
1672Tensor<typename Promote<S, T>::type, N>
1673operator/(S const & s, Tensor<T, N> const & A)
1674{
1676 B(A.get_dimension());
1677
1678 split(A, s, B);
1679
1680 return B;
1681}
1682
1683//
1684// Tensor vector product v = A u
1685//
1686template<typename S, typename T, Index N>
1688Vector<typename Promote<S, T>::type, N>
1690{
1691 return dot(A, u);
1692}
1693
1694//
1695// Vector tensor product v = u A
1696//
1697template<typename S, typename T, Index N>
1699Vector<typename Promote<S, T>::type, N>
1701{
1702 return dot(u, A);
1703}
1704
1705//
1706// Tensor dot product C = A B
1707//
1708template<typename S, typename T, Index N>
1710Tensor<typename Promote<S, T>::type, N>
1712{
1713 return dot(A, B);
1714}
1715
1716//
1717// Tensor vector product v = A u
1718//
1719template<typename S, typename T, Index N>
1721Vector<typename Promote<S, T>::type, N>
1722dot(Tensor<T, N> const & A, Vector<S, N> const & u)
1723{
1724 Index const
1725 dimension = A.get_dimension();
1726
1727 assert(u.get_dimension() == dimension);
1728
1730 v(dimension);
1731
1732 switch (dimension) {
1733
1734 default:
1735 for (Index i = 0; i < dimension; ++i) {
1736
1737 typename Promote<S, T>::type
1738 s = 0.0;
1739
1740 for (Index p = 0; p < dimension; ++p) {
1741 s += A(i, p) * u(p);
1742 }
1743 v(i) = s;
1744 }
1745 break;
1746
1747 case 3:
1748 v(0) = A(0, 0) * u(0) + A(0, 1) * u(1) + A(0, 2) * u(2);
1749 v(1) = A(1, 0) * u(0) + A(1, 1) * u(1) + A(1, 2) * u(2);
1750 v(2) = A(2, 0) * u(0) + A(2, 1) * u(1) + A(2, 2) * u(2);
1751 break;
1752
1753 case 2:
1754 v(0) = A(0, 0) * u(0) + A(0, 1) * u(1);
1755 v(1) = A(1, 0) * u(0) + A(1, 1) * u(1);
1756 break;
1757
1758 }
1759
1760 return v;
1761}
1762
1763//
1764// Vector tensor product v = u A
1765//
1766template<typename S, typename T, Index N>
1768Vector<typename Promote<S, T>::type, N>
1769dot(Vector<S, N> const & u, Tensor<T, N> const & A)
1770{
1771 Index const
1772 dimension = A.get_dimension();
1773
1774 assert(u.get_dimension() == dimension);
1775
1777 v(dimension);
1778
1779 switch (dimension) {
1780
1781 default:
1782 for (Index i = 0; i < dimension; ++i) {
1783
1784 typename Promote<S, T>::type
1785 s = 0.0;
1786
1787 for (Index p = 0; p < dimension; ++p) {
1788 s += A(p, i) * u(p);
1789 }
1790 v(i) = s;
1791 }
1792 break;
1793
1794 case 3:
1795 v(0) = A(0, 0) * u(0) + A(1, 0) * u(1) + A(2, 0) * u(2);
1796 v(1) = A(0, 1) * u(0) + A(1, 1) * u(1) + A(2, 1) * u(2);
1797 v(2) = A(0, 2) * u(0) + A(1, 2) * u(1) + A(2, 2) * u(2);
1798 break;
1799
1800 case 2:
1801 v(0) = A(0, 0) * u(0) + A(1, 0) * u(1);
1802 v(1) = A(0, 1) * u(0) + A(1, 1) * u(1);
1803 break;
1804
1805 }
1806
1807 return v;
1808}
1809
1810//
1811// Tensor tensor product C = A B
1812//
1813template<typename S, typename T, Index N>
1815Tensor<typename Promote<S, T>::type, N>
1816dot(Tensor<S, N> const & A, Tensor<T, N> const & B)
1817{
1818 Index const
1819 dimension = A.get_dimension();
1820
1821 assert(B.get_dimension() == dimension);
1822
1824 C(dimension);
1825
1826 switch (dimension) {
1827
1828 default:
1829 for (Index i = 0; i < dimension; ++i) {
1830 for (Index j = 0; j < dimension; ++j) {
1831
1832 typename Promote<S, T>::type
1833 s = 0.0;
1834
1835 for (Index p = 0; p < dimension; ++p) {
1836 s += A(i, p) * B(p, j);
1837 }
1838 C(i, j) = s;
1839 }
1840 }
1841 break;
1842
1843 case 3:
1844 C(0, 0) = A(0, 0) * B(0, 0) + A(0, 1) * B(1, 0) + A(0, 2) * B(2, 0);
1845 C(0, 1) = A(0, 0) * B(0, 1) + A(0, 1) * B(1, 1) + A(0, 2) * B(2, 1);
1846 C(0, 2) = A(0, 0) * B(0, 2) + A(0, 1) * B(1, 2) + A(0, 2) * B(2, 2);
1847
1848 C(1, 0) = A(1, 0) * B(0, 0) + A(1, 1) * B(1, 0) + A(1, 2) * B(2, 0);
1849 C(1, 1) = A(1, 0) * B(0, 1) + A(1, 1) * B(1, 1) + A(1, 2) * B(2, 1);
1850 C(1, 2) = A(1, 0) * B(0, 2) + A(1, 1) * B(1, 2) + A(1, 2) * B(2, 2);
1851
1852 C(2, 0) = A(2, 0) * B(0, 0) + A(2, 1) * B(1, 0) + A(2, 2) * B(2, 0);
1853 C(2, 1) = A(2, 0) * B(0, 1) + A(2, 1) * B(1, 1) + A(2, 2) * B(2, 1);
1854 C(2, 2) = A(2, 0) * B(0, 2) + A(2, 1) * B(1, 2) + A(2, 2) * B(2, 2);
1855 break;
1856
1857 case 2:
1858 C(0, 0) = A(0, 0) * B(0, 0) + A(0, 1) * B(1, 0);
1859 C(0, 1) = A(0, 0) * B(0, 1) + A(0, 1) * B(1, 1);
1860
1861 C(1, 0) = A(1, 0) * B(0, 0) + A(1, 1) * B(1, 0);
1862 C(1, 1) = A(1, 0) * B(0, 1) + A(1, 1) * B(1, 1);
1863 break;
1864
1865 }
1866
1867 return C;
1868}
1869
1870//
1871// Tensor tensor product C = A^T B
1872//
1873template<typename S, typename T, Index N>
1875Tensor<typename Promote<S, T>::type, N>
1876t_dot(Tensor<S, N> const & A, Tensor<T, N> const & B)
1877{
1878 Index const
1879 dimension = A.get_dimension();
1880
1881 assert(B.get_dimension() == dimension);
1882
1884 C(dimension);
1885
1886 switch (dimension) {
1887
1888 default:
1889 for (Index i = 0; i < dimension; ++i) {
1890 for (Index j = 0; j < dimension; ++j) {
1891
1892 typename Promote<S, T>::type
1893 s = 0.0;
1894
1895 for (Index p = 0; p < dimension; ++p) {
1896 s += A(p, i) * B(p, j);
1897 }
1898 C(i, j) = s;
1899 }
1900 }
1901 break;
1902
1903 case 3:
1904 C(0, 0) = A(0, 0) * B(0, 0) + A(1, 0) * B(1, 0) + A(2, 0) * B(2, 0);
1905 C(0, 1) = A(0, 0) * B(0, 1) + A(1, 0) * B(1, 1) + A(2, 0) * B(2, 1);
1906 C(0, 2) = A(0, 0) * B(0, 2) + A(1, 0) * B(1, 2) + A(2, 0) * B(2, 2);
1907
1908 C(1, 0) = A(0, 1) * B(0, 0) + A(1, 1) * B(1, 0) + A(2, 1) * B(2, 0);
1909 C(1, 1) = A(0, 1) * B(0, 1) + A(1, 1) * B(1, 1) + A(2, 1) * B(2, 1);
1910 C(1, 2) = A(0, 1) * B(0, 2) + A(1, 1) * B(1, 2) + A(2, 1) * B(2, 2);
1911
1912 C(2, 0) = A(0, 2) * B(0, 0) + A(1, 2) * B(1, 0) + A(2, 2) * B(2, 0);
1913 C(2, 1) = A(0, 2) * B(0, 1) + A(1, 2) * B(1, 1) + A(2, 2) * B(2, 1);
1914 C(2, 2) = A(0, 2) * B(0, 2) + A(1, 2) * B(1, 2) + A(2, 2) * B(2, 2);
1915 break;
1916
1917 case 2:
1918 C(0, 0) = A(0, 0) * B(0, 0) + A(1, 0) * B(1, 0);
1919 C(0, 1) = A(0, 0) * B(0, 1) + A(1, 0) * B(1, 1);
1920
1921 C(1, 0) = A(0, 1) * B(0, 0) + A(1, 1) * B(1, 0);
1922 C(1, 1) = A(0, 1) * B(0, 1) + A(1, 1) * B(1, 1);
1923 break;
1924
1925 }
1926
1927 return C;
1928}
1929
1930//
1931// Tensor matrix product C = A^T B
1932//
1933template<typename S, typename T, Index N, Index P>
1935Matrix<typename Promote<S, T>::type, N, P>
1937{
1938 Index const
1939 dimension{A.get_dimension()};
1940
1941 Index const
1942 num_cols{B.get_num_cols()};
1943
1944 assert(B.get_num_rows() == dimension);
1945
1947 C(dimension, num_cols);
1948
1949 for (Index i{0}; i < dimension; ++i) {
1950 for (Index j{0}; j < num_cols; ++j) {
1951
1952 typename Promote<S, T>::type
1953 s = 0.0;
1954
1955 for (Index p{0}; p < dimension; ++p) {
1956 s += A(p, i) * B(p, j);
1957 }
1958 C(i, j) = s;
1959 }
1960 }
1961
1962 return C;
1963}
1964
1965//
1966// Tensor vector product C = A^T B
1967//
1968template<typename S, typename T, Index N>
1970Vector<typename Promote<S, T>::type, N>
1971t_dot(Tensor<S, N> const & A, Vector<T, N> const & B)
1972{
1973 return dot(B, A);
1974}
1975
1976//
1977// Tensor tensor product C = A B^T
1978//
1979template<typename S, typename T, Index N>
1981Tensor<typename Promote<S, T>::type, N>
1982dot_t(Tensor<S, N> const & A, Tensor<T, N> const & B)
1983{
1984 Index const
1985 dimension = A.get_dimension();
1986
1987 assert(B.get_dimension() == dimension);
1988
1990 C(dimension);
1991
1992 switch (dimension) {
1993
1994 default:
1995 for (Index i = 0; i < dimension; ++i) {
1996 for (Index j = 0; j < dimension; ++j) {
1997
1998 typename Promote<S, T>::type
1999 s = 0.0;
2000
2001 for (Index p = 0; p < dimension; ++p) {
2002 s += A(i, p) * B(j, p);
2003 }
2004 C(i, j) = s;
2005 }
2006 }
2007 break;
2008
2009 case 3:
2010 C(0, 0) = A(0, 0) * B(0, 0) + A(0, 1) * B(0, 1) + A(0, 2) * B(0, 2);
2011 C(0, 1) = A(0, 0) * B(1, 0) + A(0, 1) * B(1, 1) + A(0, 2) * B(1, 2);
2012 C(0, 2) = A(0, 0) * B(2, 0) + A(0, 1) * B(2, 1) + A(0, 2) * B(2, 2);
2013
2014 C(1, 0) = A(1, 0) * B(0, 0) + A(1, 1) * B(0, 1) + A(1, 2) * B(0, 2);
2015 C(1, 1) = A(1, 0) * B(1, 0) + A(1, 1) * B(1, 1) + A(1, 2) * B(1, 2);
2016 C(1, 2) = A(1, 0) * B(2, 0) + A(1, 1) * B(2, 1) + A(1, 2) * B(2, 2);
2017
2018 C(2, 0) = A(2, 0) * B(0, 0) + A(2, 1) * B(0, 1) + A(2, 2) * B(0, 2);
2019 C(2, 1) = A(2, 0) * B(1, 0) + A(2, 1) * B(1, 1) + A(2, 2) * B(1, 2);
2020 C(2, 2) = A(2, 0) * B(2, 0) + A(2, 1) * B(2, 1) + A(2, 2) * B(2, 2);
2021 break;
2022
2023 case 2:
2024 C(0, 0) = A(0, 0) * B(0, 0) + A(0, 1) * B(0, 1);
2025 C(0, 1) = A(0, 0) * B(1, 0) + A(0, 1) * B(1, 1);
2026
2027 C(1, 0) = A(1, 0) * B(0, 0) + A(1, 1) * B(0, 1);
2028 C(1, 1) = A(1, 0) * B(1, 0) + A(1, 1) * B(1, 1);
2029 break;
2030
2031 }
2032
2033 return C;
2034}
2035
2036//
2037// Tensor tensor product C = A^T B^T
2038//
2039template<typename S, typename T, Index N>
2041Tensor<typename Promote<S, T>::type, N>
2042t_dot_t(Tensor<S, N> const & A, Tensor<T, N> const & B)
2043{
2044 Index const
2045 dimension = A.get_dimension();
2046
2047 assert(B.get_dimension() == dimension);
2048
2050 C(dimension);
2051
2052 switch (dimension) {
2053
2054 default:
2055 for (Index i = 0; i < dimension; ++i) {
2056 for (Index j = 0; j < dimension; ++j) {
2057
2058 typename Promote<S, T>::type
2059 s = 0.0;
2060
2061 for (Index p = 0; p < dimension; ++p) {
2062 s += A(p, i) * B(j, p);
2063 }
2064 C(i, j) = s;
2065 }
2066 }
2067 break;
2068
2069 case 3:
2070 C(0, 0) = A(0, 0) * B(0, 0) + A(1, 0) * B(0, 1) + A(2, 0) * B(0, 2);
2071 C(0, 1) = A(0, 0) * B(1, 0) + A(1, 0) * B(1, 1) + A(2, 0) * B(1, 2);
2072 C(0, 2) = A(0, 0) * B(2, 0) + A(1, 0) * B(2, 1) + A(2, 0) * B(2, 2);
2073
2074 C(1, 0) = A(0, 1) * B(0, 0) + A(1, 1) * B(0, 1) + A(2, 1) * B(0, 2);
2075 C(1, 1) = A(0, 1) * B(1, 0) + A(1, 1) * B(1, 1) + A(2, 1) * B(1, 2);
2076 C(1, 2) = A(0, 1) * B(2, 0) + A(1, 1) * B(2, 1) + A(2, 1) * B(2, 2);
2077
2078 C(2, 0) = A(0, 2) * B(0, 0) + A(1, 2) * B(0, 1) + A(2, 2) * B(0, 2);
2079 C(2, 1) = A(0, 2) * B(1, 0) + A(1, 2) * B(1, 1) + A(2, 2) * B(1, 2);
2080 C(2, 2) = A(0, 2) * B(2, 0) + A(1, 2) * B(2, 1) + A(2, 2) * B(2, 2);
2081 break;
2082
2083 case 2:
2084 C(0, 0) = A(0, 0) * B(0, 0) + A(1, 0) * B(0, 1);
2085 C(0, 1) = A(0, 0) * B(1, 0) + A(1, 0) * B(1, 1);
2086
2087 C(1, 0) = A(0, 1) * B(0, 0) + A(1, 1) * B(0, 1);
2088 C(1, 1) = A(0, 1) * B(1, 0) + A(1, 1) * B(1, 1);
2089 break;
2090
2091 }
2092
2093 return C;
2094}
2095
2096//
2097// Tensor tensor double dot product (contraction)
2098//
2099template<typename S, typename T, Index N>
2101typename Promote<S, T>::type
2102dotdot(Tensor<S, N> const & A, Tensor<T, N> const & B)
2103{
2104 Index const
2105 dimension = A.get_dimension();
2106
2107 assert(B.get_dimension() == dimension);
2108
2109 typename Promote<S, T>::type
2110 s = 0.0;
2111
2112 switch (dimension) {
2113
2114 default:
2115 for (Index p = 0; p < dimension; ++p) {
2116 for (Index q = 0; q < dimension; ++q) {
2117 s += A(p, q) * B(p, q);
2118 }
2119 }
2120 break;
2121
2122 case 3:
2123 s += A(0, 0) * B(0, 0) + A(0, 1) * B(0, 1) + A(0, 2) * B(0, 2);
2124 s += A(1, 0) * B(1, 0) + A(1, 1) * B(1, 1) + A(1, 2) * B(1, 2);
2125 s += A(2, 0) * B(2, 0) + A(2, 1) * B(2, 1) + A(2, 2) * B(2, 2);
2126 break;
2127
2128 case 2:
2129 s += A(0, 0) * B(0, 0) + A(0, 1) * B(0, 1);
2130 s += A(1, 0) * B(1, 0) + A(1, 1) * B(1, 1);
2131 break;
2132
2133 }
2134
2135 return s;
2136}
2137
2138//
2139// dyad
2140//
2141template<typename S, typename T, Index N>
2143Tensor<typename Promote<S, T>::type, N>
2144dyad(Vector<S, N> const & u, Vector<T, N> const & v)
2145{
2146 Index const
2147 dimension = u.get_dimension();
2148
2149 assert(v.get_dimension() == dimension);
2150
2152 A(dimension);
2153
2154 switch (dimension) {
2155
2156 default:
2157 for (Index i = 0; i < dimension; ++i) {
2158
2159 typename Promote<S, T>::type const
2160 s = u(i);
2161
2162 for (Index j = 0; j < dimension; ++j) {
2163 A(i, j) = s * v(j);
2164 }
2165 }
2166 break;
2167
2168 case 3:
2169 A(0, 0) = u(0) * v(0);
2170 A(0, 1) = u(0) * v(1);
2171 A(0, 2) = u(0) * v(2);
2172
2173 A(1, 0) = u(1) * v(0);
2174 A(1, 1) = u(1) * v(1);
2175 A(1, 2) = u(1) * v(2);
2176
2177 A(2, 0) = u(2) * v(0);
2178 A(2, 1) = u(2) * v(1);
2179 A(2, 2) = u(2) * v(2);
2180 break;
2181
2182 case 2:
2183 A(0, 0) = u(0) * v(0);
2184 A(0, 1) = u(0) * v(1);
2185
2186 A(1, 0) = u(1) * v(0);
2187 A(1, 1) = u(1) * v(1);
2188 break;
2189
2190 }
2191
2192 return A;
2193}
2194
2195//
2196// bun operator, just for Jay, and now Reese too.
2197//
2198template<typename S, typename T, Index N>
2200Tensor<typename Promote<S, T>::type, N>
2201bun(Vector<S, N> const & u, Vector<T, N> const & v)
2202{
2203 return dyad(u, v);
2204}
2205
2206//
2207// tensor product
2208//
2209template<typename S, typename T, Index N>
2211Tensor<typename Promote<S, T>::type, N>
2212tensor(Vector<S, N> const & u, Vector<T, N> const & v)
2213{
2214 return dyad(u, v);
2215}
2216
2217//
2218// diagonal tensor from vector
2219//
2220template<typename T, Index N>
2222Tensor<T, N>
2224{
2225 Index const
2226 dimension = v.get_dimension();
2227
2229 A = zero<T, N>(dimension);
2230
2231 switch (dimension) {
2232
2233 default:
2234 for (Index i = 0; i < dimension; ++i) {
2235 A(i, i) = v(i);
2236 }
2237 break;
2238
2239 case 3:
2240 A(0, 0) = v(0);
2241 A(1, 1) = v(1);
2242 A(2, 2) = v(2);
2243 break;
2244
2245 case 2:
2246 A(0, 0) = v(0);
2247 A(1, 1) = v(1);
2248 break;
2249
2250 }
2251
2252 return A;
2253}
2254
2255//
2256// diagonal of tensor in a vector
2257//
2258template<typename T, Index N>
2260Vector<T, N>
2262{
2263 Index const
2264 dimension = A.get_dimension();
2265
2267 v(dimension);
2268
2269 switch (dimension) {
2270
2271 default:
2272 for (Index i = 0; i < dimension; ++i) {
2273 v(i) = A(i, i);
2274 }
2275 break;
2276
2277 case 3:
2278 v(0) = A(0, 0);
2279 v(1) = A(1, 1);
2280 v(2) = A(2, 2);
2281 break;
2282
2283 case 2:
2284 v(0) = A(0, 0);
2285 v(1) = A(1, 1);
2286 break;
2287
2288 }
2289
2290 return v;
2291}
2292
2293//
2294// Zero 2nd-order tensor
2295// All components are zero
2296//
2297template<typename T, Index N>
2299Tensor<T, N> const
2301{
2302 return Tensor<T, N>(N, Filler::ZEROS);
2303}
2304
2308template<typename T>
2310Tensor<T, DYNAMIC> const
2311zero(Index const dimension)
2312{
2313 return Tensor<T, DYNAMIC>(dimension, Filler::ZEROS);
2314}
2315
2319template<typename T, Index N>
2321Tensor<T, N> const
2322zero(Index const dimension)
2323{
2324 return Tensor<T, N>(dimension, Filler::ZEROS);
2325}
2326
2327// Local utility functions
2328namespace {
2329
2330template<typename T, Index N>
2332void ones_in_diagonal(Tensor<T, N> & A)
2333{
2334 Index const
2335 dimension = A.get_dimension();
2336
2337 switch (dimension) {
2338
2339 default:
2340 for (Index i = 0; i < dimension; ++i) {
2341 A(i, i) = 1.0;
2342 }
2343 break;
2344
2345 case 3:
2346 A(0, 0) = 1.0;
2347 A(1, 1) = 1.0;
2348 A(2, 2) = 1.0;
2349 break;
2350
2351 case 2:
2352 A(0, 0) = 1.0;
2353 A(1, 1) = 1.0;
2354 break;
2355
2356 }
2357
2358 return;
2359}
2360
2361template<typename T, Index N>
2363void fill_levi_civita(Tensor<T, N> & A)
2364{
2365 Index const
2366 dimension = A.get_dimension();
2367
2368 for (Index i = 0; i < dimension; ++i) {
2369 for (Index j = 0; j < dimension; ++j) {
2370 A(i, j) = levi_civita<T>(i, j);
2371 }
2372 }
2373
2374 return;
2375}
2376
2377} // anonymous namespace
2378
2379//
2380// 2nd-order identity tensor
2381//
2382template<typename T, Index N>
2384Tensor<T, N> const
2386{
2388 A(N, Filler::ZEROS);
2389
2390 ones_in_diagonal(A);
2391
2392 return A;
2393}
2394
2398template<typename T>
2400Tensor<T, DYNAMIC> const
2401identity(Index const dimension)
2402{
2404 A(dimension, Filler::ZEROS);
2405
2406 ones_in_diagonal(A);
2407
2408 return A;
2409}
2410
2414template<typename T, Index N>
2416Tensor<T, N> const
2417identity(Index const dimension)
2418{
2419
2421 A(dimension, Filler::ZEROS);
2422
2423 ones_in_diagonal(A);
2424
2425 return A;
2426}
2427
2428//
2429// 2nd-order identity tensor, à la Matlab
2430//
2431template<typename T, Index N>
2433Tensor<T, N> const
2435{
2436 return identity<T, N>();
2437}
2438
2442template<typename T>
2444Tensor<T, DYNAMIC> const
2445eye(Index const dimension)
2446{
2447 return identity<T, DYNAMIC>(dimension);
2448}
2449
2453template<typename T, Index N>
2455Tensor<T, N> const
2456eye(Index const dimension)
2457{
2458 return identity<T, N>(dimension);
2459}
2460
2461//
2462// Levi-Civita symbol
2463//
2464template<typename T, Index N>
2466Tensor<T, N> const
2468{
2470 A(N, Filler::ZEROS);
2471
2472 fill_levi_civita(A);
2473
2474 return A;
2475}
2476
2480template<typename T>
2482Tensor<T, DYNAMIC> const
2483levi_civita_2(Index const dimension)
2484{
2486 A(dimension, Filler::ZEROS);
2487
2488 fill_levi_civita(A);
2489
2490 return A;
2491}
2492
2496template<typename T, Index N>
2498Tensor<T, N> const
2499levi_civita_2(Index const dimension)
2500{
2501
2503 A(dimension, Filler::ZEROS);
2504
2505 fill_levi_civita(A);
2506
2507 return A;
2508}
2509
2510//
2511// Permutation symbol
2512//
2513template<typename T, Index N>
2515Tensor<T, N> const
2517{
2518 return levi_civita_2<T, N>();
2519}
2520
2524template<typename T>
2526Tensor<T, DYNAMIC> const
2527permutation_2(Index const dimension)
2528{
2529 return levi_civita_2<T, DYNAMIC>(dimension);
2530}
2531
2532
2536template<typename T, Index N>
2538Tensor<T, N> const
2539permutation_2(Index const dimension)
2540{
2541 return levi_civita_2<T, N>(dimension);
2542}
2543
2544//
2545// Alternating symbol
2546//
2547template<typename T, Index N>
2549Tensor<T, N> const
2551{
2552 return levi_civita_2<T, N>();
2553}
2554
2558template<typename T>
2560Tensor<T, DYNAMIC> const
2561alternator_2(Index const dimension)
2562{
2563 return levi_civita_2<T, DYNAMIC>(dimension);
2564}
2565
2569template<typename T, Index N>
2571Tensor<T, N> const
2572alternator_2(Index const dimension)
2573{
2574 return levi_civita_2<T, N>(dimension);
2575}
2576
2577//
2578// 2nd-order tensor transpose
2579//
2580template<typename T, Index N>
2582Tensor<T, N>
2584{
2585 Index const
2586 dimension = A.get_dimension();
2587
2589 B = A;
2590
2591 switch (dimension) {
2592 default:
2593 for (Index i = 0; i < dimension; ++i) {
2594 for (Index j = i + 1; j < dimension; ++j) {
2595 minitensor::swap(B(i, j), B(j, i));
2596 }
2597 }
2598 break;
2599
2600 case 3:
2601 minitensor::swap(B(0, 1), B(1, 0));
2602 minitensor::swap(B(0, 2), B(2, 0));
2603 minitensor::swap(B(1, 2), B(2, 1));
2604 break;
2605
2606 case 2:
2607 minitensor::swap(B(0, 1), B(1, 0));
2608 break;
2609 }
2610
2611 return B;
2612}
2613
2614//
2615// symmetric part of 2nd-order tensor
2616// \return \f$ \frac{1}{2}(A + A^T) \f$
2617//
2618template<typename T, Index N>
2620Tensor<T, N>
2622{
2623 Index const
2624 dimension = A.get_dimension();
2625
2627 B(dimension);
2628
2629 switch (dimension) {
2630
2631 default:
2632 B = 0.5 * (A + transpose(A));
2633 break;
2634
2635 case 3:
2636 {
2637 T const & s00 = A(0, 0);
2638 T const & s11 = A(1, 1);
2639 T const & s22 = A(2, 2);
2640
2641 T const s01 = 0.5 * (A(0, 1) + A(1, 0));
2642 T const s02 = 0.5 * (A(0, 2) + A(2, 0));
2643 T const s12 = 0.5 * (A(1, 2) + A(2, 1));
2644
2645 B(0, 0) = s00;
2646 B(0, 1) = s01;
2647 B(0, 2) = s02;
2648
2649 B(1, 0) = s01;
2650 B(1, 1) = s11;
2651 B(1, 2) = s12;
2652
2653 B(2, 0) = s02;
2654 B(2, 1) = s12;
2655 B(2, 2) = s22;
2656 }
2657 break;
2658
2659 case 2:
2660 {
2661 T const & s00 = A(0, 0);
2662 T const & s11 = A(1, 1);
2663
2664 T const s01 = 0.5 * (A(0, 1) + A(1, 0));
2665
2666 B(0, 0) = s00;
2667 B(0, 1) = s01;
2668
2669 B(1, 0) = s01;
2670 B(1, 1) = s11;
2671 }
2672 break;
2673
2674 }
2675
2676 return B;
2677}
2678
2679//
2680// skew symmetric part of 2nd-order tensor
2681// \return \f$ \frac{1}{2}(A - A^T) \f$
2682//
2683template<typename T, Index N>
2685Tensor<T, N>
2687{
2688 Index const
2689 dimension = A.get_dimension();
2690
2692 B(dimension);
2693
2694 switch (dimension) {
2695
2696 default:
2697 B = 0.5 * (A - transpose(A));
2698 break;
2699
2700 case 3:
2701 {
2702 T const s01 = 0.5 * (A(0, 1) - A(1, 0));
2703 T const s02 = 0.5 * (A(0, 2) - A(2, 0));
2704 T const s12 = 0.5 * (A(1, 2) - A(2, 1));
2705
2706 B(0, 0) = 0.0;
2707 B(0, 1) = s01;
2708 B(0, 2) = s02;
2709
2710 B(1, 0) = -s01;
2711 B(1, 1) = 0.0;
2712 B(1, 2) = s12;
2713
2714 B(2, 0) = -s02;
2715 B(2, 1) = -s12;
2716 B(2, 2) = 0.0;
2717 }
2718 break;
2719
2720 case 2:
2721 {
2722 T const s01 = 0.5 * (A(0, 1) - A(1, 0));
2723
2724 B(0, 0) = 0.0;
2725 B(0, 1) = s01;
2726
2727 B(1, 0) = -s01;
2728 B(1, 1) = 0.0;
2729 }
2730 break;
2731
2732 }
2733
2734 return B;
2735}
2736
2737//
2738// skew symmetric 2nd-order tensor from vector, undefined
2739// for N!=3.
2740// \param u vector
2741//
2742template<typename T, Index N>
2744Tensor<T, N>
2746{
2747 Index const
2748 dimension = u.get_dimension();
2749
2751 A(dimension);
2752
2753 switch (dimension) {
2754
2755 case 3:
2756 A(0, 0) = 0.0;
2757 A(0, 1) = -u(2);
2758 A(0, 2) = u(1);
2759
2760 A(1, 0) = u(2);
2761 A(1, 1) = 0.0;
2762 A(1, 2) = -u(0);
2763
2764 A(2, 0) = -u(1);
2765 A(2, 1) = u(0);
2766 A(2, 2) = 0.0;
2767 break;
2768
2769 default:
2770 MT_ERROR_EXIT("Skew from vector defined for 3D only");
2771 break;
2772 }
2773
2774 return A;
2775}
2776
2777} // namespace minitensor
2778namespace minitensor {
2779
2780//
2781// tensor input
2782//
2783template<typename T, Index N>
2784std::istream &
2785operator>>(std::istream & is, Tensor<T, N> & A)
2786{
2787
2788 Index const
2789 dimension = A.get_dimension();
2790
2791 for (Index i = 0; i < dimension; ++i) {
2792 for (Index j = 0; j < dimension; ++j) {
2793 is >> A(i,j);
2794 }
2795 }
2796
2797 return is;
2798}
2799
2800//
2801// tensor output
2802//
2803template<typename T, Index N>
2804std::ostream &
2805operator<<(std::ostream & os, Tensor<T, N> const & A)
2806{
2807 Index const
2808 dimension = A.get_dimension();
2809
2810 if (dimension == 0) {
2811 return os;
2812 }
2813
2814 os << std::scientific << std::setprecision(17);
2815
2816 for (Index i = 0; i < dimension; ++i) {
2817
2818 os << std::setw(24) << A(i,0);
2819
2820 for (Index j = 1; j < dimension; ++j) {
2821 os << "," << std::setw(24) << A(i,j);
2822 }
2823
2824 os << std::endl;
2825 }
2826
2827 return os;
2828}
2829
2831} // namespace minitensor
2832
2833#endif //MiniTensor_Tensor_h
#define KOKKOS_INLINE_FUNCTION
#define MT_ERROR_EXIT(...)
KOKKOS_INLINE_FUNCTION Tensor(Source const source, Index const dimension, ArrayT &data, Index index1)
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > operator+(Matrix< S, M, N > const &A, Matrix< T, M, N > const &B)
KOKKOS_INLINE_FUNCTION Tensor(Source const source, Index const dimension, ArrayT &data, Index index1, Index index2)
KOKKOS_INLINE_FUNCTION Tensor(Source const source, ArrayT &data, Index index1, Index index2)
KOKKOS_INLINE_FUNCTION void minus(TensorBase< T, ST > const &A, TensorBase< T, ST > &B)
KOKKOS_INLINE_FUNCTION Tensor< T, N > const eye()
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > operator-(Matrix< S, M, N > const &A, Matrix< T, M, N > const &B)
std::ostream & operator<<(std::ostream &os, Matrix< T, M, N > const &A)
KOKKOS_INLINE_FUNCTION void add(TensorBase< R, SR > const &A, TensorBase< S, SS > const &B, TensorBase< T, ST > &C)
KOKKOS_INLINE_FUNCTION Tensor(T const *data_ptr, ComponentOrder const component_order)
KOKKOS_INLINE_FUNCTION Index get_num_cols() const
std::istream & operator>>(std::istream &is, Matrix< T, M, N > &A)
KOKKOS_INLINE_FUNCTION Tensor< typename Promote< S, T >::type, N > dyad(Vector< S, N > const &u, Vector< T, N > const &v)
KOKKOS_INLINE_FUNCTION Tensor< T, N > const alternator_2()
KOKKOS_INLINE_FUNCTION Tensor< T, N > const identity()
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > operator/(Matrix< T, M, N > const &A, S const &s)
KOKKOS_INLINE_FUNCTION Vector< typename Promote< S, T >::type, M > operator*(Matrix< T, M, N > const &A, Vector< S, N > const &u)
KOKKOS_INLINE_FUNCTION Tensor< T, N > const permutation_2()
KOKKOS_INLINE_FUNCTION Promote< S, T >::type dotdot(Matrix< S, M, N > const &A, Matrix< T, M, N > const &B)
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > t_dot_t(Matrix< S, P, M > const &A, Matrix< T, N, P > const &B)
KOKKOS_INLINE_FUNCTION Tensor()
KOKKOS_INLINE_FUNCTION Tensor(T const &s00, T const &s01, T const &s10, T const &s11)
KOKKOS_INLINE_FUNCTION Tensor< typename Promote< S, T >::type, N > bun(Vector< S, N > const &u, Vector< T, N > const &v)
KOKKOS_INLINE_FUNCTION Tensor(Source const source, Index const dimension, ArrayT &data, Index index1, Index index2, Index index3, Index index4)
KOKKOS_INLINE_FUNCTION Tensor< T, N > const zero()
KOKKOS_INLINE_FUNCTION Index get_dimension() const
KOKKOS_INLINE_FUNCTION Tensor(Index const dimension, T const *data_ptr, ComponentOrder const component_order)
KOKKOS_INLINE_FUNCTION Tensor(Source const source, ArrayT &data, Index index1, Index index2, Index index3, Index index4, Index index5, Index index6)
KOKKOS_INLINE_FUNCTION bool equal(TensorBase< T, ST > const &A, TensorBase< T, ST > const &B)
KOKKOS_INLINE_FUNCTION Vector< T, M > col(Matrix< T, M, N > const &A, Index const j)
KOKKOS_INLINE_FUNCTION Index get_dimension() const
KOKKOS_INLINE_FUNCTION Index get_dimension(Index const order) const
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > dot_t(Matrix< S, M, P > const &A, Matrix< T, N, P > const &B)
KOKKOS_INLINE_FUNCTION T & operator()(Index const i, Index const j)
KOKKOS_INLINE_FUNCTION Tensor(Source const source, ArrayT &data, Index index1)
KOKKOS_INLINE_FUNCTION Tensor< T, N > diag(Vector< T, N > const &v)
KOKKOS_INLINE_FUNCTION Tensor< T, N > const levi_civita_2()
KOKKOS_INLINE_FUNCTION Tensor(Source const source, Index const dimension, ArrayT &data, Index index1, Index index2, Index index3)
KOKKOS_INLINE_FUNCTION void scale(TensorBase< R, SR > const &A, S const &s, TensorBase< T, ST > &B)
KOKKOS_INLINE_FUNCTION bool operator!=(Matrix< T, M, N > const &A, Matrix< T, M, N > const &B)
KOKKOS_INLINE_FUNCTION void set_dimension(Index const dimension)
KOKKOS_INLINE_FUNCTION Matrix< T, M, N > adjoint(Matrix< T, N, M > const &A)
KOKKOS_INLINE_FUNCTION Tensor(Index const dimension)
KOKKOS_INLINE_FUNCTION Tensor(Source const source, ArrayT &data, Index index1, Index index2, Index index3, Index index4, Index index5)
KOKKOS_INLINE_FUNCTION bool operator==(Matrix< T, M, N > const &A, Matrix< T, M, N > const &B)
KOKKOS_INLINE_FUNCTION Vector< typename Promote< S, T >::type, M > dot(Matrix< T, M, N > const &A, Vector< S, N > const &u)
KOKKOS_INLINE_FUNCTION Tensor(Tensor4< T, dimension_sqrt< N >::value > const &A)
KOKKOS_INLINE_FUNCTION Index get_num_cols() const
KOKKOS_INLINE_FUNCTION Tensor(Filler const value)
KOKKOS_INLINE_FUNCTION Tensor(T const *data_ptr)
KOKKOS_INLINE_FUNCTION void split(TensorBase< R, SR > const &A, S const &s, TensorBase< T, ST > &B)
KOKKOS_INLINE_FUNCTION Index get_num_rows() const
KOKKOS_INLINE_FUNCTION Tensor(T const &s00, T const &s01, T const &s02, T const &s10, T const &s11, T const &s12, T const &s20, T const &s21, T const &s22)
static constexpr bool IS_DYNAMIC
KOKKOS_INLINE_FUNCTION Tensor(Source const source, Index const dimension, ArrayT &data, Index index1, Index index2, Index index3, Index index4, Index index5, Index index6)
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > t_dot(Matrix< S, P, M > const &A, Matrix< T, P, N > const &B)
KOKKOS_INLINE_FUNCTION Tensor< typename Promote< S, T >::type, N > tensor(Vector< S, N > const &u, Vector< T, N > const &v)
KOKKOS_INLINE_FUNCTION Tensor(Source const source, ArrayT &data, Index index1, Index index2, Index index3, Index index4)
KOKKOS_INLINE_FUNCTION Tensor(Index const dimension, Filler const value)
KOKKOS_INLINE_FUNCTION Tensor(Source const source, Index const dimension, ArrayT &data, Index index1, Index index2, Index index3, Index index4, Index index5)
virtual KOKKOS_INLINE_FUNCTION ~Tensor()
KOKKOS_INLINE_FUNCTION bool not_equal(TensorBase< T, ST > const &A, TensorBase< T, ST > const &B)
static constexpr Index ORDER
KOKKOS_INLINE_FUNCTION Tensor(Index const dimension, T const *data_ptr)
KOKKOS_INLINE_FUNCTION void set_dimension(Index const dimension, Index const order)
KOKKOS_INLINE_FUNCTION Index get_num_rows() const
KOKKOS_INLINE_FUNCTION Tensor(Source const source, ArrayT &data, Index index1, Index index2, Index index3)
KOKKOS_INLINE_FUNCTION Matrix< T, M, N > transpose(Matrix< T, N, M > const &A)
KOKKOS_INLINE_FUNCTION Tensor(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION Tensor< T, N > sym(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION void divide(TensorBase< R, SR > const &A, S const &s, TensorBase< T, ST > &B)
static KOKKOS_INLINE_FUNCTION constexpr Index get_order()
KOKKOS_INLINE_FUNCTION Tensor< T, N > skew(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION Vector< T, N > row(Matrix< T, M, N > const &A, Index const i)
KOKKOS_INLINE_FUNCTION T const & operator()(Index const i, Index const j) const
KOKKOS_INLINE_FUNCTION void subtract(TensorBase< R, SR > const &A, TensorBase< S, SS > const &B, TensorBase< T, ST > &C)
uint32_t Index
Indexing type.
constexpr Index DYNAMIC
Indicator for dynamic storage.
KOKKOS_INLINE_FUNCTION void swap(T &a, T &b)