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 "Moose.h"
13 : #include "ADRankTwoTensorForward.h"
14 : #include "ADRankThreeTensorForward.h"
15 : #include "ADRankFourTensorForward.h"
16 :
17 : #include "libmesh/libmesh.h"
18 : #include "libmesh/int_range.h"
19 :
20 : #include "metaphysicl/raw_type.h"
21 :
22 : namespace libMesh
23 : {
24 : template <typename>
25 : class TensorValue;
26 : template <typename>
27 : class TypeTensor;
28 : template <typename>
29 : class VectorValue;
30 : }
31 :
32 : // Forward declarations
33 : class MooseEnum;
34 :
35 : namespace MathUtils
36 : {
37 : template <typename T>
38 : void mooseSetToZero(T & v);
39 :
40 : /**
41 : * Helper function template specialization to set an object to zero.
42 : * Needed by DerivativeMaterialInterface
43 : */
44 : template <>
45 : void mooseSetToZero<RankThreeTensorTempl<Real>>(RankThreeTensorTempl<Real> & v);
46 : template <>
47 : void mooseSetToZero<RankThreeTensorTempl<ADReal>>(RankThreeTensorTempl<ADReal> & v);
48 : }
49 :
50 : /**
51 : * RankThreeTensor is designed to handle any N-dimensional third order tensor, r.
52 : *
53 : */
54 : template <typename T>
55 : class RankThreeTensorTempl
56 : {
57 : public:
58 : ///@{ tensor dimension and powers of the dimension
59 : static constexpr unsigned int N = Moose::dim;
60 : static constexpr unsigned int N2 = N * N;
61 : static constexpr unsigned int N3 = N * N * N;
62 : ///@}
63 :
64 : /// Initialization method
65 : enum InitMethod
66 : {
67 : initNone
68 : };
69 :
70 : /**
71 : * To fill up the 27 entries in the 3rd-order tensor, fillFromInputVector
72 : * is called with one of the following fill_methods.
73 : * See the fill*FromInputVector functions for more details
74 : */
75 : enum FillMethod
76 : {
77 : automatic,
78 : general,
79 : plane_normal
80 : };
81 :
82 : /// Default constructor; fills to zero
83 : RankThreeTensorTempl();
84 :
85 : /// Copy assignment operator must be defined if used
86 0 : RankThreeTensorTempl(const RankThreeTensorTempl<T> & a) = default;
87 :
88 : /**
89 : * Construct from other class template instantiation
90 : */
91 : template <typename T2>
92 : RankThreeTensorTempl(const RankThreeTensorTempl<T2> & copy);
93 :
94 : /// Select specific initialization pattern
95 : RankThreeTensorTempl(const InitMethod);
96 :
97 : /// Fill from vector
98 : RankThreeTensorTempl(const std::vector<T> &, FillMethod method = automatic);
99 :
100 : /// Gets the value for the index specified. Takes index = 0,1,2
101 113790 : inline T & operator()(unsigned int i, unsigned int j, unsigned int k)
102 : {
103 113790 : return _vals[((i * N + j) * N + k)];
104 : }
105 :
106 : /// Gets the value for the index specified. Takes index = 0,1,2. Used for const
107 1996 : inline T operator()(unsigned int i, unsigned int j, unsigned int k) const
108 : {
109 1996 : return _vals[((i * N + j) * N + k)];
110 : }
111 :
112 : /// Assignment-from-scalar operator.
113 : RankThreeTensorTempl<T> & operator=(const T & value);
114 :
115 : /// Zeros out the tensor.
116 : void zero();
117 :
118 : /// Print the rank three tensor
119 : void print(std::ostream & stm = Moose::out) const;
120 :
121 : friend std::ostream & operator<<(std::ostream & os, const RankThreeTensorTempl<T> & t)
122 : {
123 : t.print(os);
124 : return os;
125 : }
126 :
127 : /// copies values from "a" into this tensor
128 : RankThreeTensorTempl<T> & operator=(const RankThreeTensorTempl<T> & a);
129 :
130 : template <typename T2>
131 : RankThreeTensorTempl<T> & operator=(const RankThreeTensorTempl<T2> & a);
132 :
133 : /// b_i = r_ijk * a_jk
134 : libMesh::VectorValue<T> operator*(const RankTwoTensorTempl<T> & a) const;
135 :
136 : /// b_ij = r_ijk * a_k
137 : RankTwoTensorTempl<T> operator*(const libMesh::VectorValue<T> & a) const;
138 :
139 : /// r_ijk*a
140 : RankThreeTensorTempl<T> operator*(const T a) const;
141 :
142 : /// r_ijk *= a
143 : RankThreeTensorTempl<T> & operator*=(const T a);
144 :
145 : /// r_ijk/a
146 : RankThreeTensorTempl<T> operator/(const T a) const;
147 :
148 : /// r_ijk /= a for all i, j, k
149 : RankThreeTensorTempl<T> & operator/=(const T a);
150 :
151 : /// r_ijk += a_ijk for all i, j, k
152 : RankThreeTensorTempl<T> & operator+=(const RankThreeTensorTempl<T> & a);
153 :
154 : /// r_ijkl + a_ijk
155 : RankThreeTensorTempl<T> operator+(const RankThreeTensorTempl<T> & a) const;
156 :
157 : /// r_ijk -= a_ijk
158 : RankThreeTensorTempl<T> & operator-=(const RankThreeTensorTempl<T> & a);
159 :
160 : /// r_ijk - a_ijk
161 : RankThreeTensorTempl<T> operator-(const RankThreeTensorTempl<T> & a) const;
162 :
163 : /// -r_ijk
164 : RankThreeTensorTempl<T> operator-() const;
165 :
166 : /// \sqrt(r_ijk*r_ijk)
167 : T L2norm() const;
168 :
169 : /**
170 : * Rotate the tensor using
171 : * r_ijk = R_im R_in R_ko r_mno
172 : */
173 : template <class T2>
174 : void rotate(const T2 & R);
175 :
176 : /**
177 : * Rotate the tensor using
178 : * r_ijk = R_im R_in R_ko r_mno
179 : */
180 : void rotate(const libMesh::TensorValue<T> & R);
181 :
182 : /// Static method for use in validParams for getting the "fill_method"
183 : static MooseEnum fillMethodEnum();
184 :
185 : /**
186 : * fillFromInputVector takes some number of inputs to fill
187 : * the Rank-3 tensor.
188 : * @param input the numbers that will be placed in the tensor
189 : * @param fill_method this can be:
190 : * general (use fillGeneralFromInputVector)
191 : * more fill_methods to be implemented soon!
192 : */
193 : void fillFromInputVector(const std::vector<T> & input, FillMethod fill_method = automatic);
194 :
195 : /**
196 : * Fills RankThreeTensor from plane normal vectors
197 : * ref. Kuhl et. al. Int. J. Solids Struct. 38(2001) 2933-2952
198 : * @param input plane normal vector
199 : */
200 : void fillFromPlaneNormal(const libMesh::VectorValue<T> & input);
201 :
202 : /**
203 : * Creates fourth order tensor D_{ijkl}=A_{mij}*b_{mn}*A_{nkl} where A is rank 3 and b is rank 2
204 : * @param a RankThreeTensor A in the equation above
205 : */
206 : RankFourTensorTempl<T> mixedProductRankFour(const RankTwoTensorTempl<T> & a) const;
207 :
208 : /**
209 : * Creates a vector from the double contraction of a rank three and rank two tensor.
210 : *
211 : * c_i = A_{ijk} * b_{jk}
212 : */
213 : libMesh::VectorValue<T> doubleContraction(const RankTwoTensorTempl<T> & b) const;
214 :
215 : protected:
216 : /// The values of the rank-three tensor stored by index=((i * LIBMESH_DIM + j) * LIBMESH_DIM + k)
217 : T _vals[N3];
218 :
219 : void fillGeneralFromInputVector(const std::vector<T> & input);
220 :
221 : template <class T2>
222 : friend void dataStore(std::ostream &, RankThreeTensorTempl<T2> &, void *);
223 :
224 : template <class T2>
225 : friend void dataLoad(std::istream &, RankThreeTensorTempl<T2> &, void *);
226 :
227 : template <class T2>
228 : friend class RankTwoTensorTempl;
229 :
230 : template <class T2>
231 : friend class RankThreeTensorTempl;
232 :
233 : template <class T2>
234 : friend class RankFourTensorTempl;
235 : };
236 :
237 : namespace MetaPhysicL
238 : {
239 : template <typename T>
240 : struct RawType<RankThreeTensorTempl<T>>
241 : {
242 : typedef RankThreeTensorTempl<typename RawType<T>::value_type> value_type;
243 :
244 2 : static value_type value(const RankThreeTensorTempl<T> & in)
245 : {
246 2 : value_type ret;
247 8 : for (const auto i : libMesh::make_range(RankThreeTensorTempl<T>::N))
248 24 : for (const auto j : libMesh::make_range(RankThreeTensorTempl<T>::N))
249 72 : for (const auto k : libMesh::make_range(RankThreeTensorTempl<T>::N))
250 54 : ret(i, j, k) = raw_value(in(i, j, k));
251 :
252 2 : return ret;
253 : }
254 : };
255 : }
256 :
257 : template <typename T>
258 : template <typename T2>
259 : RankThreeTensorTempl<T>::RankThreeTensorTempl(const RankThreeTensorTempl<T2> & copy)
260 : {
261 : for (const auto i : libMesh::make_range(N3))
262 : _vals[i] = copy._vals[i];
263 : }
264 :
265 : template <typename T>
266 : template <class T2>
267 : void
268 1 : RankThreeTensorTempl<T>::rotate(const T2 & R)
269 : {
270 1 : unsigned int index = 0;
271 4 : for (const auto i : libMesh::make_range(N))
272 12 : for (const auto j : libMesh::make_range(N))
273 36 : for (const auto k : libMesh::make_range(N))
274 : {
275 27 : unsigned int index2 = 0;
276 27 : T sum = 0.0;
277 108 : for (const auto m : libMesh::make_range(N))
278 : {
279 81 : T a = R(i, m);
280 324 : for (const auto n : libMesh::make_range(N))
281 : {
282 243 : T ab = a * R(j, n);
283 972 : for (const auto o : libMesh::make_range(N))
284 729 : sum += ab * R(k, o) * _vals[index2++];
285 : }
286 : }
287 27 : _vals[index++] = sum;
288 : }
289 1 : }
290 :
291 : template <typename T>
292 : RankThreeTensorTempl<T>
293 1 : operator*(T a, const RankThreeTensorTempl<T> & b)
294 : {
295 1 : return b * a;
296 : }
297 :
298 : ///r=v*A where r is rank 2, v is vector and A is rank 3
299 : template <typename T>
300 : RankTwoTensorTempl<T>
301 1 : operator*(const libMesh::VectorValue<T> & p, const RankThreeTensorTempl<T> & b)
302 : {
303 : static_assert(RankThreeTensorTempl<T>::N == RankTwoTensorTempl<T>::N,
304 : "RankTwoTensor and RankThreeTensor have to have the same dimension N.");
305 1 : RankTwoTensorTempl<T> result;
306 :
307 4 : for (const auto i : libMesh::make_range(RankThreeTensorTempl<T>::N))
308 12 : for (const auto j : libMesh::make_range(RankThreeTensorTempl<T>::N))
309 36 : for (const auto k : libMesh::make_range(RankThreeTensorTempl<T>::N))
310 27 : result(i, j) += p(k) * b(k, i, j);
311 :
312 1 : return result;
313 0 : }
314 :
315 : template <typename T>
316 : template <typename T2>
317 : RankThreeTensorTempl<T> &
318 2 : RankThreeTensorTempl<T>::operator=(const RankThreeTensorTempl<T2> & a)
319 : {
320 8 : for (const auto i : libMesh::make_range(N))
321 24 : for (const auto j : libMesh::make_range(N))
322 72 : for (const auto k : libMesh::make_range(N))
323 54 : (*this)(i, j, k) = a(i, j, k);
324 :
325 2 : return *this;
326 : }
|
