suanPan
tensor.h
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2017-2024 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
29#ifndef TENSOR_H
30#define TENSOR_H
31
32#include <Toolbox/utility.h>
33
34template<typename T> class Quaternion;
35
36namespace tensor {
37 mat isotropic_stiffness(double, double);
38 mat orthotropic_stiffness(const vec&, const vec&);
39
43
44 static const vec unit_tensor2{1., 1., 1., 0., 0., 0.};
45
46 namespace stress {
47 // applies to 3D tensor only, either principal or not
48 double invariant1(const vec&);
49 // applies to 3D tensor only, either principal or not
50 double invariant2(const vec&);
51 // applies to 3D tensor only, either principal or not
52 double invariant3(const vec&);
53
54 // compute lode angle based on input stress
55 double lode(vec);
56 // compute derivative of lode angle based on input stress
57 vec lode_der(vec);
58
59 static const vec norm_weight{1., 1., 1., 2., 2., 2.};
60 } // namespace stress
61 namespace strain {
62 // applies to 3D tensor only, either principal or not
63 double invariant1(const vec&);
64 // applies to 3D tensor only, either principal or not
65 double invariant2(const vec&);
66 // applies to 3D tensor only, either principal or not
67 double invariant3(const vec&);
68
69 // compute lode angle based on input deviatoric strain
70 double lode(vec);
71
72 static const vec norm_weight{1., 1., 1., .5, .5, .5};
73 } // namespace strain
74 double trace2(const vec&);
75 double trace3(const vec&);
76 double mean3(const vec&);
77 vec dev(const vec&);
78 vec dev(vec&&);
79
80 mat dev(const mat&);
81 mat dev(mat&&);
82
83 namespace strain {
84 mat to_green(mat&&);
85 mat to_green(const mat&);
86 mat to_tensor(const vec&);
87 vec to_voigt(const mat&);
88 double norm(const vec&);
89 double norm(vec&&);
90 double double_contraction(const vec&, const vec&);
91 double double_contraction(vec&&, vec&&);
92 } // namespace strain
93 namespace stress {
94 mat to_tensor(const vec&);
95 vec to_voigt(const mat&);
96 double norm(const vec&);
97 double norm(vec&&);
98 double double_contraction(const vec&, const vec&);
99 double double_contraction(vec&&, vec&&);
100 } // namespace stress
101
102 namespace base {
103 class Base3D {
104 const vec3 g1, g2, g3;
105 mat33 g;
106
107 public:
108 Base3D(const vec3&, const vec3&, const vec3&);
109 [[nodiscard]] std::tuple<vec3, vec3, vec3> to_inverse() const;
110 };
111
112 vec3 unit_norm(const vec3&, const vec3&);
113 } // namespace base
114
115 mat diff_unit(const vec&);
116
117 mat diff_triad(const vec3&, const vec3&, const vec3&);
118} // namespace tensor
119
120namespace transform {
121 void hoffman_projection(const vec&, mat&, mat&);
122 mat hill_projection(double, double, double, double, double, double);
123
124 double atan2(const vec&);
127
128 mat eigen_to_tensor_base(const mat&);
129 mat eigen_to_tensile_stress(const vec&, const mat&);
130 mat eigen_to_tensile_derivative(const vec&, const mat&);
131
132 template<typename T> Mat<T> skew_symm(const Mat<T>& R) {
133 suanpan_assert([&] { if(R.n_elem != 3) throw invalid_argument("need 3 element vector"); });
134
135 Mat<T> S(3, 3, fill::zeros);
136
137 S(0, 1) = -(S(1, 0) = R(2));
138 S(2, 0) = -(S(0, 2) = R(1));
139 S(1, 2) = -(S(2, 1) = R(0));
140
141 return S;
142 }
143
144 template<typename T> concept HasEval = requires(const T& x) { { x.eval() } -> std::convertible_to<mat>; };
145
146 template<HasEval T> mat skew_symm(const T& R) { return skew_symm(R.eval()); }
147
148 template<typename T> Mat<T> rodrigues(const Mat<T>& R) { return arma::expmat(transform::skew_symm(R)); }
149
150 template<typename T> Quaternion<T> to_quaternion(const Mat<T>& R) {
151 if(3 == R.n_elem) {
152 const auto angle = arma::norm(R);
153
154 if(suanpan::approx_equal(angle, 0.)) return {0., 0., 0., 0.};
155
156 return {std::cos(.5 * angle), std::sin(.5 * angle) / angle * R};
157 }
158
159 if(9 == R.n_elem) {
160 const auto tr_r = arma::trace(R);
161 const auto max_r = arma::max(R.diag());
162 if(tr_r >= max_r) {
163 const auto q0 = .5 * std::sqrt(tr_r + 1.);
164 const auto q1 = .25 / q0 * (R(2, 1) - R(1, 2));
165 const auto q2 = .25 / q0 * (R(0, 2) - R(2, 0));
166 const auto q3 = .25 / q0 * (R(1, 0) - R(0, 1));
167 return {q0, q1, q2, q3};
168 }
169
170 for(auto I = 0; I < 3; ++I)
171 if(suanpan::approx_equal(R(I, I), max_r)) {
172 const auto J = (I + 1) % 3;
173 const auto K = (J + 1) % 3;
174 vec q(3, fill::none);
175 q(I) = std::sqrt(.5 * max_r + .25 * (1. - tr_r));
176 const double q0 = .25 / q(I) * (R(K, J) - R(J, K));
177 q(J) = .25 / q(I) * (R(J, I) + R(I, J));
178 q(K) = .25 / q(I) * (R(K, I) + R(I, K));
179
180 return {q0, std::move(q)};
181 }
182 }
183
184 throw invalid_argument("need either rotation vector or matrix");
185 }
186
187 template<typename T> Col<T> to_pseudo(const Mat<T>& R) {
188 const Mat<T> S = arma::real(arma::logmat(R));
189
190 return {S(2, 1), S(0, 2), S(1, 0)};
191 }
192
193 namespace strain {
194 double angle(const vec&);
195 mat trans(double);
196 vec principal(const vec&);
197 vec rotate(const vec&, double);
198 } // namespace strain
199 namespace stress {
200 double angle(const vec&);
201 mat trans(double);
202 vec principal(const vec&);
203 vec rotate(const vec&, double);
204 } // namespace stress
205 namespace beam {
206 mat global_to_local(double, double, double);
207 mat global_to_local(const vec&, double);
208 } // namespace beam
209 namespace triangle {
210 vec to_area_coordinate(const vec&, const mat&);
211 }
212} // namespace transform
213
214namespace suanpan {
215 template<typename T> T ramp(const T in) { return in > T(0) ? in : T(0); }
216} // namespace suanpan
217
218#endif
219
An Quaternion class.
Definition: Quaternion.hpp:34
Definition: tensor.h:103
std::tuple< vec3, vec3, vec3 > to_inverse() const
Definition: tensor.cpp:283
Base3D(const vec3 &, const vec3 &, const vec3 &)
Definition: tensor.cpp:274
Definition: tensor.h:144
Mat< T > stress(T X, T Y, unsigned S)
Definition: shape.h:475
Mat< T > strain(T X, T Y, T V, unsigned S)
Definition: shape.h:523
Col< T > beam(T int_pts, unsigned order, double length)
Definition: shape.h:151
T triangle(const Mat< T > &EC)
Definition: shape.h:107
Definition: MatrixModifier.hpp:36
T ramp(const T in)
Definition: tensor.h:215
std::enable_if_t<!std::numeric_limits< T >::is_integer, bool > approx_equal(T x, T y, int ulp=2)
Definition: utility.h:60
vec3 unit_norm(const vec3 &, const vec3 &)
Definition: tensor.cpp:289
mat to_green(mat &&)
Definition: tensor.cpp:292
mat to_tensor(const vec &)
Definition: tensor.cpp:316
double invariant3(const vec &)
compute the third invariant of the given 3D strain tensor, could be either normal or deviatoric strai...
Definition: tensor.cpp:112
double norm(const vec &)
Definition: tensor.cpp:370
double invariant2(const vec &)
compute the second invariant of the given 3D strain tensor, could be either normal or deviatoric stra...
Definition: tensor.cpp:100
double invariant1(const vec &)
compute the first invariant of the given 3D strain tensor, could be either normal or deviatoric strai...
Definition: tensor.cpp:89
double double_contraction(const vec &, const vec &)
Definition: tensor.cpp:382
vec to_voigt(const mat &)
Definition: tensor.cpp:341
double lode(vec)
Definition: tensor.cpp:154
vec to_voigt(const mat &)
Definition: tensor.cpp:411
double invariant2(const vec &)
compute the second invariant of the given 3D stress tensor, could be either normal or deviatoric stre...
Definition: tensor.cpp:135
double double_contraction(const vec &, const vec &)
Definition: tensor.cpp:452
vec lode_der(vec)
Definition: tensor.cpp:172
double lode(vec)
Definition: tensor.cpp:163
double invariant3(const vec &)
compute the third invariant of the given 3D stress tensor, could be either normal or deviatoric stres...
Definition: tensor.cpp:147
double invariant1(const vec &)
compute the first invariant of the given 3D stress tensor, could be either normal or deviatoric stres...
Definition: tensor.cpp:124
double norm(const vec &)
Definition: tensor.cpp:440
mat to_tensor(const vec &)
Definition: tensor.cpp:386
double norm(vec &&)
Definition: tensor.cpp:446
Definition: tensor.h:36
double trace3(const vec &)
Only accepts 3D tensor!
Definition: tensor.cpp:196
mat diff_unit(const vec &)
Definition: tensor.cpp:231
mat unit_deviatoric_tensor4v2()
Definition: tensor.cpp:67
mat unit_deviatoric_tensor4()
Definition: tensor.cpp:57
double mean3(const vec &)
Definition: tensor.cpp:202
mat diff_triad(const vec3 &, const vec3 &, const vec3 &)
Definition: tensor.cpp:236
vec dev(const vec &)
Definition: tensor.cpp:204
mat unit_symmetric_tensor4()
Definition: tensor.cpp:75
mat isotropic_stiffness(double, double)
Definition: tensor.cpp:20
mat orthotropic_stiffness(const vec &, const vec &)
Definition: tensor.cpp:34
double trace2(const vec &)
Only accepts 2D tensor!
Definition: tensor.cpp:185
mat global_to_local(double, double, double)
Definition: tensor.cpp:708
vec principal(const vec &)
Definition: tensor.cpp:660
mat trans(double)
Definition: tensor.cpp:646
double angle(const vec &)
Definition: tensor.cpp:640
vec rotate(const vec &, double)
Definition: tensor.cpp:672
vec principal(const vec &)
Definition: tensor.cpp:694
mat trans(double)
Definition: tensor.cpp:680
double angle(const vec &)
Definition: tensor.cpp:678
vec rotate(const vec &, double)
Definition: tensor.cpp:706
vec to_area_coordinate(const vec &, const mat &)
Definition: tensor.cpp:625
Definition: tensor.h:120
mat compute_jacobian_nominal_to_principal(const mat &)
Definition: tensor.cpp:511
Mat< T > skew_symm(const Mat< T > &R)
Definition: tensor.h:132
mat hill_projection(double, double, double, double, double, double)
Definition: tensor.cpp:488
Col< T > to_pseudo(const Mat< T > &R)
Definition: tensor.h:187
double atan2(const vec &)
Definition: tensor.cpp:509
mat compute_jacobian_principal_to_nominal(const mat &)
Definition: tensor.cpp:545
mat eigen_to_tensile_stress(const vec &, const mat &)
Definition: tensor.cpp:596
void hoffman_projection(const vec &, mat &, mat &)
Generate two projection matrix based on the given yield stress according to the Hoffman yielding crit...
Definition: tensor.cpp:462
mat eigen_to_tensor_base(const mat &)
Definition: tensor.cpp:579
mat eigen_to_tensile_derivative(const vec &, const mat &)
Definition: tensor.cpp:603
Quaternion< T > to_quaternion(const Mat< T > &R)
Definition: tensor.h:150
Mat< T > rodrigues(const Mat< T > &R)
Definition: tensor.h:148
void suanpan_assert(const std::function< void()> &F)
Definition: suanPan.h:296