v0.15.0
Loading...
Searching...
No Matches
EshelbianAux.hpp
Go to the documentation of this file.
1/**
2 * @file EshelbianAux.hpp
3 * @brief Auxilary functions for Eshelbian plasticity
4 * @date 2024-08-30
5 *
6 * @copyright Copyright (c) 2024
7 *
8 */
9
10#include <Lie.hpp>
11
12namespace EshelbianPlasticity {
13
14inline auto diff_deviator(FTensor::Ddg<double, 3, 3> &&t_diff_stress) {
15 FTensor::Index<'i', 3> i;
16 FTensor::Index<'j', 3> j;
17 FTensor::Index<'k', 3> k;
18 FTensor::Index<'l', 3> l;
19
20 FTensor::Ddg<double, 3, 3> t_diff_deviator;
21 t_diff_deviator(i, j, k, l) = t_diff_stress(i, j, k, l);
22
23 constexpr double third = boost::math::constants::third<double>();
24
25 t_diff_deviator(0, 0, 0, 0) -= third;
26 t_diff_deviator(0, 0, 1, 1) -= third;
27
28 t_diff_deviator(1, 1, 0, 0) -= third;
29 t_diff_deviator(1, 1, 1, 1) -= third;
30
31 t_diff_deviator(2, 2, 0, 0) -= third;
32 t_diff_deviator(2, 2, 1, 1) -= third;
33
34 t_diff_deviator(0, 0, 2, 2) -= third;
35 t_diff_deviator(1, 1, 2, 2) -= third;
36 t_diff_deviator(2, 2, 2, 2) -= third;
37
38 return t_diff_deviator;
39}
40
41constexpr static auto size_symm = (3 * (3 + 1)) / 2;
42
43inline auto diff_tensor() {
44 FTensor::Index<'i', 3> i;
45 FTensor::Index<'j', 3> j;
46 FTensor::Index<'k', 3> k;
47 FTensor::Index<'l', 3> l;
50 t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
51 return t_diff;
52};
53
54inline auto symm_L_tensor() {
55 FTensor::Index<'i', 3> i;
56 FTensor::Index<'j', 3> j;
59 t_L(i, j, L) = 0;
60 t_L(0, 0, 0) = 1;
61 t_L(1, 1, 3) = 1;
62 t_L(2, 2, 5) = 1;
63 t_L(1, 0, 1) = 1;
64 t_L(2, 0, 2) = 1;
65 t_L(2, 1, 4) = 1;
66 return t_L;
67}
68
69inline auto voigt_to_symm() {
73 t_L(L, M) = 0;
74 t_L(0, 0) = 1;
75 t_L(1, 3) = 1;
76 t_L(2, 5) = 1;
77 t_L(3, 1) = 1;
78 t_L(4, 2) = 1;
79 t_L(5, 4) = 1;
80 return t_L;
81}
82
83inline auto diff_symmetrize() {
84
89
91
92 t_diff(i, j, k, l) = 0;
93 t_diff(0, 0, 0, 0) = 1;
94 t_diff(1, 1, 1, 1) = 1;
95
96 t_diff(1, 0, 1, 0) = 0.5;
97 t_diff(1, 0, 0, 1) = 0.5;
98
99 t_diff(0, 1, 0, 1) = 0.5;
100 t_diff(0, 1, 1, 0) = 0.5;
101
102 if constexpr (SPACE_DIM == 3) {
103 t_diff(2, 2, 2, 2) = 1;
104
105 t_diff(2, 0, 2, 0) = 0.5;
106 t_diff(2, 0, 0, 2) = 0.5;
107 t_diff(0, 2, 0, 2) = 0.5;
108 t_diff(0, 2, 2, 0) = 0.5;
109
110 t_diff(2, 1, 2, 1) = 0.5;
111 t_diff(2, 1, 1, 2) = 0.5;
112 t_diff(1, 2, 1, 2) = 0.5;
113 t_diff(1, 2, 2, 1) = 0.5;
114 }
115
116 return t_diff;
117}
118
119template <class T> struct TensorTypeExtractor {
120 typedef typename std::remove_pointer<T>::type Type;
121};
122template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
123 typedef typename std::remove_pointer<T>::type Type;
124};
125
126struct isEq {
127 static inline auto check(const double &a, const double &b) {
128 constexpr double eps = std::numeric_limits<float>::epsilon();
129 return std::abs(a - b) / absMax < eps;
130 }
131 static inline double absMax = 1.;
132};
133
134inline auto is_eq(const double &a, const double &b) {
135 return isEq::check(a, b);
136};
137
138template <int DIM> inline auto get_uniq_nb(double *ptr) {
139 std::array<double, DIM> tmp;
140 std::copy(ptr, &ptr[DIM], tmp.begin());
141 std::sort(tmp.begin(), tmp.end());
143 std::max(std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1])), 1.);
144 return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
145}
146
147template <typename A, typename B>
148inline auto sort_eigen_vals(A &eig, B &eigen_vec) {
149
150 constexpr int dim = 3;
151
153 std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
154 constexpr double eps = std::numeric_limits<double>::epsilon();
155 isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
156
157 int i = 0, j = 1, k = 2;
158
159 if (is_eq(eig(0), eig(1))) {
160 i = 0;
161 j = 2;
162 k = 1;
163 } else if (is_eq(eig(0), eig(2))) {
164 i = 0;
165 j = 1;
166 k = 2;
167 } else if (is_eq(eig(1), eig(2))) {
168 i = 1;
169 j = 0;
170 k = 2;
171 }
172
174 eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
175
176 eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
177
178 eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
179
180 FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
181
182 {
183 FTENSOR_INDEX(dim, i);
184 FTENSOR_INDEX(dim, j);
185 eig(i) = eig_c(i);
186 eigen_vec(i, j) = eigen_vec_c(i, j);
187 }
188}
189}
constexpr double third
Lie algebra implementation.
#define FTENSOR_INDEX(DIM, I)
constexpr double a
static const double eps
constexpr int SPACE_DIM
Kronecker Delta class symmetric.
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto get_uniq_nb(double *ptr)
auto diff_deviator(FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
auto sort_eigen_vals(A &eig, B &eigen_vec)
static constexpr auto size_symm
auto is_eq(const double &a, const double &b)
Tensors class implemented by Walter Landry.
Definition FTensor.hpp:51
std::remove_pointer< T >::type Type
static auto check(const double &a, const double &b)