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 "RankTwoTensor.h"
13 : #include "SymmetricRankTwoTensor.h"
14 :
15 : // forward declarations
16 : template <typename>
17 : class FactorizedRankTwoTensorTempl;
18 :
19 : /**
20 : * FactorizedRankTwoTensorTempl is designed to perform the spectral decomposition of an
21 : * underlying symmetric second order tensor and reuse its bases for future operations if possible.
22 : *
23 : * FactorizedRankTwoTensorTempl templates on the underlying second order tensor.
24 : * IMPORTANT: the underlying second order tensor must be symmetric. A check is only performed in
25 : * debug mode.
26 : *
27 : * Only operations that reuses the known factorization are provided. Otherwise, you will need to
28 : * first retrieve the underlying RankTwoTensorTempl to perform the operation.
29 : *
30 : * TODO? Although I am calling it factorization, it only really refers to eigenvalue decompositon at
31 : * this point. I am not sure if in the future we need to add other types of similarity
32 : * transformation.
33 : */
34 : template <typename T>
35 : class FactorizedRankTwoTensorTempl
36 : {
37 : public:
38 : /// For generic programming
39 : typedef typename T::value_type value_type;
40 :
41 : /// No default constructor
42 : FactorizedRankTwoTensorTempl() = delete;
43 :
44 : /// Copy constructor
45 6 : FactorizedRankTwoTensorTempl(const FactorizedRankTwoTensorTempl<T> & A) = default;
46 :
47 : /// Constructor if the factorization isn't known a priori
48 23 : FactorizedRankTwoTensorTempl(const T & A)
49 23 : {
50 23 : if (!A.isSymmetric())
51 0 : mooseError("The tensor is not symmetric.");
52 23 : A.symmetricEigenvaluesEigenvectors(_eigvals, _eigvecs);
53 23 : }
54 :
55 : // Construct from the factorization
56 : // Assume that regardless of the type of T, the eigenvectors are always stored as as a full second
57 : // order tensor, e.g. the eigenvector matrix of a symmetric second order tensor isn't symmetric in
58 : // general. Hence we don't take advantage of orthonormal and/or unitary second order tensors.
59 8 : FactorizedRankTwoTensorTempl(const std::vector<typename T::value_type> & eigvals,
60 : const RankTwoTensorTempl<typename T::value_type> & eigvecs)
61 8 : : _eigvals(eigvals), _eigvecs(eigvecs)
62 : {
63 8 : }
64 :
65 : // @{ getters
66 : template <typename T2 = T>
67 4 : T2 get() const
68 : {
69 4 : return static_cast<T2>(assemble());
70 : }
71 : template <typename T2 = T>
72 180 : T2 get()
73 : {
74 180 : return static_cast<T2>(assemble());
75 : }
76 70 : const std::vector<typename T::value_type> & eigvals() const { return _eigvals; }
77 36 : std::vector<typename T::value_type> eigvals() { return _eigvals; }
78 50 : const RankTwoTensorTempl<typename T::value_type> & eigvecs() const { return _eigvecs; }
79 63 : RankTwoTensorTempl<typename T::value_type> eigvecs() { return _eigvecs; }
80 : // @}
81 :
82 : void print(std::ostream & stm = Moose::out) const;
83 :
84 : // Returns _A rotated by R.
85 : FactorizedRankTwoTensorTempl<T>
86 : rotated(const RankTwoTensorTempl<typename T::value_type> & R) const;
87 :
88 : /// Returns the transpose of _A.
89 : FactorizedRankTwoTensorTempl<T> transpose() const;
90 :
91 : // @{ Assignment operators
92 : FactorizedRankTwoTensorTempl<T> & operator=(const FactorizedRankTwoTensorTempl<T> & A);
93 : FactorizedRankTwoTensorTempl<T> & operator=(const T & A);
94 : // @}
95 :
96 : /// performs _A *= a in place, also updates eigen values
97 : FactorizedRankTwoTensorTempl<T> & operator*=(const typename T::value_type & a);
98 :
99 : /// returns _A * a, also updates eigen values
100 : template <typename T2>
101 : FactorizedRankTwoTensorTempl<T> operator*(const T2 & a) const;
102 :
103 : /// performs _A /= a in place, also updates eigen values
104 : FactorizedRankTwoTensorTempl<T> & operator/=(const typename T::value_type & a);
105 :
106 : /// returns _A / a, also updates eigen values
107 : template <typename T2>
108 : FactorizedRankTwoTensorTempl<T> operator/(const T2 & a) const;
109 :
110 : /// Defines logical equality with another second order tensor
111 : bool operator==(const T & A) const;
112 :
113 : /// Defines logical equality with another FactorizedRankTwoTensorTempl<T>
114 : bool operator==(const FactorizedRankTwoTensorTempl<T> & A) const;
115 :
116 : /// inverse of _A
117 : FactorizedRankTwoTensorTempl<T> inverse() const;
118 :
119 : /// add identity times a to _A
120 : void addIa(const typename T::value_type & a);
121 :
122 : /// trace of _A
123 : typename T::value_type trace() const;
124 :
125 : /// determinant of _A
126 : typename T::value_type det() const;
127 :
128 : private:
129 : // Assemble the tensor from the factorization
130 184 : RankTwoTensorTempl<typename T::value_type> assemble() const
131 : {
132 : typedef RankTwoTensorTempl<typename T::value_type> R2T;
133 368 : return _eigvals[0] * R2T::selfOuterProduct(_eigvecs.column(0)) +
134 368 : _eigvals[1] * R2T::selfOuterProduct(_eigvecs.column(1)) +
135 920 : _eigvals[2] * R2T::selfOuterProduct(_eigvecs.column(2));
136 : }
137 :
138 : // The eigen values of _A;
139 : std::vector<typename T::value_type> _eigvals;
140 :
141 : // The eigen vectors of _A;
142 : RankTwoTensorTempl<typename T::value_type> _eigvecs;
143 : };
144 :
145 : namespace MathUtils
146 : {
147 : #define FactorizedRankTwoTensorOperatorMapBody(operator) \
148 : { \
149 : std::vector<typename T::value_type> op_eigvals; \
150 : for (const auto & eigval : A.eigvals()) \
151 : op_eigvals.push_back(operator); \
152 : return FactorizedRankTwoTensorTempl<T>(op_eigvals, A.eigvecs()); \
153 : }
154 :
155 : #define FactorizedRankTwoTensorOperatorMapUnary(operatorname, operator) \
156 : template <typename T> \
157 : FactorizedRankTwoTensorTempl<T> operatorname(const FactorizedRankTwoTensorTempl<T> & A) \
158 : { \
159 : FactorizedRankTwoTensorOperatorMapBody(operator) \
160 : } \
161 : static_assert(true, "")
162 :
163 : #define FactorizedRankTwoTensorOperatorMapBinary(operatorname, operator) \
164 : template <typename T, typename T2> \
165 : FactorizedRankTwoTensorTempl<T> operatorname(const FactorizedRankTwoTensorTempl<T> & A, \
166 : const T2 & arg) \
167 : { \
168 : if constexpr (libMesh::ScalarTraits<T2>::value) \
169 : { \
170 : FactorizedRankTwoTensorOperatorMapBody(operator) \
171 : } \
172 : } \
173 : static_assert(true, "")
174 :
175 : #define FactorizedRankTwoTensorOperatorMapDerivativeBody(operator, derivative) \
176 : { \
177 : std::vector<typename T::value_type> op_eigvals, op_derivs; \
178 : for (const auto & eigval : A.eigvals()) \
179 : { \
180 : op_eigvals.push_back(operator); \
181 : op_derivs.push_back(derivative); \
182 : } \
183 : \
184 : RankFourTensorTempl<typename T::value_type> P, Gab, Gba; \
185 : typedef RankTwoTensorTempl<typename T::value_type> R2T; \
186 : \
187 : for (auto a : make_range(3)) \
188 : { \
189 : const auto Ma = R2T::selfOuterProduct(A.eigvecs().column(a)); \
190 : P += op_derivs[a] * Ma.outerProduct(Ma); \
191 : } \
192 : \
193 : for (auto a : make_range(3)) \
194 : for (auto b : make_range(a)) \
195 : { \
196 : const auto Ma = R2T::selfOuterProduct(A.eigvecs().column(a)); \
197 : const auto Mb = R2T::selfOuterProduct(A.eigvecs().column(b)); \
198 : \
199 : usingTensorIndices(i_, j_, k_, l_); \
200 : Gab = Ma.template times<i_, k_, j_, l_>(Mb) + Ma.template times<i_, l_, j_, k_>(Mb); \
201 : Gba = Mb.template times<i_, k_, j_, l_>(Ma) + Mb.template times<i_, l_, j_, k_>(Ma); \
202 : \
203 : Real theta_ab; \
204 : if (!MooseUtils::absoluteFuzzyEqual(A.eigvals()[a], A.eigvals()[b])) \
205 : theta_ab = 0.5 * (op_eigvals[a] - op_eigvals[b]) / (A.eigvals()[a] - A.eigvals()[b]); \
206 : else \
207 : theta_ab = 0.25 * (op_derivs[a] + op_derivs[b]); \
208 : \
209 : P += theta_ab * (Gab + Gba); \
210 : } \
211 : return P; \
212 : }
213 :
214 : #define FactorizedRankTwoTensorOperatorMapDerivativeUnary(derivativename, operator, derivative) \
215 : template <typename T> \
216 : RankFourTensorTempl<typename T::value_type> derivativename( \
217 : const FactorizedRankTwoTensorTempl<T> & A) \
218 : { \
219 : FactorizedRankTwoTensorOperatorMapDerivativeBody(operator, derivative) \
220 : } \
221 : static_assert(true, "")
222 :
223 : #define FactorizedRankTwoTensorOperatorMapDerivativeBinary(derivativename, operator, derivative) \
224 : template <typename T, typename T2> \
225 : RankFourTensorTempl<typename T::value_type> derivativename( \
226 : const FactorizedRankTwoTensorTempl<T> & A, const T2 & arg) \
227 : { \
228 : if constexpr (libMesh::ScalarTraits<T2>::value) \
229 : { \
230 : FactorizedRankTwoTensorOperatorMapDerivativeBody(operator, derivative) \
231 : } \
232 : } \
233 : static_assert(true, "")
234 :
235 : // TODO: While the macros are here, in the future we could instantiate other operator maps like
236 : // trignometry functions.
237 : // @{
238 6 : FactorizedRankTwoTensorOperatorMapUnary(log, std::log(eigval));
239 14 : FactorizedRankTwoTensorOperatorMapDerivativeUnary(dlog, std::log(eigval), 1 / eigval);
240 :
241 6 : FactorizedRankTwoTensorOperatorMapUnary(exp, std::exp(eigval));
242 14 : FactorizedRankTwoTensorOperatorMapDerivativeUnary(dexp, std::exp(eigval), std::exp(eigval));
243 :
244 6 : FactorizedRankTwoTensorOperatorMapUnary(sqrt, std::sqrt(eigval));
245 14 : FactorizedRankTwoTensorOperatorMapDerivativeUnary(dsqrt,
246 : std::sqrt(eigval),
247 : std::pow(eigval, -1. / 2.) / 2.);
248 :
249 6 : FactorizedRankTwoTensorOperatorMapUnary(cbrt, std::cbrt(eigval));
250 14 : FactorizedRankTwoTensorOperatorMapDerivativeUnary(dcbrt,
251 : std::cbrt(eigval),
252 : std::pow(eigval, -2. / 3.) / 3.);
253 :
254 6 : FactorizedRankTwoTensorOperatorMapBinary(pow, std::pow(eigval, arg));
255 14 : FactorizedRankTwoTensorOperatorMapDerivativeBinary(dpow,
256 : std::pow(eigval, arg),
257 : arg * std::pow(eigval, arg - 1));
258 : // @}
259 : } // end namespace MathUtils
260 :
261 : template <typename T>
262 : template <typename T2>
263 : FactorizedRankTwoTensorTempl<T>
264 1 : FactorizedRankTwoTensorTempl<T>::operator*(const T2 & a) const
265 : {
266 : if constexpr (libMesh::ScalarTraits<T2>::value)
267 : {
268 1 : FactorizedRankTwoTensorTempl<T> A = *this;
269 4 : for (auto & eigval : A._eigvals)
270 3 : eigval *= a;
271 2 : return A;
272 1 : }
273 : }
274 :
275 : template <typename T>
276 : template <typename T2>
277 : FactorizedRankTwoTensorTempl<T>
278 1 : FactorizedRankTwoTensorTempl<T>::operator/(const T2 & a) const
279 : {
280 : if constexpr (libMesh::ScalarTraits<T2>::value)
281 : {
282 1 : FactorizedRankTwoTensorTempl<T> A = *this;
283 4 : for (auto & eigval : A._eigvals)
284 3 : eigval /= a;
285 2 : return A;
286 1 : }
287 : }
288 :
289 : typedef FactorizedRankTwoTensorTempl<RankTwoTensor> FactorizedRankTwoTensor;
290 : typedef FactorizedRankTwoTensorTempl<ADRankTwoTensor> ADFactorizedRankTwoTensor;
291 : typedef FactorizedRankTwoTensorTempl<SymmetricRankTwoTensor> FactorizedSymmetricRankTwoTensor;
292 : typedef FactorizedRankTwoTensorTempl<ADSymmetricRankTwoTensor> ADFactorizedSymmetricRankTwoTensor;
|
