1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributor(s):
7 // Eric Bechet
8 //
9
10 #ifndef STENSOR43_H
11 #define STENSOR43_H
12
13 #include "STensor33.h"
14 #include "fullMatrix.h"
15 #include "Numeric.h"
16
17 // concrete class for general 3rd-order tensor in three-dimensional space
18
19 class STensor43 {
20 protected:
21 // 0000 0001 0002 0010 ... 2211 2212 2220 2221 2222
22 double _val[81];
23
24 public:
getIndex(int i,int j,int k,int l)25 inline int getIndex(int i, int j, int k, int l) const
26 {
27 static int _index[3][3][3][3] = {
28 {{{0, 1, 2}, {3, 4, 5}, {6, 7, 8}},
29 {{9, 10, 11}, {12, 13, 14}, {15, 16, 17}},
30 {{18, 19, 20}, {21, 22, 23}, {24, 25, 26}}},
31 {{{27, 28, 29}, {30, 31, 32}, {33, 34, 35}},
32 {{36, 37, 38}, {39, 40, 41}, {42, 43, 44}},
33 {{45, 46, 47}, {48, 49, 50}, {51, 52, 53}}},
34 {{{54, 55, 56}, {57, 58, 59}, {60, 61, 62}},
35 {{63, 64, 65}, {66, 67, 68}, {69, 70, 71}},
36 {{72, 73, 74}, {75, 76, 77}, {78, 79, 80}}}};
37 return _index[i][j][k][l];
38 }
STensor43(const STensor43 & other)39 STensor43(const STensor43 &other)
40 {
41 for(int i = 0; i < 81; i++) _val[i] = other._val[i];
42 }
43 // default constructor, null tensor
44 STensor43(const double v = 0.0)
45 {
46 for(int i = 0; i < 3; i++)
47 for(int j = 0; j < 3; j++)
48 for(int k = 0; k < 3; k++)
49 for(int l = 0; l < 3; l++)
50 if((i == k) && (j == l))
51 _val[getIndex(i, j, k, l)] = v;
52 else
53 _val[getIndex(i, j, k, l)] = 0.0;
54 }
55 // Symmetric identity tensor
STensor43(const double vik,const double vil)56 STensor43(const double vik, const double vil)
57 {
58 for(int i = 0; i < 3; i++)
59 for(int j = 0; j < 3; j++)
60 for(int k = 0; k < 3; k++)
61 for(int l = 0; l < 3; l++) {
62 _val[getIndex(i, j, k, l)] = 0.;
63 if((i == k) && (j == l)) _val[getIndex(i, j, k, l)] += 0.5 * vik;
64 if((i == l) && (j == k)) _val[getIndex(i, j, k, l)] += 0.5 * vil;
65 }
66 }
operator()67 inline double &operator()(int i, int j, int k, int l)
68 {
69 return _val[getIndex(i, j, k, l)];
70 }
operator()71 inline double operator()(int i, int j, int k, int l) const
72 {
73 return _val[getIndex(i, j, k, l)];
74 }
75 STensor43 operator+(const STensor43 &other) const
76 {
77 STensor43 res(*this);
78 for(int i = 0; i < 81; i++) res._val[i] += other._val[i];
79 return res;
80 }
81 STensor43 &operator+=(const STensor43 &other)
82 {
83 for(int i = 0; i < 81; i++) _val[i] += other._val[i];
84 return *this;
85 }
86 STensor43 &operator-=(const STensor43 &other)
87 {
88 for(int i = 0; i < 81; i++) _val[i] -= other._val[i];
89 return *this;
90 }
91 STensor43 &operator=(const STensor43 &other)
92 {
93 for(int i = 0; i < 81; i++) _val[i] = other._val[i];
94 return *this;
95 }
96 STensor43 &operator*=(const double &other)
97 {
98 for(int i = 0; i < 81; i++) _val[i] *= other;
99 return *this;
100 }
101
transpose(int n,int m)102 STensor43 transpose(int n, int m) const
103 {
104 STensor43 ithis;
105 if((n == 0 && m == 1) || (n == 1 && m == 0)) {
106 for(int i = 0; i < 3; i++)
107 for(int j = 0; j < 3; j++)
108 for(int k = 0; k < 3; k++)
109 for(int l = 0; l < 3; l++) ithis(i, j, k, l) = (*this)(j, i, k, l);
110 return ithis;
111 }
112 if((n == 0 && m == 2) || (n == 2 && m == 0)) {
113 for(int i = 0; i < 3; i++)
114 for(int j = 0; j < 3; j++)
115 for(int k = 0; k < 3; k++)
116 for(int l = 0; l < 3; l++) ithis(i, j, k, l) = (*this)(k, j, i, l);
117 return ithis;
118 }
119 if((n == 0 && m == 3) || (n == 3 && m == 0)) {
120 for(int i = 0; i < 3; i++)
121 for(int j = 0; j < 3; j++)
122 for(int k = 0; k < 3; k++)
123 for(int l = 0; l < 3; l++) ithis(i, j, k, l) = (*this)(l, j, k, i);
124 return ithis;
125 }
126 if((n == 1 && m == 2) || (n == 2 && m == 1)) {
127 for(int i = 0; i < 3; i++)
128 for(int j = 0; j < 3; j++)
129 for(int k = 0; k < 3; k++)
130 for(int l = 0; l < 3; l++) ithis(i, j, k, l) = (*this)(i, k, j, l);
131 return ithis;
132 }
133 if((n == 1 && m == 3) || (n == 3 && m == 1)) {
134 for(int i = 0; i < 3; i++)
135 for(int j = 0; j < 3; j++)
136 for(int k = 0; k < 3; k++)
137 for(int l = 0; l < 3; l++) ithis(i, j, k, l) = (*this)(i, l, k, j);
138 return ithis;
139 }
140 if((n == 2 && m == 3) || (n == 3 && m == 2)) {
141 for(int i = 0; i < 3; i++)
142 for(int j = 0; j < 3; j++)
143 for(int k = 0; k < 3; k++)
144 for(int l = 0; l < 3; l++) ithis(i, j, k, l) = (*this)(i, j, l, k);
145 return ithis;
146 }
147 return ithis += (*this);
148 }
149 /* STensor43& operator *= (const STensor43 &other)
150 {
151 // to be implemented
152 return *this;
153 }*/
154 void print(const char *) const;
data()155 const double *data() const { return _val; }
data()156 double *data() { return _val; }
157
axpy(const double a,const STensor43 & other)158 void axpy(const double a, const STensor43 &other)
159 {
160 for(int i = 0; i < 81; i++) _val[i] += a * other._val[i];
161 }
162 };
163
164 // tensor product
tensprod(const STensor3 & a,const STensor3 & b,STensor43 & c)165 inline void tensprod(const STensor3 &a, const STensor3 &b, STensor43 &c)
166 {
167 for(int i = 0; i < 3; i++)
168 for(int j = 0; j < 3; j++)
169 for(int k = 0; k < 3; k++)
170 for(int l = 0; l < 3; l++) c(i, j, k, l) = a(i, j) * b(k, l);
171 }
172
tensprod(const SVector3 & a,const STensor33 & b,STensor43 & c)173 inline void tensprod(const SVector3 &a, const STensor33 &b, STensor43 &c)
174 {
175 for(int i = 0; i < 3; i++)
176 for(int j = 0; j < 3; j++)
177 for(int k = 0; k < 3; k++)
178 for(int l = 0; l < 3; l++) c(i, j, k, l) = a(i) * b(j, k, l);
179 }
tensprod(const STensor33 & a,const SVector3 & b,STensor43 & c)180 inline void tensprod(const STensor33 &a, const SVector3 &b, STensor43 &c)
181 {
182 for(int i = 0; i < 3; i++)
183 for(int j = 0; j < 3; j++)
184 for(int k = 0; k < 3; k++)
185 for(int l = 0; l < 3; l++) c(i, j, k, l) = a(i, j, k) * b(l);
186 }
187
dot(const STensor43 & a,const STensor43 & b)188 inline double dot(const STensor43 &a, const STensor43 &b)
189 {
190 double prod = 0;
191 for(int i = 0; i < 3; i++)
192 for(int j = 0; j < 3; j++)
193 for(int k = 0; k < 3; k++)
194 for(int l = 0; l < 3; l++) prod += a(i, j, k, l) * b(i, j, k, l);
195 return prod;
196 }
197
198 // full contracted product
199 inline STensor43 operator*(const STensor43 &t, double m)
200 {
201 STensor43 val(t);
202 val *= m;
203 return val;
204 }
205 inline STensor43 operator*(double m, const STensor43 &t)
206 {
207 STensor43 val(t);
208 val *= m;
209 return val;
210 }
211
212 inline STensor33 operator*(const STensor43 &t, const SVector3 &m)
213 {
214 STensor33 val(0.);
215 for(int i = 0; i < 3; i++)
216 for(int j = 0; j < 3; j++)
217 for(int k = 0; k < 3; k++)
218 for(int l = 0; l < 3; l++) val(i, j, k) += t(i, j, k, l) * m(l);
219 return val;
220 }
221 inline STensor33 operator*(const SVector3 &m, const STensor43 &t)
222 {
223 STensor33 val(0.);
224 for(int i = 0; i < 3; i++)
225 for(int j = 0; j < 3; j++)
226 for(int k = 0; k < 3; k++)
227 for(int l = 0; l < 3; l++) val(j, k, l) += m(i) * t(i, j, k, l);
228 return val;
229 }
230
231 inline STensor3 operator*(const STensor43 &t, const STensor3 &m)
232 {
233 STensor3 val(0.);
234 for(int i = 0; i < 3; i++)
235 for(int j = 0; j < 3; j++)
236 for(int k = 0; k < 3; k++)
237 for(int l = 0; l < 3; l++) val(i, j) += t(i, j, k, l) * m(l, k);
238 return val;
239 }
240 inline STensor3 operator*(const STensor3 &m, const STensor43 &t)
241 {
242 STensor3 val(0.);
243 for(int i = 0; i < 3; i++)
244 for(int j = 0; j < 3; j++)
245 for(int k = 0; k < 3; k++)
246 for(int l = 0; l < 3; l++) val(k, l) += m(j, i) * t(i, j, k, l);
247 return val;
248 }
249
250 inline SVector3 operator*(const STensor43 &t, const STensor33 &m)
251 {
252 SVector3 val(0.);
253 for(int i = 0; i < 3; i++)
254 for(int j = 0; j < 3; j++)
255 for(int k = 0; k < 3; k++)
256 for(int l = 0; l < 3; l++) val(i) += t(i, j, k, l) * m(l, k, j);
257 return val;
258 }
259 inline SVector3 operator*(const STensor33 &m, const STensor43 &t)
260 {
261 SVector3 val(0.);
262 for(int i = 0; i < 3; i++)
263 for(int j = 0; j < 3; j++)
264 for(int k = 0; k < 3; k++)
265 for(int l = 0; l < 3; l++) val(l) += m(k, j, i) * t(i, j, k, l);
266 return val;
267 }
268
269 inline double operator*(const STensor43 &m, const STensor43 &t)
270 {
271 double val(0.);
272 for(int i = 0; i < 3; i++)
273 for(int j = 0; j < 3; j++)
274 for(int k = 0; k < 3; k++)
275 for(int l = 0; l < 3; l++) val += m(i, j, k, l) * t(l, k, j, i);
276 return val;
277 }
278
279 #endif
280