Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "RankThreeTensor.h"
13 : #include "RankTwoTensor.h"
14 : #include "RankFourTensor.h"
15 :
16 : // MOOSE includes
17 : #include "MooseEnum.h"
18 : #include "MooseException.h"
19 : #include "MooseUtils.h"
20 : #include "MatrixTools.h"
21 : #include "PermutationTensor.h"
22 :
23 : #include "libmesh/utility.h"
24 : #include "libmesh/vector_value.h"
25 : #include "libmesh/tensor_value.h"
26 :
27 : // C++ includes
28 : #include <iomanip>
29 : #include <ostream>
30 :
31 : namespace MathUtils
32 : {
33 : template <>
34 : void mooseSetToZero<RankThreeTensor>(RankThreeTensor & v);
35 : template <>
36 : void mooseSetToZero<ADRankThreeTensor>(ADRankThreeTensor & v);
37 : }
38 :
39 : template <typename T>
40 : MooseEnum
41 2 : RankThreeTensorTempl<T>::fillMethodEnum()
42 : {
43 8 : return MooseEnum("general plane_normal");
44 : }
45 :
46 : template <typename T>
47 640 : RankThreeTensorTempl<T>::RankThreeTensorTempl()
48 : {
49 : mooseAssert(N == 3, "RankThreeTensor is currently only tested for 3 dimensions.");
50 :
51 14896 : for (auto i : make_range(N3))
52 14364 : _vals[i] = 0;
53 532 : }
54 :
55 : template <typename T>
56 2 : RankThreeTensorTempl<T>::RankThreeTensorTempl(const InitMethod init)
57 : {
58 2 : switch (init)
59 : {
60 2 : case initNone:
61 2 : break;
62 :
63 0 : default:
64 0 : mooseError("Unknown RankThreeTensor initialization pattern.");
65 : }
66 2 : }
67 :
68 : template <typename T>
69 42 : RankThreeTensorTempl<T>::RankThreeTensorTempl(const std::vector<T> & input, FillMethod fill_method)
70 : {
71 42 : fillFromInputVector(input, fill_method);
72 36 : }
73 :
74 : template <typename T>
75 : void
76 132 : RankThreeTensorTempl<T>::zero()
77 : {
78 3696 : for (auto i : make_range(N3))
79 3564 : _vals[i] = 0;
80 132 : }
81 :
82 : template <typename T>
83 : RankThreeTensorTempl<T> &
84 2 : RankThreeTensorTempl<T>::operator=(const T & value)
85 : {
86 56 : for (auto i : make_range(N3))
87 54 : _vals[i] = value;
88 2 : return *this;
89 : }
90 :
91 : template <typename T>
92 : RankThreeTensorTempl<T> &
93 6 : RankThreeTensorTempl<T>::operator=(const RankThreeTensorTempl<T> & a)
94 : {
95 168 : for (auto i : make_range(N3))
96 162 : _vals[i] = a._vals[i];
97 :
98 6 : return *this;
99 : }
100 :
101 : template <typename T>
102 : libMesh::VectorValue<T>
103 2 : RankThreeTensorTempl<T>::operator*(const RankTwoTensorTempl<T> & a) const
104 : {
105 2 : libMesh::VectorValue<T> result;
106 :
107 8 : for (auto i : make_range(N))
108 : {
109 6 : T sum = 0;
110 6 : unsigned int i1 = i * N2;
111 24 : for (unsigned int j1 = 0; j1 < N2; j1 += N)
112 72 : for (auto k : make_range(N))
113 54 : sum += _vals[i1 + j1 + k] * a._coords[j1 + k];
114 6 : result(i) = sum;
115 : }
116 :
117 2 : return result;
118 0 : }
119 :
120 : template <typename T>
121 : RankTwoTensorTempl<T>
122 0 : RankThreeTensorTempl<T>::operator*(const libMesh::VectorValue<T> & a) const
123 : {
124 0 : RankTwoTensorTempl<T> result;
125 :
126 0 : for (auto i : make_range(N))
127 0 : for (auto j : make_range(N))
128 : {
129 0 : T sum = 0;
130 0 : unsigned int i1 = i * N2;
131 0 : unsigned int j1 = j * N;
132 0 : for (auto k : make_range(N))
133 0 : sum += _vals[i1 + j1 + k] * a(k);
134 0 : result(i, j) = sum;
135 : }
136 :
137 0 : return result;
138 0 : }
139 :
140 : template <typename T>
141 : RankThreeTensorTempl<T>
142 4 : RankThreeTensorTempl<T>::operator*(const T b) const
143 : {
144 4 : RankThreeTensorTempl<T> result;
145 :
146 112 : for (auto i : make_range(N3))
147 108 : result._vals[i] = _vals[i] * b;
148 :
149 4 : return result;
150 0 : }
151 :
152 : template <typename T>
153 : RankThreeTensorTempl<T> &
154 2 : RankThreeTensorTempl<T>::operator*=(const T a)
155 : {
156 56 : for (auto i : make_range(N3))
157 54 : _vals[i] *= a;
158 :
159 2 : return *this;
160 : }
161 :
162 : template <typename T>
163 : RankThreeTensorTempl<T>
164 2 : RankThreeTensorTempl<T>::operator/(const T b) const
165 : {
166 2 : RankThreeTensorTempl<T> result;
167 :
168 56 : for (auto i : make_range(N3))
169 54 : result._vals[i] = _vals[i] / b;
170 :
171 2 : return result;
172 0 : }
173 :
174 : template <typename T>
175 : RankThreeTensorTempl<T> &
176 2 : RankThreeTensorTempl<T>::operator/=(const T a)
177 : {
178 56 : for (auto i : make_range(N3))
179 54 : _vals[i] /= a;
180 :
181 2 : return *this;
182 : }
183 :
184 : template <typename T>
185 : RankThreeTensorTempl<T> &
186 2 : RankThreeTensorTempl<T>::operator+=(const RankThreeTensorTempl<T> & a)
187 : {
188 56 : for (auto i : make_range(N3))
189 54 : _vals[i] += a._vals[i];
190 :
191 2 : return *this;
192 : }
193 :
194 : template <typename T>
195 : RankThreeTensorTempl<T>
196 2 : RankThreeTensorTempl<T>::operator+(const RankThreeTensorTempl<T> & b) const
197 : {
198 2 : RankThreeTensorTempl<T> result;
199 :
200 56 : for (auto i : make_range(N3))
201 54 : result._vals[i] = _vals[i] + b._vals[i];
202 :
203 2 : return result;
204 0 : }
205 :
206 : template <typename T>
207 : RankThreeTensorTempl<T> &
208 2 : RankThreeTensorTempl<T>::operator-=(const RankThreeTensorTempl<T> & a)
209 : {
210 56 : for (auto i : make_range(N3))
211 54 : _vals[i] -= a._vals[i];
212 :
213 2 : return *this;
214 : }
215 :
216 : template <typename T>
217 : RankThreeTensorTempl<T>
218 2 : RankThreeTensorTempl<T>::operator-(const RankThreeTensorTempl<T> & b) const
219 : {
220 2 : RankThreeTensorTempl<T> result;
221 :
222 56 : for (auto i : make_range(N3))
223 54 : result._vals[i] = _vals[i] - b._vals[i];
224 :
225 2 : return result;
226 0 : }
227 :
228 : template <typename T>
229 : RankThreeTensorTempl<T>
230 2 : RankThreeTensorTempl<T>::operator-() const
231 : {
232 2 : RankThreeTensorTempl<T> result;
233 :
234 56 : for (auto i : make_range(N3))
235 54 : result._vals[i] = -_vals[i];
236 :
237 2 : return result;
238 0 : }
239 :
240 : template <typename T>
241 : T
242 2 : RankThreeTensorTempl<T>::L2norm() const
243 : {
244 2 : T l2 = 0;
245 :
246 56 : for (auto i : make_range(N3))
247 54 : l2 += Utility::pow<2>(_vals[i]);
248 :
249 : using std::sqrt;
250 2 : return sqrt(l2);
251 0 : }
252 :
253 : template <typename T>
254 : void
255 128 : RankThreeTensorTempl<T>::fillFromInputVector(const std::vector<T> & input, FillMethod fill_method)
256 : {
257 128 : zero();
258 :
259 128 : if (fill_method == automatic)
260 : {
261 40 : if (input.size() == 27)
262 36 : fill_method = general;
263 4 : else if (input.size() == 3)
264 2 : fill_method = plane_normal;
265 : else
266 2 : mooseError("Unsupported automatic fill method, use 27 values for 'general' and 3 for "
267 : "'plane_normal', the supplied size was ",
268 2 : input.size(),
269 : ".");
270 : }
271 :
272 126 : if (fill_method == general)
273 122 : fillGeneralFromInputVector(input);
274 :
275 4 : else if (fill_method == plane_normal)
276 : {
277 4 : if (input.size() != 3)
278 2 : mooseError("To use fillFromPlaneNormal, your input must have size 3, the supplied size was ",
279 2 : input.size(),
280 : ".");
281 2 : fillFromPlaneNormal(libMesh::VectorValue<T>(input[0], input[1], input[2]));
282 : }
283 :
284 : else
285 : // This is un-reachable unless a FillMethod is added and the if statement is not updated
286 0 : mooseError("fillFromInputVector called with unknown fill_method of ", fill_method);
287 122 : }
288 :
289 : template <typename T>
290 : void
291 4 : RankThreeTensorTempl<T>::fillFromPlaneNormal(const libMesh::VectorValue<T> & input)
292 : {
293 4 : unsigned int index = 0;
294 16 : for (auto i : make_range(N))
295 : {
296 12 : const T a = input(i);
297 48 : for (auto j : make_range(N))
298 : {
299 36 : const T b = input(j);
300 144 : for (auto k : make_range(N))
301 : {
302 108 : const T c = input(k);
303 108 : T sum = 0;
304 108 : sum = -2.0 * a * b * c;
305 108 : if (i == j)
306 36 : sum += c;
307 108 : if (i == k)
308 36 : sum += b;
309 108 : _vals[index++] = sum / 2.0;
310 : }
311 : }
312 : }
313 4 : }
314 :
315 : template <typename T>
316 : RankFourTensorTempl<T>
317 2 : RankThreeTensorTempl<T>::mixedProductRankFour(const RankTwoTensorTempl<T> & a) const
318 : {
319 2 : RankFourTensorTempl<T> result;
320 :
321 2 : unsigned int index = 0;
322 8 : for (auto i : make_range(N))
323 24 : for (auto j : make_range(N))
324 72 : for (auto k : make_range(N))
325 216 : for (auto l : make_range(N))
326 : {
327 648 : for (auto m : make_range(N))
328 1944 : for (auto n : make_range(N))
329 1458 : result._vals[index] += (*this)(m, i, j) * a(m, n) * (*this)(n, k, l);
330 162 : index++;
331 : }
332 :
333 2 : return result;
334 0 : }
335 :
336 : template <typename T>
337 : void
338 2 : RankThreeTensorTempl<T>::rotate(const libMesh::TensorValue<T> & R)
339 : {
340 2 : RankThreeTensorTempl<T> old = *this;
341 :
342 2 : unsigned int index = 0;
343 8 : for (auto i : make_range(N))
344 24 : for (auto j : make_range(N))
345 72 : for (auto k : make_range(N))
346 : {
347 54 : T sum = 0.0;
348 54 : unsigned int index2 = 0;
349 216 : for (auto m : make_range(N))
350 : {
351 162 : T a = R(i, m);
352 648 : for (auto n : make_range(N))
353 : {
354 486 : T ab = a * R(j, n);
355 1944 : for (auto o : make_range(N))
356 1458 : sum += ab * R(k, o) * old._vals[index2++];
357 : }
358 : }
359 54 : _vals[index++] = sum;
360 : }
361 2 : }
362 :
363 : template <typename T>
364 : void
365 122 : RankThreeTensorTempl<T>::fillGeneralFromInputVector(const std::vector<T> & input)
366 : {
367 122 : if (input.size() != 27)
368 2 : mooseError(
369 : "To use fillGeneralFromInputVector, your input must have size 27, the supplied size was ",
370 2 : input.size(),
371 : ".");
372 :
373 3360 : for (auto i : make_range(N3))
374 3240 : _vals[i] = input[i];
375 120 : }
376 :
377 : template <typename T>
378 : libMesh::VectorValue<T>
379 2 : RankThreeTensorTempl<T>::doubleContraction(const RankTwoTensorTempl<T> & b) const
380 : {
381 2 : libMesh::VectorValue<T> result;
382 :
383 8 : for (auto i : make_range(N))
384 60 : for (auto j : make_range(N2))
385 54 : result(i) += _vals[i * N2 + j] * b._coords[j];
386 :
387 2 : return result;
388 0 : }
389 :
390 : template <typename T>
391 : void
392 2 : RankThreeTensorTempl<T>::print(std::ostream & stm) const
393 : {
394 8 : for (auto i : make_range(N))
395 : {
396 6 : stm << "a(" << i << ", j, k) = \n";
397 24 : for (auto j : make_range(N))
398 : {
399 72 : for (auto k : make_range(N))
400 54 : stm << std::setw(15) << (*this)(i, j, k) << ' ';
401 18 : stm << "\n";
402 : }
403 6 : stm << "\n";
404 : }
405 :
406 2 : stm << std::flush;
407 2 : }
|
