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 1 : RankThreeTensorTempl<T>::fillMethodEnum()
42 : {
43 1 : return MooseEnum("general plane_normal");
44 : }
45 :
46 : template <typename T>
47 494 : RankThreeTensorTempl<T>::RankThreeTensorTempl()
48 : {
49 : mooseAssert(N == 3, "RankThreeTensor is currently only tested for 3 dimensions.");
50 :
51 12320 : for (auto i : make_range(N3))
52 11880 : _vals[i] = 0;
53 440 : }
54 :
55 : template <typename T>
56 1 : RankThreeTensorTempl<T>::RankThreeTensorTempl(const InitMethod init)
57 : {
58 1 : switch (init)
59 : {
60 1 : case initNone:
61 1 : break;
62 :
63 0 : default:
64 0 : mooseError("Unknown RankThreeTensor initialization pattern.");
65 : }
66 1 : }
67 :
68 : template <typename T>
69 21 : RankThreeTensorTempl<T>::RankThreeTensorTempl(const std::vector<T> & input, FillMethod fill_method)
70 : {
71 21 : fillFromInputVector(input, fill_method);
72 18 : }
73 :
74 : template <typename T>
75 : void
76 66 : RankThreeTensorTempl<T>::zero()
77 : {
78 1848 : for (auto i : make_range(N3))
79 1782 : _vals[i] = 0;
80 66 : }
81 :
82 : template <typename T>
83 : RankThreeTensorTempl<T> &
84 1 : RankThreeTensorTempl<T>::operator=(const T & value)
85 : {
86 28 : for (auto i : make_range(N3))
87 27 : _vals[i] = value;
88 1 : return *this;
89 : }
90 :
91 : template <typename T>
92 : RankThreeTensorTempl<T> &
93 3 : RankThreeTensorTempl<T>::operator=(const RankThreeTensorTempl<T> & a)
94 : {
95 84 : for (auto i : make_range(N3))
96 81 : _vals[i] = a._vals[i];
97 :
98 3 : return *this;
99 : }
100 :
101 : template <typename T>
102 : libMesh::VectorValue<T>
103 1 : RankThreeTensorTempl<T>::operator*(const RankTwoTensorTempl<T> & a) const
104 : {
105 1 : libMesh::VectorValue<T> result;
106 :
107 4 : for (auto i : make_range(N))
108 : {
109 3 : T sum = 0;
110 3 : unsigned int i1 = i * N2;
111 12 : for (unsigned int j1 = 0; j1 < N2; j1 += N)
112 36 : for (auto k : make_range(N))
113 27 : sum += _vals[i1 + j1 + k] * a._coords[j1 + k];
114 3 : result(i) = sum;
115 : }
116 :
117 1 : 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 2 : RankThreeTensorTempl<T>::operator*(const T b) const
143 : {
144 2 : RankThreeTensorTempl<T> result;
145 :
146 56 : for (auto i : make_range(N3))
147 54 : result._vals[i] = _vals[i] * b;
148 :
149 2 : return result;
150 0 : }
151 :
152 : template <typename T>
153 : RankThreeTensorTempl<T> &
154 1 : RankThreeTensorTempl<T>::operator*=(const T a)
155 : {
156 28 : for (auto i : make_range(N3))
157 27 : _vals[i] *= a;
158 :
159 1 : return *this;
160 : }
161 :
162 : template <typename T>
163 : RankThreeTensorTempl<T>
164 1 : RankThreeTensorTempl<T>::operator/(const T b) const
165 : {
166 1 : RankThreeTensorTempl<T> result;
167 :
168 28 : for (auto i : make_range(N3))
169 27 : result._vals[i] = _vals[i] / b;
170 :
171 1 : return result;
172 0 : }
173 :
174 : template <typename T>
175 : RankThreeTensorTempl<T> &
176 1 : RankThreeTensorTempl<T>::operator/=(const T a)
177 : {
178 28 : for (auto i : make_range(N3))
179 27 : _vals[i] /= a;
180 :
181 1 : return *this;
182 : }
183 :
184 : template <typename T>
185 : RankThreeTensorTempl<T> &
186 1 : RankThreeTensorTempl<T>::operator+=(const RankThreeTensorTempl<T> & a)
187 : {
188 28 : for (auto i : make_range(N3))
189 27 : _vals[i] += a._vals[i];
190 :
191 1 : return *this;
192 : }
193 :
194 : template <typename T>
195 : RankThreeTensorTempl<T>
196 1 : RankThreeTensorTempl<T>::operator+(const RankThreeTensorTempl<T> & b) const
197 : {
198 1 : RankThreeTensorTempl<T> result;
199 :
200 28 : for (auto i : make_range(N3))
201 27 : result._vals[i] = _vals[i] + b._vals[i];
202 :
203 1 : return result;
204 0 : }
205 :
206 : template <typename T>
207 : RankThreeTensorTempl<T> &
208 1 : RankThreeTensorTempl<T>::operator-=(const RankThreeTensorTempl<T> & a)
209 : {
210 28 : for (auto i : make_range(N3))
211 27 : _vals[i] -= a._vals[i];
212 :
213 1 : return *this;
214 : }
215 :
216 : template <typename T>
217 : RankThreeTensorTempl<T>
218 1 : RankThreeTensorTempl<T>::operator-(const RankThreeTensorTempl<T> & b) const
219 : {
220 1 : RankThreeTensorTempl<T> result;
221 :
222 28 : for (auto i : make_range(N3))
223 27 : result._vals[i] = _vals[i] - b._vals[i];
224 :
225 1 : return result;
226 0 : }
227 :
228 : template <typename T>
229 : RankThreeTensorTempl<T>
230 1 : RankThreeTensorTempl<T>::operator-() const
231 : {
232 1 : RankThreeTensorTempl<T> result;
233 :
234 28 : for (auto i : make_range(N3))
235 27 : result._vals[i] = -_vals[i];
236 :
237 1 : return result;
238 0 : }
239 :
240 : template <typename T>
241 : T
242 1 : RankThreeTensorTempl<T>::L2norm() const
243 : {
244 1 : T l2 = 0;
245 :
246 28 : for (auto i : make_range(N3))
247 27 : l2 += Utility::pow<2>(_vals[i]);
248 :
249 1 : return std::sqrt(l2);
250 0 : }
251 :
252 : template <typename T>
253 : void
254 64 : RankThreeTensorTempl<T>::fillFromInputVector(const std::vector<T> & input, FillMethod fill_method)
255 : {
256 64 : zero();
257 :
258 64 : if (fill_method == automatic)
259 : {
260 20 : if (input.size() == 27)
261 18 : fill_method = general;
262 2 : else if (input.size() == 3)
263 1 : fill_method = plane_normal;
264 : else
265 1 : mooseError("Unsupported automatic fill method, use 27 values for 'general' and 3 for "
266 : "'plane_normal', the supplied size was ",
267 1 : input.size(),
268 : ".");
269 : }
270 :
271 63 : if (fill_method == general)
272 61 : fillGeneralFromInputVector(input);
273 :
274 2 : else if (fill_method == plane_normal)
275 : {
276 2 : if (input.size() != 3)
277 1 : mooseError("To use fillFromPlaneNormal, your input must have size 3, the supplied size was ",
278 1 : input.size(),
279 : ".");
280 1 : fillFromPlaneNormal(libMesh::VectorValue<T>(input[0], input[1], input[2]));
281 : }
282 :
283 : else
284 : // This is un-reachable unless a FillMethod is added and the if statement is not updated
285 0 : mooseError("fillFromInputVector called with unknown fill_method of ", fill_method);
286 61 : }
287 :
288 : template <typename T>
289 : void
290 2 : RankThreeTensorTempl<T>::fillFromPlaneNormal(const libMesh::VectorValue<T> & input)
291 : {
292 2 : unsigned int index = 0;
293 8 : for (auto i : make_range(N))
294 : {
295 6 : const T a = input(i);
296 24 : for (auto j : make_range(N))
297 : {
298 18 : const T b = input(j);
299 72 : for (auto k : make_range(N))
300 : {
301 54 : const T c = input(k);
302 54 : T sum = 0;
303 54 : sum = -2.0 * a * b * c;
304 54 : if (i == j)
305 18 : sum += c;
306 54 : if (i == k)
307 18 : sum += b;
308 54 : _vals[index++] = sum / 2.0;
309 : }
310 : }
311 : }
312 2 : }
313 :
314 : template <typename T>
315 : RankFourTensorTempl<T>
316 1 : RankThreeTensorTempl<T>::mixedProductRankFour(const RankTwoTensorTempl<T> & a) const
317 : {
318 1 : RankFourTensorTempl<T> result;
319 :
320 1 : unsigned int index = 0;
321 4 : for (auto i : make_range(N))
322 12 : for (auto j : make_range(N))
323 36 : for (auto k : make_range(N))
324 108 : for (auto l : make_range(N))
325 : {
326 324 : for (auto m : make_range(N))
327 972 : for (auto n : make_range(N))
328 729 : result._vals[index] += (*this)(m, i, j) * a(m, n) * (*this)(n, k, l);
329 81 : index++;
330 : }
331 :
332 1 : return result;
333 0 : }
334 :
335 : template <typename T>
336 : void
337 1 : RankThreeTensorTempl<T>::rotate(const libMesh::TensorValue<T> & R)
338 : {
339 1 : RankThreeTensorTempl<T> old = *this;
340 :
341 1 : unsigned int index = 0;
342 4 : for (auto i : make_range(N))
343 12 : for (auto j : make_range(N))
344 36 : for (auto k : make_range(N))
345 : {
346 27 : T sum = 0.0;
347 27 : unsigned int index2 = 0;
348 108 : for (auto m : make_range(N))
349 : {
350 81 : T a = R(i, m);
351 324 : for (auto n : make_range(N))
352 : {
353 243 : T ab = a * R(j, n);
354 972 : for (auto o : make_range(N))
355 729 : sum += ab * R(k, o) * old._vals[index2++];
356 : }
357 : }
358 27 : _vals[index++] = sum;
359 : }
360 1 : }
361 :
362 : template <typename T>
363 : void
364 61 : RankThreeTensorTempl<T>::fillGeneralFromInputVector(const std::vector<T> & input)
365 : {
366 61 : if (input.size() != 27)
367 1 : mooseError(
368 : "To use fillGeneralFromInputVector, your input must have size 27, the supplied size was ",
369 1 : input.size(),
370 : ".");
371 :
372 1680 : for (auto i : make_range(N3))
373 1620 : _vals[i] = input[i];
374 60 : }
375 :
376 : template <typename T>
377 : libMesh::VectorValue<T>
378 1 : RankThreeTensorTempl<T>::doubleContraction(const RankTwoTensorTempl<T> & b) const
379 : {
380 1 : libMesh::VectorValue<T> result;
381 :
382 4 : for (auto i : make_range(N))
383 30 : for (auto j : make_range(N2))
384 27 : result(i) += _vals[i * N2 + j] * b._coords[j];
385 :
386 1 : return result;
387 0 : }
388 :
389 : template <typename T>
390 : void
391 1 : RankThreeTensorTempl<T>::print(std::ostream & stm) const
392 : {
393 4 : for (auto i : make_range(N))
394 : {
395 3 : stm << "a(" << i << ", j, k) = \n";
396 12 : for (auto j : make_range(N))
397 : {
398 36 : for (auto k : make_range(N))
399 27 : stm << std::setw(15) << (*this)(i, j, k) << ' ';
400 9 : stm << "\n";
401 : }
402 3 : stm << "\n";
403 : }
404 :
405 1 : stm << std::flush;
406 1 : }
|
