v0.10.0
ElasticOperators.hpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 #define MAT_TO_TENSOR4_3D(MAT) \
16  &MAT(0, 0), &MAT(1, 0), &MAT(2, 0), &MAT(3, 0), &MAT(4, 0), &MAT(5, 0), \
17  &MAT(6, 0), &MAT(7, 0), &MAT(8, 0), &MAT(9, 0), &MAT(10, 0), \
18  &MAT(11, 0), &MAT(12, 0), &MAT(13, 0), &MAT(14, 0), &MAT(15, 0), \
19  &MAT(16, 0), &MAT(17, 0), &MAT(18, 0), &MAT(19, 0), &MAT(20, 0), \
20  &MAT(21, 0), &MAT(22, 0), &MAT(23, 0), &MAT(24, 0), &MAT(25, 0), \
21  &MAT(26, 0), &MAT(27, 0), &MAT(28, 0), &MAT(29, 0), &MAT(30, 0), \
22  &MAT(31, 0), &MAT(32, 0), &MAT(33, 0), &MAT(34, 0), &MAT(35, 0), \
23  &MAT(36, 0), &MAT(37, 0), &MAT(38, 0), &MAT(39, 0), &MAT(40, 0), \
24  &MAT(41, 0), &MAT(42, 0), &MAT(43, 0), &MAT(44, 0), &MAT(45, 0), \
25  &MAT(46, 0), &MAT(47, 0), &MAT(48, 0), &MAT(49, 0), &MAT(50, 0), \
26  &MAT(51, 0), &MAT(52, 0), &MAT(53, 0), &MAT(54, 0), &MAT(55, 0), \
27  &MAT(56, 0), &MAT(57, 0), &MAT(58, 0), &MAT(59, 0), &MAT(60, 0), \
28  &MAT(61, 0), &MAT(62, 0), &MAT(63, 0), &MAT(64, 0), &MAT(65, 0), \
29  &MAT(66, 0), &MAT(67, 0), &MAT(68, 0), &MAT(69, 0), &MAT(70, 0), \
30  &MAT(71, 0), &MAT(72, 0), &MAT(73, 0), &MAT(74, 0), &MAT(75, 0), \
31  &MAT(76, 0), &MAT(77, 0), &MAT(78, 0), &MAT(79, 0), &MAT(80, 0)
32 
33 template <typename T> inline auto to_non_symm(T &symm) {
34  Tensor2<double, 3, 3> non_symm;
35  Number<0> N0;
36  Number<1> N1;
37  Number<2> N2;
38  non_symm(N0, N0) = symm(N0, N0);
39  non_symm(N1, N1) = symm(N1, N1);
40  non_symm(N2, N2) = symm(N2, N2);
41  non_symm(N0, N1) = non_symm(N1, N0) = symm(N0, N1);
42  non_symm(N0, N2) = non_symm(N2, N0) = symm(N0, N2);
43  non_symm(N2, N1) = non_symm(N1, N2) = symm(N2, N1);
44  return non_symm;
45 };
46 
47 template <typename T> inline auto to_symm(T &non_symm) {
48  Tensor2_symmetric<double, 3> symm;
49  Number<0> N0;
50  Number<1> N1;
51  Number<2> N2;
52  symm(N0, N0) = non_symm(N0, N0);
53  symm(N1, N1) = non_symm(N1, N1);
54  symm(N2, N2) = non_symm(N2, N2);
55  symm(N0, N1) = non_symm(N0, N1);
56  symm(N0, N2) = non_symm(N0, N2);
57  symm(N1, N2) = non_symm(N1, N2);
58  return symm;
59 };
60 
61 // template <size_t DIM>
62 // inline auto getPointerPackVector(array<double, MAX_DOFS_ON_ENTITY> &nf);
63 // template <>
64 // inline auto getPointerPackVector<3>(array<double, MAX_DOFS_ON_ENTITY> &nf) {
65 // return Tensor1<PackPtr<double *, 3>, 3>{&nf[0], &nf[1], &nf[2]};
66 // };
67 // template <>
68 // inline auto getPointerPackVector<3>(array<double, MAX_DOFS_ON_ENTITY> &nf) {
69 // return Tensor1<PackPtr<double *, 2>, 2>{&nf[0], &nf[1]};
70 // };
71 
72 inline bool is_eq(const double &a, const double &b) {
73  constexpr double eps = std::numeric_limits<double>::epsilon();
74  return abs(a - b) < eps;
75 };
76 inline bool is_diff(const double &a, const double &b) {
77  constexpr double eps = std::numeric_limits<double>::epsilon();
78  return abs(a - b) > eps;
79 };
80 
81 template <typename T>
82 inline Ddg<double, 3, 3> tensor4_to_ddg(Tensor4<T, 3, 3, 3, 3> &t) {
83  Ddg<double, 3, 3> loc_tens;
84 
85  Number<0> N0;
86  Number<1> N1;
87  Number<2> N2;
88 
89  loc_tens(N0, N0, N0, N0) = t(N0, N0, N0, N0); // 0
90  loc_tens(N0, N0, N0, N1) = t(N0, N0, N0, N1); // 1
91  loc_tens(N0, N0, N0, N2) = t(N0, N0, N0, N2); // 2
92  loc_tens(N0, N0, N1, N1) = t(N0, N0, N1, N1); // 3
93  loc_tens(N0, N0, N1, N2) = t(N0, N0, N1, N2); // 4
94  loc_tens(N0, N0, N2, N2) = t(N0, N0, N2, N2); // 5
95  loc_tens(N0, N1, N0, N0) = t(N0, N1, N0, N0); // 6
96  loc_tens(N0, N1, N0, N1) = t(N0, N1, N0, N1); // 7
97  loc_tens(N0, N1, N0, N2) = t(N0, N1, N0, N2); // 8
98  loc_tens(N0, N1, N1, N1) = t(N0, N1, N1, N1); // 9
99  loc_tens(N0, N1, N1, N2) = t(N0, N1, N1, N2); // 10
100  loc_tens(N0, N1, N2, N2) = t(N0, N1, N2, N2); // 11
101  loc_tens(N0, N2, N0, N0) = t(N0, N2, N0, N0); // 12
102  loc_tens(N0, N2, N0, N1) = t(N0, N2, N0, N1); // 13
103  loc_tens(N0, N2, N0, N2) = t(N0, N2, N0, N2); // 14
104  loc_tens(N0, N2, N1, N1) = t(N0, N2, N1, N1); // 15
105  loc_tens(N0, N2, N1, N2) = t(N0, N2, N1, N2); // 16
106  loc_tens(N0, N2, N2, N2) = t(N0, N2, N2, N2); // 17
107  loc_tens(N1, N1, N0, N0) = t(N1, N1, N0, N0); // 18
108  loc_tens(N1, N1, N0, N1) = t(N1, N1, N0, N1); // 19
109  loc_tens(N1, N1, N0, N2) = t(N1, N1, N0, N2); // 20
110  loc_tens(N1, N1, N1, N1) = t(N1, N1, N1, N1); // 21
111  loc_tens(N1, N1, N1, N2) = t(N1, N1, N1, N2); // 22
112  loc_tens(N1, N1, N2, N2) = t(N1, N1, N2, N2); // 23
113  loc_tens(N1, N2, N0, N0) = t(N1, N2, N0, N0); // 24
114  loc_tens(N1, N2, N0, N1) = t(N1, N2, N0, N1); // 25
115  loc_tens(N1, N2, N0, N2) = t(N1, N2, N0, N2); // 26
116  loc_tens(N1, N2, N1, N1) = t(N1, N2, N1, N1); // 27
117  loc_tens(N1, N2, N1, N2) = t(N1, N2, N1, N2); // 28
118  loc_tens(N1, N2, N2, N2) = t(N1, N2, N2, N2); // 29
119  loc_tens(N2, N2, N0, N0) = t(N2, N2, N0, N0); // 30
120  loc_tens(N2, N2, N0, N1) = t(N2, N2, N0, N1); // 31
121  loc_tens(N2, N2, N0, N2) = t(N2, N2, N0, N2); // 32
122  loc_tens(N2, N2, N1, N1) = t(N2, N2, N1, N1); // 33
123  loc_tens(N2, N2, N1, N2) = t(N2, N2, N1, N2); // 34
124  loc_tens(N2, N2, N2, N2) = t(N2, N2, N2, N2); // 35
125 
126  return loc_tens;
127 };
128 
129 template <typename T> inline auto get_tensor4_from_ddg(Ddg<T, 3, 3> &ddg) {
130  Tensor4<double, 3, 3, 3, 3> tens4;
131  for (int ii = 0; ii != 3; ++ii)
132  for (int jj = 0; jj != 3; ++jj)
133  for (int kk = 0; kk != 3; ++kk)
134  for (int ll : {0, 1, 2})
135  tens4(ii, jj, kk, ll) = ddg(ii, jj, kk, ll);
136  return tens4;
137 };
138 
139 template <typename T> inline bool is_ddg(Tensor4<T, 3, 3, 3, 3> tensor) {
140  auto test_ddg = tensor4_to_ddg(tensor);
141  for (int ii = 0; ii != 3; ++ii)
142  for (int jj = 0; jj != 3; ++jj)
143  for (int kk = 0; kk != 3; ++kk)
144  for (int ll : {0, 1, 2}) {
145  const double diff =
146  abs(test_ddg(ii, jj, kk, ll) - tensor(ii, jj, kk, ll));
147  if (diff > 1e-6)
148  return false;
149  }
150 
151  return true;
152 };
153 
154 template <typename T> inline auto exp_map(const Tensor2_symmetric<T, 3> &X) {
155  Tensor2<double, 3, 3> exp_tensor(1., 0., 0., 0., 1., 0., 0., 0., 1.);
156  Tensor2<double, 3, 3> X_nt;
157 
158  constexpr int max_it = 1000;
159  constexpr double eps = 1e-12;
160 
161  int fac_n = 1;
162 
163  auto Xcopy = to_non_symm(X);
164  auto X_n = to_non_symm(X);
165  exp_tensor(i, j) += Xcopy(i, j);
166 
167  for (int n = 2; n != max_it; ++n) {
168  fac_n *= n;
169  X_nt(i, j) = X_n(i, j);
170  X_n(i, k) = X_nt(i, j) * Xcopy(j, k);
171  const double inv_nf = 1. / fac_n;
172  exp_tensor(i, j) += inv_nf * X_n(i, j);
173  const double e = inv_nf * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
174  if (e < eps)
175  return to_symm(exp_tensor);
176  }
177 
178  return to_symm(exp_tensor);
179 };
180 
181 inline auto get_rotation_R(Tensor1<double, 3> t1_omega, double tt) {
182 
184 
185  auto get_rotation = [&](FTensor1 &t_omega) {
187 
188  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
189  t_R(i, j) = t_kd(i, j);
190 
191  const double angle = sqrt(t_omega(i) * t_omega(i));
192  if (std::abs(angle) < 1e-18)
193  return t_R;
194 
196  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
197  const double a = sin(angle) / angle;
198  const double ss_2 = sin(angle / 2.);
199  const double b = 2. * ss_2 * ss_2 / (angle * angle);
200  t_R(i, j) += a * t_Omega(i, j);
201  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
202 
203  return t_R;
204  };
205 
206  Tensor1<double, 3> my_c_omega;
207  my_c_omega(i) = t1_omega(i) * tt;
208  auto t_R = get_rotation(my_c_omega);
209 
210  return get_rotation(my_c_omega);
211 };
212 
213 template <typename T>
214 inline auto get_log_map(const Tensor2_symmetric<T, 3> &X) {
215 
216  constexpr int max_it = 1000;
217  constexpr double eps = 1e-12;
218 
219  auto log_tensor = to_non_symm(X);
220  log_tensor(i, j) -= kronecker_delta(i, j);
221 
222  auto Xcopy = to_non_symm(log_tensor);
223  auto X_n = to_non_symm(log_tensor);
224  Tensor2<double, 3, 3> X_nt_tmp;
225 
226  for (int n = 2; n != max_it; ++n) {
227 
228  X_nt_tmp(i, j) = X_n(i, j);
229  X_n(i, k) = X_nt_tmp(i, j) * Xcopy(j, k);
230  const double inv_n = 1. / n;
231  log_tensor(i, j) += inv_n * pow(-1, n + 1) * X_n(i, j);
232  const double e = inv_n * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
233  if (e < eps)
234  break;
235  }
236 
237  return to_symm(log_tensor);
238 };
239 
240 template <typename T> inline auto get_diff_log_map(Tensor2_symmetric<T, 3> &X) {
241  Tensor2<double, 3, 3> X_nt_tmp;
242  Tensor4<double, 3, 3, 3, 3> diff_log_map;
243  Tensor4<double, 3, 3, 3, 3> diff_log_map_tmp;
244 
245  constexpr int max_it = 1000;
246  constexpr double eps = 1e-12;
247 
248  auto get_Xn = [&](auto &log_tensor) {
249  vector<Tensor2<double, 3, 3>> Xn;
250  Xn.emplace_back();
251  auto &cu_Xn = Xn.back();
252  cu_Xn(i, j) = kronecker_delta(i, j);
253  for (int n = 2; n != max_it; ++n) {
254  const double cof = 1. / n;
255  auto X_nt_tmp = Xn.back();
256  Xn.emplace_back();
257  auto &c_Xn = Xn.back();
258  c_Xn(i, k) = X_nt_tmp(i, j) * log_tensor(j, k);
259  const double e = cof * std::abs(sqrt(c_Xn(i, j) * c_Xn(i, j)));
260  if (e < eps)
261  return Xn;
262  }
263  return Xn;
264  };
265 
266  auto log_tensor = to_non_symm(X);
267  log_tensor(i, j) -= kronecker_delta(i, j);
268  diff_log_map(i, j, k, l) = 0.;
269 
270  auto Xn_vec = get_Xn(log_tensor);
271 
272  for (int n = 1; n != Xn_vec.size() + 1; ++n) {
273  const double cof = pow(-1, n + 1) / n;
274  for (int m = 1; m != n + 1; ++m) {
275  auto &Xn_1 = Xn_vec[m - 1];
276  auto &Xn_2 = Xn_vec[n - m];
277  diff_log_map(i, j, k, l) += cof * Xn_1(i, k) * Xn_2(l, j);
278  }
279  }
280  // FIXME: make it more effient
281  auto isddg = tensor4_to_ddg((*cache).Is);
282  diff_log_map_tmp(i, j, k, l) = isddg(i, j, m, n) * diff_log_map(m, n, k, l);
283  return tensor4_to_ddg(diff_log_map_tmp);
284 };
285 
286 template <typename T>
287 inline MoFEMErrorCode
288 get_eigen_val_and_proj_lapack(const Tensor2_symmetric<T, 3> &X,
289  Tensor1<T, 3> &lambda,
290  Tensor2<T, 3, 3> &eig_vec) {
291 
293  Tensor1<double, 3> eig;
294  Tensor2<double, 3, 3> eigen_vec_lap;
295  // eig_vec(i, j) = X(i, j);
296  eigen_vec_lap(0, 0) = X(0, 0);
297  eigen_vec_lap(1, 1) = X(1, 1);
298  eigen_vec_lap(2, 2) = X(2, 2);
299  eigen_vec_lap(0, 1) = eigen_vec_lap(1, 0) = X(0, 1);
300  eigen_vec_lap(2, 1) = eigen_vec_lap(1, 2) = X(2, 1);
301  eigen_vec_lap(2, 0) = eigen_vec_lap(0, 2) = X(2, 0);
302 
303  constexpr int dim = 3;
304  int n = dim, lda = dim, info, lwork = -1;
305 
306  double wkopt;
307  info = lapack_dsyev('V', 'U', n, &(eigen_vec_lap(0, 0)), lda, &eig(0), &wkopt,
308  lwork);
309  if (info != 0)
310  SETERRQ1(PETSC_COMM_SELF, 1,
311  "something is wrong with lapack_dsyev info = %d", info);
312 
313  lwork = (int)wkopt;
314  double work[lwork];
315  info = lapack_dsyev('V', 'U', n, &(eigen_vec_lap(0, 0)), lda, &eig(0), work,
316  lwork);
317  if (info != 0)
318  SETERRQ1(PETSC_COMM_SELF, 1,
319  "something is wrong with lapack_dsyev info = %d", info);
320  lambda(i) = eig(i);
321  eig_vec(i, j) = eigen_vec_lap(i, j);
322  // FIXME: perturb
323  // const double h = 1e-12;
324  // if (is_eq(lambda[0], lambda[1]))
325  // lambda[0] *= (1 + h);
326  // if (is_eq(lambda[1], lambda[2]))
327  // lambda[2] *= (1 - h);
328 
330 }
331 
332 template <typename T1, typename T2, typename T3>
333 inline auto get_ln_X_lapack(const Tensor2_symmetric<T1, 3> &X,
334  Tensor1<T2, 3> &lambda,
335  Tensor2<T3, 3, 3> &eig_vec) {
336  Tensor2_symmetric<double, 3> lnX;
337  Tensor2_symmetric<double, 3> diag(0., 0., 0., 0., 0., 0.);
338  for (int ii = 0; ii != 3; ++ii)
339  diag(ii, ii) = log(lambda(ii));
340 
341  lnX(i, j) = eig_vec(l, i) ^ (diag(l, k) * eig_vec(k, j));
342 
343  return lnX;
344 };
345 
346 template <typename T>
347 inline MoFEMErrorCode
348 get_eigen_val_and_proj_souza(const Tensor2_symmetric<T, 3> &X, double *e,
349  Tensor2<double, 3, 3> *E) {
350 
352  constexpr auto t_kd = Kronecker_Delta<int>();
353  constexpr double third = boost::math::constants::third<double>();
354  Tensor2<double, 3, 3> X2;
355  X2(i, j) = X(i, k) * X(k, j);
356 
357  const double I1 = X(i, i);
358  const double I2 = 0.5 * (I1 * I1 - X(j, i) * X(i, j));
359  double I3;
361 
362  const double R = (-2 * I1 * I1 * I1 + 9 * I1 * I2 - 27 * I3) / 54.;
363  const double Q = (I1 * I1 - 3 * I2) / 9.;
364  const double theta = acos(R / sqrt(Q * Q * Q));
365  double sqrtQ = sqrt(Q);
366  if (isnormal(sqrtQ)) {
367  e[0] = -2 * sqrtQ * cos(theta * third) + I1 * third;
368  e[1] = -2 * sqrtQ * cos((theta + 2 * M_PI) * third) + I1 * third;
369  e[2] = -2 * sqrtQ * cos((theta - 2 * M_PI) * third) + I1 * third;
370  } else {
371  e[0] = I1 * third;
372  e[1] = I1 * third;
373  e[2] = I1 * third;
374  }
375 
376  if (is_diff(e[0], e[1]) && is_diff(e[1], e[2]) && is_diff(e[0], e[2])) {
377  for (int ii = 0; ii != 3; ++ii) {
378  const double &x = e[ii];
379  const double coef = x / (2 * x * x * x - I1 * x * x + I3);
380  E[ii](i, j) = coef * (X2(i, j) - (I1 - x) * X(i, j) +
381  (I3 / x) * kronecker_delta(i, j));
382  }
383 
384  } else if (is_eq(e[0], e[1]) && is_eq(e[1], e[2])) {
385  E[0](i, j) = t_kd(i, j);
386  E[1](i, j) = 0.;
387  E[2](i, j) = 0.;
388  } else if (is_eq(e[1], e[2])) {
389  E[0](i, j) =
390  (e[0] / (2. * e[0] * e[0] * e[0] - I1 * e[0] * e[0] + I3)) *
391  (X(i, k) * X(k, j) - (I1 - e[0]) * X(i, j) + (I3 / e[0]) * t_kd(i, j));
392  E[1](i, j) = t_kd(i, j) - E[0](i, j);
393  E[2](i, j) = 0.;
394  } else if (is_eq(e[0], e[1])) {
395  E[2](i, j) =
396  (e[2] / (2. * e[2] * e[2] * e[2] - I1 * e[2] * e[2] + I3)) *
397  (X(i, k) * X(k, j) - (I1 - e[2]) * X(i, j) + (I3 / e[2]) * t_kd(i, j));
398  E[0](i, j) = t_kd(i, j) - E[2](i, j);
399  E[1](i, j) = 0.;
400  } else {
401  E[1](i, j) =
402  (e[1] / (2. * e[1] * e[1] * e[1] - I1 * e[1] * e[1] + I3)) *
403  (X(i, k) * X(k, j) - (I1 - e[1]) * X(i, j) + (I3 / e[1]) * t_kd(i, j));
404  E[2](i, j) = t_kd(i, j) - E[1](i, j);
405  E[0](i, j) = 0.;
406  }
407 
409 }
410 
411 template <typename T> inline auto get_ln_X(const Tensor2_symmetric<T, 3> &X) {
412 
413  // Tensor2<double, 3, 3> lnX(0., 0., 0., 0., 0., 0., 0., 0., 0.);
414  // Tensor2<double, 3, 3> E[3];
415  // double lambda[3];
416  // auto nX = to_non_symm(X);
417  // CHKERR calcEigenvalues3x3(nX, lambda);
418  // CHKERR calcEigenproj3x3(nX, lambda, E);
419  // // CHKERR get_eigen_val_and_proj_souza(X, lambda, E);
420 
421  // for (int ii = 0; ii != 3; ii++)
422  // lnX(i, j) += E[ii](i, j) * log(lambda[ii]);
423 
424  auto test = get_log_map(X);
425  return get_log_map(X);
426  // return to_symm(lnX);
427 };
428 
429 template <typename T>
430 inline auto calculate_d_05lnX_dX(Tensor2_symmetric<T, 3> &X) {
431 
432  auto my_ln_x = get_diff_log_map(X);
433  my_ln_x(i, j, k, l) *= 0.5;
434  return my_ln_x;
435 
436  Tensor2<double, 3, 3> E[3];
437  double eig[3];
438  auto nX = to_non_symm(X);
439  CHKERR calcEigenvalues3x3(nX, eig);
440  CHKERR calcEigenproj3x3(nX, eig, E);
441  Tensor4<double, 3, 3, 3, 3> D_Xa;
442  CHKERR calcDerivativeLog3x3(nX, eig, E, D_Xa);
443  D_Xa(i, j, k, l) *= 0.5;
444  return tensor4_to_ddg(D_Xa);
445 };
446 
447 template <typename T, typename TP, typename T3>
448 inline auto
449 calculate_nominal_moduli_TL(Tensor2_symmetric<T, 3> &X, Tensor2<TP, 3, 3> &F,
450  Tensor2_symmetric<TP, 3> &stress,
451  Tensor1<T, 3> &lambda, Tensor2<T, 3, 3> &eig_vec,
452  Tensor4<T3, 3, 3, 3, 3> &TLs3,
453  bool compute_tangent) {
454 
455  // Tensor4<double, 3, 3, 3, 3> TLs3;
456  Tensor4<double, 3, 3, 3, 3> P3;
457  Tensor2<double, 3, 3> invF;
458 
459  double det_f;
460  CHKERR determinantTensor3by3(F, det_f);
461  CHKERR invertTensor3by3(F, det_f, invF);
462 
463  // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
464 
465  // VectorDouble lambda(3);
466  VectorDouble e(3);
467  VectorDouble f(3);
468  VectorDouble d(3);
469  MatrixDouble Vi(3, 3);
470  MatrixDouble Ksi(3, 3);
471  MatrixDouble Zeta(3, 3);
472  Vi.clear();
473  Ksi.clear();
474  Zeta.clear();
475  double eta = 0;
476 
477  Tensor1<double, 3> _N[3];
478  Tensor1<double, 3> _n[3];
479 
480  for (int dd : {0, 1, 2}) {
481  for (int cc : {0, 1, 2})
482  _N[dd](cc) = eig_vec(dd, cc);
483  }
484 
485  for (int dd : {0, 1, 2})
486  _n[dd](i) = F(i, j) * _N[dd](j);
487 
488  for (int ii = 0; ii != 3; ++ii) {
489  e(ii) = 0.5 * log(lambda(ii));
490  d(ii) = 1. / lambda(ii);
491  f(ii) = -2. / (lambda(ii) * lambda(ii));
492  }
493 
494  const double &lam0 = lambda(0);
495  const double &lam1 = lambda(1);
496  const double &lam2 = lambda(2);
497 
498  if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
499  // all different
500  for (int ii = 0; ii != 3; ++ii)
501  for (int jj = 0; jj != 3; ++jj) {
502  if (ii != jj) {
503 
504  Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
505  Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
506  for (int kk = 0; kk != 3; ++kk) {
507 
508  if (kk != ii && kk != jj)
509  eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
510  (lambda(ii) - lambda(kk)));
511  }
512  }
513  }
514 
515  } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
516  // all the same
517  for (int ii = 0; ii != 3; ++ii)
518  for (int jj = 0; jj != 3; ++jj)
519  if (ii != jj) {
520 
521  Vi(ii, jj) = 0.5 * d(ii);
522  Ksi(ii, jj) = 0.125 * f(ii);
523  }
524  eta = 0.125 * f(0);
525  } else if (is_eq(lam0, lam1)) {
526  int i_ = 0;
527  int j_ = 1;
528  int k_ = 2;
529  Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
530  Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
531 
532  for (int ii = 0; ii != 3; ++ii)
533  for (int jj = 0; jj != 3; ++jj)
534  if (ii != jj) {
535  if (jj == k_ || ii == k_) {
536  Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
537  Ksi(ii, jj) =
538  (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
539  }
540  }
541  eta = Ksi(k_, i_);
542  } else if (is_eq(lam0, lam2)) {
543  int i_ = 0;
544  int j_ = 2;
545  int k_ = 1;
546  Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
547  Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
548 
549  for (int ii = 0; ii != 3; ++ii)
550  for (int jj = 0; jj != 3; ++jj)
551  if (ii != jj) {
552  if (jj == k_ || ii == k_) {
553  Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
554  Ksi(ii, jj) =
555  (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
556  }
557  }
558  eta = Ksi(k_, i_);
559  } else if (is_eq(lam1, lam2)) {
560  int i_ = 1;
561  int j_ = 2;
562  int k_ = 0;
563  Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
564  Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
565 
566  for (int ii = 0; ii != 3; ++ii)
567  for (int jj = 0; jj != 3; ++jj)
568  if (ii != jj) {
569  if (jj == k_ || ii == k_) {
570  Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
571  Ksi(ii, jj) =
572  (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
573  }
574  }
575  eta = Ksi(k_, i_);
576  }
577 
578  P3(i, j, k, l) = 0.;
579 
580  Tensor2<double, 3, 3> Mii;
581  Tensor2<double, 3, 3> Mij;
582  Tensor2<double, 3, 3> Mik;
583  Tensor2<double, 3, 3> Mjk;
584  Tensor2<double, 3, 3> Mjj;
585 
586  // FIXME: TODO: put it nicely in one LOOP
587  for (int ii = 0; ii != 3; ++ii) {
588  Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
589  P3(i, j, k, l) += 0.5 * d(ii) * _N[ii](i) * _N[ii](j) * Mii(k, l);
590  for (int jj = 0; jj != 3; ++jj)
591  if (ii != jj) {
592  Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
593  P3(i, j, k, l) += Vi(ii, jj) * _N[ii](i) * _N[jj](j) * Mij(k, l);
594  }
595  }
596 
597  if (!compute_tangent)
598  return P3;
599 
600  TLs3(i, j, k, l) = 0.;
601  for (int ii = 0; ii != 3; ++ii)
602  for (int jj = 0; jj != 3; ++jj)
603  Zeta(ii, jj) = (stress(i, j) * (_N[ii](i) * _N[jj](j)));
604 
605  for (int ii = 0; ii != 3; ++ii) {
606  Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
607  TLs3(i, j, k, l) += 0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j) * Mii(k, l);
608  }
609 
610  for (int ii = 0; ii != 3; ++ii)
611  for (int jj = 0; jj != 3; ++jj)
612  if (jj != ii)
613  for (int kk = 0; kk != 3; ++kk)
614  if (kk != ii && kk != jj) {
615  Mik(i, j) = _n[ii](i) * _N[kk](j) + _n[kk](i) * _N[ii](j);
616  Mjk(i, j) = _n[jj](i) * _N[kk](j) + _n[kk](i) * _N[jj](j);
617  TLs3(i, j, k, l) += 2. * eta * Zeta(ii, jj) * Mik(i, j) * Mjk(k, l);
618  }
619 
620  for (int ii = 0; ii != 3; ++ii)
621  for (int jj = 0; jj != 3; ++jj)
622  if (jj != ii) {
623  Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
624  Mjj(i, j) = _n[jj](i) * _N[jj](j) + _n[jj](i) * _N[jj](j);
625  const double k2 = 2. * Ksi(ii, jj);
626  const double zp = k2 * Zeta(ii, jj);
627  const double zp2 = k2 * Zeta(jj, jj);
628  TLs3(i, j, k, l) += zp * Mij(i, j) * Mjj(k, l) +
629  zp * Mjj(i, j) * Mij(k, l) +
630  zp2 * Mij(i, j) * Mij(k, l);
631  }
632 
633  Tensor2<double, 3, 3> Piola;
634  Tensor2<double, 3, 3> invFPiola;
635  Tensor4<double, 3, 3, 3, 3> Ttmp3;
636  auto nstress = to_non_symm(stress);
637  Piola(k, l) = nstress(m, n) * P3(m, n, k, l);
638  invFPiola(i, j) = invF(i, k) * Piola(k, j);
639  Ttmp3(i, j, k, l) = 0.5 * invFPiola(i, k) * kronecker_delta(j, l) +
640  0.5 * invFPiola(i, l) * kronecker_delta(j, k);
641  TLs3(i, j, k, l) += Ttmp3(i, j, k, l);
642 
643  // TLs3(i, j, k, l) += invF(i, j) * (stress(m, n) * P3(m, n, k, l));
644 
645  // return TLs3;
646  return P3;
647 };
648 
649 namespace OpElasticTools {
650 
651 //! [Common data]
652 struct CommonData {
653  Ddg<double, 3, 3> tD;
654 
655  boost::shared_ptr<MatrixDouble> mGradPtr;
656  boost::shared_ptr<MatrixDouble> mStrainPtr;
657  boost::shared_ptr<MatrixDouble> mStressPtr;
658 
659  boost::shared_ptr<MatrixDouble> mDeformationGradient;
660  boost::shared_ptr<MatrixDouble> matC;
661  boost::shared_ptr<VectorDouble> detF;
662 
663  boost::shared_ptr<MatrixDouble> matEigenVal;
664  boost::shared_ptr<MatrixDouble> matEigenVector;
665 
666  boost::shared_ptr<MatrixDouble> materialTangent;
667  boost::shared_ptr<MatrixDouble> dE_dF_mat;
668 
669  SmartPetscObj<Vec> stateVecPlast;
670  SmartPetscObj<Vec> stateVecCont;
671 
672  boost::shared_ptr<MatrixDouble> contactDispPtr;
673 
675 
676  double timeStepSize;
677 };
678 
679 template <typename T>
680 inline Tensor2<double, 3, 3>
681 getFiniteDiff(const int &kk, const int &ll, CommonData &common_data,
682  Tensor2<T, 3, 3> &F, Tensor2<double, 3, 3> &t_stress) {
684  Tensor2<double, 3, 3> out;
685 
686  Tensor2<double, 3, 3> strainPl;
687  Tensor2<double, 3, 3> strainMin;
688  strainPl(i, j) = F(i, j);
689  strainMin(i, j) = F(i, j);
690  const double h = 1e-8;
691  strainPl(kk, ll) += h;
692  strainMin(kk, ll) -= h;
693  Tensor4<PackPtr<double *, 1>, 3, 3, 3, 3> dummy;
694 
695  auto calculate_stress = [&](auto &F) {
696  Tensor2<double, 3, 3> Piola;
697  Tensor2_symmetric<double, 3> stress_sym;
698  Tensor2_symmetric<double, 3> C;
699  Tensor2_symmetric<double, 3> strain;
700 
701  C(i, j) = F(m, i) ^ F(m, j);
702 
703  Tensor1<double, 3> lambda;
704  Tensor2<double, 3, 3> eig_vec;
705  CHKERR get_eigen_val_and_proj_lapack(C, lambda, eig_vec);
706 
707  auto lnC = get_ln_X_lapack(C, lambda, eig_vec);
708  // auto lnC = get_log_map(C);
709  strain(i, j) = 0.5 * lnC(i, j);
710 
711  // if (with_plasticity)
712  // strain(i, j) -= cur_plastic_strains(i, j); //FIXME:
713 
714  stress_sym(i, j) = common_data.tD(i, j, k, l) * strain(k, l);
715  auto stress_tmp = to_non_symm(stress_sym);
716 
717  // Tensor4<double, 3, 3, 3, 3> D2;
718  // constexpr auto t_kd = Kronecker_Delta<int>();
719  // Tensor4<double, 3, 3, 3, 3> dC_dF;
720  // dC_dF(i, j, k, l) = (t_kd(i, l) * F(k, j)) + (t_kd(j, l) * F(k, i));
721  // auto t_L = get_diff_log_map(C);
722  // t_L(i, j, k, l) *= 0.5;
723  // D2(i, j, k, l) = t_L(i, j, m, n) * dC_dF(m, n, k, l);
724 
725  auto D2 = calculate_nominal_moduli_TL(C, F, stress_sym, lambda, eig_vec,
726  dummy, false);
727 
728  Piola(k, l) = stress_tmp(i, j) * D2(i, j, k, l);
729 
730  return Piola;
731  };
732 
733  auto StressPlusH = calculate_stress(strainPl);
734  // auto StressMinusH = calculate_stress(strainMin);
735  out(i, j) = (StressPlusH(i, j) - t_stress(i, j)) / (h);
736 
737  // out(i, j) = (StressPlusH(i, j) - StressMinusH(i, j)) / (2 * h);
738 
739  return out;
740 }
741 
742 template <typename T>
743 inline Tensor4<double, 3, 3, 3, 3>
744 getDTensorFunction(CommonData &common_data, Tensor2<T, 3, 3> &F,
745  Tensor2<double, 3, 3> &t_stress) {
746  Tensor4<double, 3, 3, 3, 3> my_D;
747  for (int k = 0; k != 3; ++k) {
748  for (int l = 0; l != 3; ++l) {
749  auto d_stress = getFiniteDiff(k, l, common_data, F, t_stress);
750  for (int i = 0; i != 3; ++i) {
751  for (int j = 0; j != 3; ++j) {
752  my_D(i, j, k, l) = d_stress(i, j);
753  }
754  }
755  }
756  }
757  return my_D;
758 };
759 
760 //! [Operators definitions]
761 typedef boost::function<Tensor1<double, 3>(const double, const double,
762  const double)>
763  VectorFun;
764 
765 struct OpStrain : public DomainEleOp {
766  OpStrain(const std::string field_name,
767  boost::shared_ptr<CommonData> common_data_ptr);
768  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
769 
770 private:
771  boost::shared_ptr<CommonData> commonDataPtr;
772 };
773 
774 struct OpLogStrain : public DomainEleOp {
775  OpLogStrain(const std::string field_name,
776  boost::shared_ptr<CommonData> common_data_ptr);
777  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
778 
779 private:
780  boost::shared_ptr<CommonData> commonDataPtr;
781 };
782 
783 struct OpStress : public DomainEleOp {
784  OpStress(const std::string field_name,
785  boost::shared_ptr<CommonData> common_data_ptr);
786  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
787 
788 private:
789  boost::shared_ptr<CommonData> commonDataPtr;
790 };
791 
793  OpCalculateLogStressTangent(const std::string field_name,
794  boost::shared_ptr<CommonData> common_data_ptr);
795  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
796 
797 private:
798  boost::shared_ptr<CommonData> commonDataPtr;
799 };
800 
802  OpRotatingFrameRhs(const std::string field_name,
803  boost::shared_ptr<CommonData> common_data_ptr);
804  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
805 
806 private:
807  boost::shared_ptr<CommonData> commonDataPtr;
808 };
809 
811  OpRotatingFrameBoundaryRhs(const std::string field_name,
812  boost::shared_ptr<CommonData> common_data_ptr);
813  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
814 
815 private:
816  boost::shared_ptr<CommonData> commonDataPtr;
817 };
818 
820  OpRotatingFrameLhs(const std::string row_field_name,
821  const std::string col_field_name,
822  boost::shared_ptr<CommonData> common_data_ptr);
823  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
824  EntityType col_type, EntData &row_data,
825  EntData &col_data);
826 
827 private:
829  boost::shared_ptr<CommonData> commonDataPtr;
830 };
831 
833  OpRotatingFrameBoundaryLhs(const std::string row_field_name,
834  const std::string col_field_name,
835  boost::shared_ptr<CommonData> common_data_ptr);
836  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
837  EntityType col_type, EntData &row_data,
838  EntData &col_data);
839 
840 private:
841  boost::shared_ptr<CommonData> commonDataPtr;
843 };
844 
845 template <bool LOGSTRAIN, bool REACTIONS>
846 struct OpInternalForceRhs : public DomainEleOp {
847  OpInternalForceRhs(const std::string field_name,
848  boost::shared_ptr<CommonData> common_data_ptr,
849  MoFEM::Interface &m_field);
850  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
851 
852 private:
853  boost::shared_ptr<CommonData> commonDataPtr;
855 };
856 
857 struct OpForceRhs : public DomainEleOp {
858  OpForceRhs(const std::string field_name,
859  boost::shared_ptr<CommonData> common_data_ptr,
860  VectorFun body_force);
861  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
862 
863 private:
865  boost::shared_ptr<CommonData> commonDataPtr;
866 };
867 
869  OpCentrifugalForce(const std::string field_name,
870  boost::shared_ptr<CommonData> common_data_ptr);
871  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
872 
873 private:
874  boost::shared_ptr<CommonData> commonDataPtr;
875 };
876 
877 struct OpStiffnessMatrixLhs : public DomainEleOp {
878  OpStiffnessMatrixLhs(const std::string row_field_name,
879  const std::string col_field_name,
880  boost::shared_ptr<CommonData> common_data_ptr);
881  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
882  EntityType col_type, EntData &row_data,
883  EntData &col_data);
884 
885 private:
887  boost::shared_ptr<CommonData> commonDataPtr;
888 };
889 
891  OpLogStrainMatrixLhs(const std::string row_field_name,
892  const std::string col_field_name,
893  boost::shared_ptr<CommonData> common_data_ptr);
894  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
895  EntityType col_type, EntData &row_data,
896  EntData &col_data);
897 
898 private:
900  boost::shared_ptr<CommonData> commonDataPtr;
901 };
902 
903 template <bool LOGSTRAIN> struct OpPostProcElastic : public DomainEleOp {
904  OpPostProcElastic(const std::string field_name,
905  moab::Interface &post_proc_mesh,
906  std::vector<EntityHandle> &map_gauss_pts,
907  boost::shared_ptr<CommonData> common_data_ptr);
908  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
909 
910 private:
912  std::vector<EntityHandle> &mapGaussPts;
913  boost::shared_ptr<CommonData> commonDataPtr;
914 };
915 //! [Operators definitions]
916 
918  OpSaveReactionForces(MoFEM::Interface &m_field, const std::string field_name,
919  moab::Interface &post_proc_mesh,
920  std::vector<EntityHandle> &map_gauss_pts,
921  boost::shared_ptr<CommonData> common_data_ptr);
922  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
923 
924 private:
927  std::vector<EntityHandle> &mapGaussPts;
928  boost::shared_ptr<CommonData> commonDataPtr;
929 };
930 
932  OpBoundaryMassMatrix_UU(std::string disp_name, size_t comp,
933  const Range &essential_ents)
934  : BoundaryEleOp(disp_name, disp_name, BoundaryEleOp::OPROWCOL),
935  cOmp(comp), essentialEnts(essential_ents) {
936  sYmm = true;
937  }
938 
939  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
940  EntityType col_type, EntData &row_data,
941  EntData &col_data) {
943  const int nb_row_dofs = row_data.getIndices().size();
944  const int nb_col_dofs = col_data.getIndices().size();
945 
946  if (nb_row_dofs && nb_col_dofs) {
947 
949  if (essentialEnts.find(ent) == essentialEnts.end()) {
951  }
952 
953  auto get_tensor3 = [](MatrixDouble &m, const int r, const int c) {
955  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
956  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
957  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
958  };
959 
960  locMat.resize(nb_row_dofs, nb_col_dofs, false);
961  locMat.clear();
962 
963  const int nb_integration_pts = getGaussPts().size2();
964 
965  auto t_row_val = row_data.getFTensor0N();
966 
967  auto t_w = getFTensor0IntegrationWeight();
968  const double vol = getMeasure();
969 
970  for (int gg = 0; gg != nb_integration_pts; ++gg) {
971  const double a = vol * t_w;
972 
973  for (int rr = 0; rr != nb_row_dofs / 3; ++rr) {
974 
975  auto t_col_val = col_data.getFTensor0N(gg, 0);
976 
977  for (int cc = 0; cc != nb_col_dofs / 3; ++cc) {
978  auto t_subLocMat = get_tensor3(locMat, 3 * rr, 3 * cc);
979  t_subLocMat(cOmp, cOmp) += (a * t_row_val * t_col_val);
980 
981  ++t_col_val;
982  }
983  ++t_row_val;
984  }
985  ++t_w;
986  }
987 
988  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data,
989  &locMat(0, 0), ADD_VALUES);
990 
991  if (row_side != col_side || row_type != col_type) {
992  locTransMat.resize(nb_col_dofs, nb_row_dofs, false);
993  noalias(locTransMat) = trans(locMat);
994  CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
995  &locTransMat(0, 0), ADD_VALUES);
996  }
997  }
999  }
1000 
1001 private:
1004  size_t cOmp;
1005 };
1006 
1008  OpEssentialRHS_U(std::string disp_name,
1009  boost::shared_ptr<CommonData> common_data_ptr, size_t comp,
1010  VectorDouble &disp, const Range &essential_ents)
1011  : BoundaryEleOp(disp_name, BoundaryEleOp::OPROW),
1012  essentialEnts(essential_ents), cOmp(comp), disP(disp),
1013  commonDataPtr(common_data_ptr) {}
1014 
1015  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data) {
1017  const int nb_row_dofs = row_data.getIndices().size();
1018 
1019  if (nb_row_dofs) {
1020 
1022  if (essentialEnts.find(ent) == essentialEnts.end()) {
1024  }
1025 
1026  auto get_tensor1 = [](VectorDouble &m, const int r) {
1027  return FTensor::Tensor1<double *, 3>(&m(r + 0), &m(r + 1), &m(r + 2));
1028  };
1029 
1030  auto t_dispValue = getFTensor1FromMat<3>(*commonDataPtr->contactDispPtr);
1031 
1032  locRes.resize(nb_row_dofs, false);
1033  locRes.clear();
1034 
1035  const int nb_integration_pts = getGaussPts().size2();
1036  auto t_row_val = row_data.getFTensor0N();
1037  auto t_w = getFTensor0IntegrationWeight();
1038  const double vol = getMeasure();
1039 
1040  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1041  const double a = vol * t_w;
1042 
1043  for (int rr = 0; rr != nb_row_dofs / 3; ++rr) {
1044  auto t_subLocRes = get_tensor1(locRes, 3 * rr);
1045  t_subLocRes(cOmp) += a * t_row_val * (t_dispValue(cOmp) - disP(0));
1046  ++t_row_val;
1047  }
1048  ++t_dispValue;
1049  ++t_w;
1050  }
1051 
1052  CHKERR VecSetValues(getFEMethod()->ts_F, row_data, &*locRes.begin(),
1053  ADD_VALUES);
1054  }
1056  }
1057 
1058 private:
1060  size_t cOmp;
1061  boost::shared_ptr<CommonData> commonDataPtr;
1064 };
1065 
1066 }; // namespace OpElasticTools
const double T
auto get_ln_X_lapack(const Tensor2_symmetric< T1, 3 > &X, Tensor1< T2, 3 > &lambda, Tensor2< T3, 3, 3 > &eig_vec)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
MoFEMErrorCode calcEigenvalues3x3(const Tensor2< double, 3, 3 > &Tensor, double *eigen)
boost::shared_ptr< CommonData > commonDataPtr
OpPostProcElastic(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
[Stiffness]
Definition: ElasticOps.hpp:327
moab::Interface & postProcMesh
Definition: ElasticOps.hpp:94
Deprecated interface functions.
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDeformationGradient
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:571
boost::shared_ptr< VectorDouble > detF
auto get_diff_log_map(Tensor2_symmetric< T, 3 > &X)
boost::shared_ptr< CommonData > commonDataPtr
SmartPetscObj< Vec > stateVecPlast
boost::shared_ptr< MatrixDouble > matEigenVal
auto calculate_nominal_moduli_TL(Tensor2_symmetric< T, 3 > &X, Tensor2< TP, 3, 3 > &F, Tensor2_symmetric< TP, 3 > &stress, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec, Tensor4< T3, 3, 3, 3, 3 > &TLs3, bool compute_tangent)
OpStrain(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Operators definitions]
Definition: ElasticOps.hpp:100
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
FTensor::Index< 'j', 2 > j
Definition: ElasticOps.hpp:27
OpSaveReactionForces(MoFEM::Interface &m_field, const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
auto get_log_map(const Tensor2_symmetric< T, 3 > &X)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
[Stiffness]
auto get_ln_X(const Tensor2_symmetric< T, 3 > &X)
OpForceRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, VectorFun body_force)
[Internal force]
Definition: ElasticOps.hpp:205
boost::shared_ptr< MatrixDouble > dE_dF_mat
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
Kronecker Delta class.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpInternalForceRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Calculate stress]
Definition: ElasticOps.hpp:152
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
Definition: ElasticOps.hpp:338
boost::shared_ptr< MatrixDouble > contactDispPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpLogStrain(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: ElasticOps.hpp:22
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
static Index< 'l', 3 > l
OpLogStrainMatrixLhs(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
auto to_symm(T &non_symm)
OpEssentialRHS_U(std::string disp_name, boost::shared_ptr< CommonData > common_data_ptr, size_t comp, VectorDouble &disp, const Range &essential_ents)
OpRotatingFrameBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
OpCalculateLogStressTangent(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
SmartPetscObj< Vec > stateVecCont
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
Definition: ElasticOps.hpp:134
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: ElasticOps.hpp:20
OpCentrifugalForce(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode get_eigen_val_and_proj_lapack(const Tensor2_symmetric< T, 3 > &X, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec)
Tensor4< double, 3, 3, 3, 3 > getDTensorFunction(CommonData &common_data, Tensor2< T, 3, 3 > &F, Tensor2< double, 3, 3 > &t_stress)
std::vector< EntityHandle > & mapGaussPts
Definition: ElasticOps.hpp:95
static Index< 'n', 3 > n
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:552
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate strain]
OpBoundaryMassMatrix_UU(std::string disp_name, size_t comp, const Range &essential_ents)
bool is_diff(const double &a, const double &b)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:273
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
const int dim
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode calcDerivativeLog3x3(const Tensor2< double, 3, 3 > &X, double *eigen, Tensor2< double, 3, 3 > *eigenProj, Tensor4< double, 3, 3, 3, 3 > &D_X)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
static Index< 'm', 3 > m
boost::shared_ptr< MatrixDouble > matC
bool sYmm
If true assume that matrix is symmetric structure.
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
Ddg< double, 3, 3 > tensor4_to_ddg(Tensor4< T, 3, 3, 3, 3 > &t)
MoFEMErrorCode get_eigen_val_and_proj_souza(const Tensor2_symmetric< T, 3 > &X, double *e, Tensor2< double, 3, 3 > *E)
auto calculate_d_05lnX_dX(Tensor2_symmetric< T, 3 > &X)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
auto get_rotation_R(Tensor1< double, 3 > t1_omega, double tt)
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Body force]
boost::function< FTensor::Tensor1< double, 2 >const double, const double)> VectorFun
[Operators definitions]
Definition: ElasticOps.hpp:33
static Number< 2 > N2
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:41
std::vector< EntityHandle > & mapGaussPts
FTensor::Index< 'k', 2 > k
Definition: ElasticOps.hpp:28
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:70
bool is_eq(const double &a, const double &b)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
OpRotatingFrameLhs(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Tensor2< double, 3, 3 > getFiniteDiff(const int &kk, const int &ll, CommonData &common_data, Tensor2< T, 3, 3 > &F, Tensor2< double, 3, 3 > &t_stress)
OpStiffnessMatrixLhs(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Body force]
Definition: ElasticOps.hpp:262
auto getFTensor0IntegrationWeight()
Get integration weights.
static Index< 'j', 3 > j
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:83
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:50
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ElasticOps.hpp:26
OpRotatingFrameBoundaryLhs(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Body force]
Definition: ElasticOps.hpp:212
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
const double r
rate factor
boost::shared_ptr< CommonData > commonDataPtr
static Number< 1 > N1
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate strain]
Definition: ElasticOps.hpp:109
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: ElasticOps.hpp:21
const double R
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
auto exp_map(const Tensor2_symmetric< T, 3 > &X)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > matEigenVector
bool is_ddg(Tensor4< T, 3, 3, 3, 3 > tensor)
FTensor::Ddg< double, 2, 2 > tD
Definition: ElasticOps.hpp:19
static Index< 'k', 3 > k
boost::shared_ptr< CommonData > commonDataPtr
FTensor::Index< 'l', 2 > l
Definition: ElasticOps.hpp:29
Definition: ElasticOps.hpp:86
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > materialTangent
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
Data on single entity (This is passed as argument to DataOperator::doWork)
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:96
auto to_non_symm(T &symm)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
constexpr double third
auto get_tensor4_from_ddg(Ddg< T, 3, 3 > &ddg)
static const double eps
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
static Number< 0 > N0
OpRotatingFrameRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:59
MoFEMErrorCode calcEigenproj3x3(const Tensor2< double, 3, 3 > &Tensor, double *eigen, Tensor2< double, 3, 3 > *eigenProj)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Internal force]
Definition: ElasticOps.hpp:158
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
[Stiffness]
Definition: ElasticOps.hpp:271
OpStress(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Calculate strain]
Definition: ElasticOps.hpp:125