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 12 : FactorizedRankTwoTensorTempl(const FactorizedRankTwoTensorTempl<T> & A) = default;
46 :
47 : /// Constructor if the factorization isn't known a priori
48 46 : FactorizedRankTwoTensorTempl(const T & A)
49 46 : {
50 46 : if (!A.isSymmetric())
51 0 : mooseError("The tensor is not symmetric.");
52 46 : A.symmetricEigenvaluesEigenvectors(_eigvals, _eigvecs);
53 46 : }
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 16 : FactorizedRankTwoTensorTempl(const std::vector<typename T::value_type> & eigvals,
60 : const RankTwoTensorTempl<typename T::value_type> & eigvecs)
61 16 : : _eigvals(eigvals), _eigvecs(eigvecs)
62 : {
63 16 : }
64 :
65 : // @{ getters
66 : template <typename T2 = T>
67 8 : T2 get() const
68 : {
69 8 : return static_cast<T2>(assemble());
70 : }
71 : template <typename T2 = T>
72 360 : T2 get()
73 : {
74 360 : return static_cast<T2>(assemble());
75 : }
76 140 : const std::vector<typename T::value_type> & eigvals() const { return _eigvals; }
77 72 : std::vector<typename T::value_type> eigvals() { return _eigvals; }
78 100 : const RankTwoTensorTempl<typename T::value_type> & eigvecs() const { return _eigvecs; }
79 126 : 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 368 : RankTwoTensorTempl<typename T::value_type> assemble() const
131 : {
132 : typedef RankTwoTensorTempl<typename T::value_type> R2T;
133 736 : return _eigvals[0] * R2T::selfOuterProduct(_eigvecs.column(0)) +
134 736 : _eigvals[1] * R2T::selfOuterProduct(_eigvecs.column(1)) +
135 1840 : _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 : using std::log, std::exp, std::sqrt, std::cbrt, std::pow; \
150 : std::vector<typename T::value_type> op_eigvals; \
151 : for (const auto & eigval : A.eigvals()) \
152 : op_eigvals.push_back(operator); \
153 : return FactorizedRankTwoTensorTempl<T>(op_eigvals, A.eigvecs()); \
154 : }
155 :
156 : #define FactorizedRankTwoTensorOperatorMapUnary(operatorname, operator) \
157 : template <typename T> \
158 : FactorizedRankTwoTensorTempl<T> operatorname(const FactorizedRankTwoTensorTempl<T> & A) \
159 : { \
160 : FactorizedRankTwoTensorOperatorMapBody(operator) \
161 : } \
162 : static_assert(true, "")
163 :
164 : #define FactorizedRankTwoTensorOperatorMapBinary(operatorname, operator) \
165 : template <typename T, typename T2> \
166 : FactorizedRankTwoTensorTempl<T> operatorname(const FactorizedRankTwoTensorTempl<T> & A, \
167 : const T2 & arg) \
168 : { \
169 : if constexpr (libMesh::ScalarTraits<T2>::value) \
170 : { \
171 : FactorizedRankTwoTensorOperatorMapBody(operator) \
172 : } \
173 : } \
174 : static_assert(true, "")
175 :
176 : #define FactorizedRankTwoTensorOperatorMapDerivativeBody(operator, derivative) \
177 : { \
178 : using std::log, std::exp, std::sqrt, std::cbrt, std::pow; \
179 : std::vector<typename T::value_type> op_eigvals, op_derivs; \
180 : for (const auto & eigval : A.eigvals()) \
181 : { \
182 : op_eigvals.push_back(operator); \
183 : op_derivs.push_back(derivative); \
184 : } \
185 : \
186 : RankFourTensorTempl<typename T::value_type> P, Gab, Gba; \
187 : typedef RankTwoTensorTempl<typename T::value_type> R2T; \
188 : \
189 : for (auto a : make_range(3)) \
190 : { \
191 : const auto Ma = R2T::selfOuterProduct(A.eigvecs().column(a)); \
192 : P += op_derivs[a] * Ma.outerProduct(Ma); \
193 : } \
194 : \
195 : for (auto a : make_range(3)) \
196 : for (auto b : make_range(a)) \
197 : { \
198 : const auto Ma = R2T::selfOuterProduct(A.eigvecs().column(a)); \
199 : const auto Mb = R2T::selfOuterProduct(A.eigvecs().column(b)); \
200 : \
201 : usingTensorIndices(i_, j_, k_, l_); \
202 : Gab = Ma.template times<i_, k_, j_, l_>(Mb) + Ma.template times<i_, l_, j_, k_>(Mb); \
203 : Gba = Mb.template times<i_, k_, j_, l_>(Ma) + Mb.template times<i_, l_, j_, k_>(Ma); \
204 : \
205 : Real theta_ab; \
206 : if (!MooseUtils::absoluteFuzzyEqual(A.eigvals()[a], A.eigvals()[b])) \
207 : theta_ab = 0.5 * (op_eigvals[a] - op_eigvals[b]) / (A.eigvals()[a] - A.eigvals()[b]); \
208 : else \
209 : theta_ab = 0.25 * (op_derivs[a] + op_derivs[b]); \
210 : \
211 : P += theta_ab * (Gab + Gba); \
212 : } \
213 : return P; \
214 : }
215 :
216 : #define FactorizedRankTwoTensorOperatorMapDerivativeUnary(derivativename, operator, derivative) \
217 : template <typename T> \
218 : RankFourTensorTempl<typename T::value_type> derivativename( \
219 : const FactorizedRankTwoTensorTempl<T> & A) \
220 : { \
221 : FactorizedRankTwoTensorOperatorMapDerivativeBody(operator, derivative) \
222 : } \
223 : static_assert(true, "")
224 :
225 : #define FactorizedRankTwoTensorOperatorMapDerivativeBinary(derivativename, operator, derivative) \
226 : template <typename T, typename T2> \
227 : RankFourTensorTempl<typename T::value_type> derivativename( \
228 : const FactorizedRankTwoTensorTempl<T> & A, const T2 & arg) \
229 : { \
230 : if constexpr (libMesh::ScalarTraits<T2>::value) \
231 : { \
232 : FactorizedRankTwoTensorOperatorMapDerivativeBody(operator, derivative) \
233 : } \
234 : } \
235 : static_assert(true, "")
236 :
237 : // TODO: While the macros are here, in the future we could instantiate other operator maps like
238 : // trignometry functions.
239 : // @{
240 10 : FactorizedRankTwoTensorOperatorMapUnary(log, log(eigval));
241 28 : FactorizedRankTwoTensorOperatorMapDerivativeUnary(dlog, log(eigval), 1 / eigval);
242 :
243 10 : FactorizedRankTwoTensorOperatorMapUnary(exp, exp(eigval));
244 28 : FactorizedRankTwoTensorOperatorMapDerivativeUnary(dexp, exp(eigval), exp(eigval));
245 :
246 10 : FactorizedRankTwoTensorOperatorMapUnary(sqrt, sqrt(eigval));
247 28 : FactorizedRankTwoTensorOperatorMapDerivativeUnary(dsqrt, sqrt(eigval), pow(eigval, -1. / 2.) / 2.);
248 :
249 10 : FactorizedRankTwoTensorOperatorMapUnary(cbrt, cbrt(eigval));
250 28 : FactorizedRankTwoTensorOperatorMapDerivativeUnary(dcbrt, cbrt(eigval), pow(eigval, -2. / 3.) / 3.);
251 :
252 10 : FactorizedRankTwoTensorOperatorMapBinary(pow, pow(eigval, arg));
253 28 : FactorizedRankTwoTensorOperatorMapDerivativeBinary(dpow,
254 : pow(eigval, arg),
255 : arg * pow(eigval, arg - 1));
256 : // @}
257 : } // end namespace MathUtils
258 :
259 : template <typename T>
260 : template <typename T2>
261 : FactorizedRankTwoTensorTempl<T>
262 2 : FactorizedRankTwoTensorTempl<T>::operator*(const T2 & a) const
263 : {
264 : if constexpr (libMesh::ScalarTraits<T2>::value)
265 : {
266 2 : FactorizedRankTwoTensorTempl<T> A = *this;
267 8 : for (auto & eigval : A._eigvals)
268 6 : eigval *= a;
269 4 : return A;
270 2 : }
271 : }
272 :
273 : template <typename T>
274 : template <typename T2>
275 : FactorizedRankTwoTensorTempl<T>
276 2 : FactorizedRankTwoTensorTempl<T>::operator/(const T2 & a) const
277 : {
278 : if constexpr (libMesh::ScalarTraits<T2>::value)
279 : {
280 2 : FactorizedRankTwoTensorTempl<T> A = *this;
281 8 : for (auto & eigval : A._eigvals)
282 6 : eigval /= a;
283 4 : return A;
284 2 : }
285 : }
286 :
287 : typedef FactorizedRankTwoTensorTempl<RankTwoTensor> FactorizedRankTwoTensor;
288 : typedef FactorizedRankTwoTensorTempl<ADRankTwoTensor> ADFactorizedRankTwoTensor;
289 : typedef FactorizedRankTwoTensorTempl<SymmetricRankTwoTensor> FactorizedSymmetricRankTwoTensor;
290 : typedef FactorizedRankTwoTensorTempl<ADSymmetricRankTwoTensor> ADFactorizedSymmetricRankTwoTensor;
|
