v0.14.0
EshelbianOperators.cpp
Go to the documentation of this file.
1 /**
2  * \file EshelbianOperators.cpp
3  * \example EshelbianOperators.cpp
4  *
5  * \brief Implementation of operators
6  */
7 
8 #include <MoFEM.hpp>
9 using namespace MoFEM;
10 
11 #include <EshelbianPlasticity.hpp>
12 #include <boost/math/constants/constants.hpp>
13 
14 #include <lapack_wrap.h>
15 
16 #include <MatrixFunction.hpp>
17 
18 namespace EshelbianPlasticity {
19 
20 inline auto diff_deviator(FTensor::Ddg<double, 3, 3> &&t_diff_stress) {
25 
26  FTensor::Ddg<double, 3, 3> t_diff_deviator;
27  t_diff_deviator(i, j, k, l) = t_diff_stress(i, j, k, l);
28 
29  constexpr double third = boost::math::constants::third<double>();
30 
31  t_diff_deviator(0, 0, 0, 0) -= third;
32  t_diff_deviator(0, 0, 1, 1) -= third;
33 
34  t_diff_deviator(1, 1, 0, 0) -= third;
35  t_diff_deviator(1, 1, 1, 1) -= third;
36 
37  t_diff_deviator(2, 2, 0, 0) -= third;
38  t_diff_deviator(2, 2, 1, 1) -= third;
39 
40  t_diff_deviator(0, 0, 2, 2) -= third;
41  t_diff_deviator(1, 1, 2, 2) -= third;
42  t_diff_deviator(2, 2, 2, 2) -= third;
43 
44  return t_diff_deviator;
45 }
46 
47 constexpr static auto size_symm = (3 * (3 + 1)) / 2;
48 
49 inline auto diff_tensor() {
56  t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
57  return t_diff;
58 };
59 
60 inline auto symm_L_tensor() {
65  t_L(i, j, L) = 0;
66  t_L(0, 0, 0) = 1;
67  t_L(1, 1, 3) = 1;
68  t_L(2, 2, 5) = 1;
69  t_L(1, 0, 1) = 1;
70  t_L(2, 0, 2) = 1;
71  t_L(2, 1, 4) = 1;
72  return t_L;
73 }
74 
75 template <class T> struct TensorTypeExtractor {
77 };
78 template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
80 };
81 
82 template <typename T>
84  RotSelector rotSelector = LARGE_ROT) {
85 
86  using D = typename TensorTypeExtractor<T>::Type;
87 
92 
93  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
94  t_R(i, j) = t_kd(i, j);
95 
97  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
98 
99  if (rotSelector == SMALL_ROT) {
100  t_R(i, j) += t_Omega(i, j);
101  return t_R;
102  }
103 
104  const auto angle =
105  sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
106  const auto a = sin(angle) / angle;
107  t_R(i, j) += a * t_Omega(i, j);
108  if (rotSelector == MODERATE_ROT)
109  return t_R;
110 
111  const auto ss_2 = sin(angle / 2.) / angle;
112  const auto b = 2. * ss_2 * ss_2;
113  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
114 
115  return t_R;
116 }
117 
118 template <typename T>
120  RotSelector rotSelector = LARGE_ROT) {
121 
122  using D = typename TensorTypeExtractor<T>::Type;
123 
128 
130 
131  const auto angle =
132  sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
133  if (rotSelector == SMALL_ROT) {
134  t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, j, k);
135  return t_diff_R;
136  }
137 
138  const auto ss = sin(angle);
139  const auto a = ss / angle;
140  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
141 
143  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
144 
145  const auto angle2 = angle * angle;
146  const auto cc = cos(angle);
147  const auto diff_a = (angle * cc - ss) / angle2;
148  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
149  if (rotSelector == MODERATE_ROT)
150  return t_diff_R;
151 
152  const auto ss_2 = sin(angle / 2.);
153  const auto cc_2 = cos(angle / 2.);
154  const auto b = 2. * ss_2 * ss_2 / angle2;
155  t_diff_R(i, j, k) +=
156  b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
157  FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
158  const auto diff_b =
159  (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
160  t_diff_R(i, j, k) +=
161  diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
162 
163  return t_diff_R;
164 }
165 
166 template <typename T>
168  RotSelector rotSelector = LARGE_ROT) {
169 
170  using D = typename TensorTypeExtractor<T>::Type;
171 
174 
175  constexpr double eps = 1e-10;
176  for (int l = 0; l != 3; ++l) {
178  t_omega_c(i) = t_omega(i);
179  t_omega_c(l) += std::complex<double>(0, eps);
180  auto t_diff_R_c = get_diff_rotation_form_vector(t_omega_c, rotSelector);
181  for (int i = 0; i != 3; ++i) {
182  for (int j = 0; j != 3; ++j) {
183  for (int k = 0; k != 3; ++k) {
184  t_diff2_R(i, j, k, l) = t_diff_R_c(i, j, k).imag() / eps;
185  }
186  }
187  }
188  }
189 
190  return t_diff2_R;
191 }
192 
193 struct isEq {
194  static inline auto check(const double &a, const double &b) {
195  constexpr double eps = std::numeric_limits<float>::epsilon();
196  return std::abs(a - b) / absMax < eps;
197  }
198  static double absMax;
199 };
200 
201 double isEq::absMax = 1;
202 
203 inline auto is_eq(const double &a, const double &b) {
204  return isEq::check(a, b);
205 };
206 
207 template <int DIM> inline auto get_uniq_nb(double *ptr) {
208  std::array<double, DIM> tmp;
209  std::copy(ptr, &ptr[DIM], tmp.begin());
210  std::sort(tmp.begin(), tmp.end());
211  isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
212  constexpr double eps = std::numeric_limits<float>::epsilon();
213  isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
214  return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
215 }
216 
217 template <int DIM>
220 
221  isEq::absMax =
222  std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
223 
224  int i = 0, j = 1, k = 2;
225 
226  if (is_eq(eig(0), eig(1))) {
227  i = 0;
228  j = 2;
229  k = 1;
230  } else if (is_eq(eig(0), eig(2))) {
231  i = 0;
232  j = 1;
233  k = 2;
234  } else if (is_eq(eig(1), eig(2))) {
235  i = 1;
236  j = 0;
237  k = 2;
238  }
239 
240  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
241  eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
242 
243  eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
244 
245  eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
246 
247  FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
248 
249  {
252  eig(i) = eig_c(i);
253  eigen_vec(i, j) = eigen_vec_c(i, j);
254  }
255 }
256 
257 MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
258  EntityType type,
259  EntData &data) {
261 
263  auto t_L = symm_L_tensor();
264 
265  int nb_integration_pts = getGaussPts().size2();
272 
273  dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
274  dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts, false);
275  dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts, false);
276  dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts, false);
277  dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts, false);
278  dataAtPts->leviKirchhoffOmegaAtPts.resize(9, nb_integration_pts, false);
279 
280  dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts, false);
281  dataAtPts->adjointPdUAtPts.resize(size_symm, nb_integration_pts, false);
282  dataAtPts->adjointPdUdPAtPts.resize(9 * size_symm, nb_integration_pts, false);
283  dataAtPts->adjointPdUdOmegaAtPts.resize(3 * size_symm, nb_integration_pts,
284  false);
285  dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
286  dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
287  dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts, false);
288  dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts, false);
289  dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
290  dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
291  dataAtPts->nbUniq.resize(nb_integration_pts, false);
292 
293  dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
294 
295  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
296  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
297  auto t_h_dlog_u =
298  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
299  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
300  auto t_levi_kirchoff_dP =
301  getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
302  auto t_levi_kirchoff_domega =
303  getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
304  auto t_approx_P_adjont_dstretch =
305  getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
306  auto t_approx_P_adjont_log_du =
307  getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
308  auto t_approx_P_adjont_log_du_dP =
309  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
310  auto t_approx_P_adjont_log_du_domega =
311  getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
312  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
313  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
314  auto t_diff_R = getFTensor3FromMat<3, 3, 3>(dataAtPts->diffRotMatAtPts);
315  auto t_log_u =
316  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
317  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
318  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
319  auto t_diff_u =
320  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
321  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
322  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
323 
324  auto &nbUniq = dataAtPts->nbUniq;
325  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
326  auto t_log_stretch_total =
327  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
329 
330  for (int gg = 0; gg != nb_integration_pts; ++gg) {
331 
337  FTensor::Tensor2<double, 3, 3> t_approx_P_intermidiate;
338 
339  t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
340 
341  auto calulate_rotation = [&]() {
342  auto t0_diff =
343  get_diff_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
344  auto t0 = get_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
345  t_diff_R(i, j, k) = t0_diff(i, j, k);
346  t_R(i, j) = t0(i, j);
347  };
348 
349  auto calulate_streach = [&]() {
350  eigen_vec(i, j) = t_log_u(i, j);
351  CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
352  // rare case when two eigen values are equal
353  nbUniq[gg] = get_uniq_nb<3>(&eig(0));
354  if (nbUniq[gg] < 3) {
355  sort_eigen_vals<3>(eig, eigen_vec);
356  }
357  t_eigen_vals(i) = eig(i);
358  t_eigen_vecs(i, j) = eigen_vec(i, j);
359  auto t_u_tmp =
360  EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f);
361  t_u(i, j) = t_u_tmp(i, j);
362  auto t_diff_u_tmp =
363  EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f,
364  EshelbianCore::d_f, nbUniq[gg]);
365  t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
366  t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
367  };
368 
369  calulate_rotation();
370  calulate_streach();
371 
372  switch (EshelbianCore::gradApperoximator) {
373  case LARGE_ROT:
374 
375  // rotation can be large, moderate or small
376  t_Ru(i, m) = t_R(i, l) * t_u(l, m);
377  t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
378  t_h_domega(i, j, k) = (t_diff_R(i, l, k) * t_u(l, m)) * t_h1(m, j);
379  t_h_dlog_u(i, j, L) = (t_R(i, l) * t_Ldiff_u(l, m, L)) * t_h1(m, j);
380 
381  // conservation of angular momentum at current configuration
382  t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
383  t_approx_P_adjont_dstretch(l, m) =
384  t_approx_P_intermidiate(i, m) * t_R(i, l);
385 
386  t_levi_kirchoff(k) =
387  levi_civita(l, m, k) * (t_approx_P_adjont_dstretch(l, m));
388  t_levi_kirchoff_dP(i, j, k) =
389  (levi_civita(l, m, k) * t_R(i, l)) * t_h1(m, j);
390  t_levi_kirchoff_domega(k, n) =
391  levi_civita(l, m, k) *
392  (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
393 
394  t_approx_P_adjont_log_du(L) =
395  t_Ldiff_u(l, m, L) * t_approx_P_adjont_dstretch(l, m);
396  t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
397  t_approx_P_adjont_log_du_domega(n, L) =
398  t_Ldiff_u(l, m, L) *
399  (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
400 
401  break;
402  case MODERATE_ROT:
403 
404  {
405 
406  // assume small current rotation
408  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
409  t_Ru(i, m) = t_Omega(i, m) + t_u(i, m);
410  t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
411  t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, m, k) * t_h1(m, j);
412  t_h_dlog_u(i, j, L) = t_Ldiff_u(i, m, L) * t_h1(m, j);
413 
414  // conservation of angular momentum at intermediate configuration
415  t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
416  t_approx_P_adjont_dstretch(i, m) = t_approx_P_intermidiate(i, m);
417 
418  t_levi_kirchoff(k) =
419  levi_civita(i, m, k) * (t_approx_P_adjont_dstretch(i, m));
420  t_levi_kirchoff_dP(i, j, k) = levi_civita(i, m, k) * t_h1(m, j);
421  t_levi_kirchoff_domega(k, n) = 0;
422 
423  t_approx_P_adjont_log_du(L) =
424  t_Ldiff_u(i, m, L) * t_approx_P_adjont_dstretch(i, m);
425  t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
426  t_approx_P_adjont_log_du_domega(n, L) = 0;
427 
428  }
429 
430  break;
431 
432  case SMALL_ROT:
433 
434  {
435 
436  // assume small current rotation
438  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
439  t_h(i, j) = t_Omega(i, j) + t_u(i, j);
440  t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, j, k);
441  t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
442 
443  // conservation of angular momentum at reference condition
444  t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
445  t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
446  t_levi_kirchoff_domega(k, l) = 0;
447 
448  t_approx_P_adjont_dstretch(i, j) = t_approx_P(i, j);
449  t_approx_P_adjont_log_du(L) =
450  t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
451  t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
452  t_approx_P_adjont_log_du_domega(m, L) = 0;
453 
454  }
455 
456  break;
457  }
458 
460  t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
461  eigen_vec(i, j) = t_C_h1(i, j);
462  CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
463  switch (EshelbianCore::stretchSelector) {
465  break;
467  for (int ii = 0; ii != 3; ++ii) {
468  eig(ii) = std::max(eig(ii),
469  2 * EshelbianCore::inv_f(EshelbianCore::f(
470  std::numeric_limits<double>::min_exponent)));
471  }
472  break;
473  }
474 
475  auto t_log_u2_h1_tmp =
476  EigenMatrix::getMat(eig, eigen_vec, EshelbianCore::inv_f);
477 
478  switch (EshelbianCore::gradApperoximator) {
479  case LARGE_ROT:
480  case MODERATE_ROT:
481  t_log_stretch_total(i, j) = t_log_u2_h1_tmp(i, j) / 2 + t_log_u(i, j);
482  break;
483  case SMALL_ROT:
484  t_log_stretch_total(i, j) = t_log_u(i, j);
485  break;
486  };
487 
488  ++t_h;
489  ++t_h_domega;
490  ++t_h_dlog_u;
491  ++t_levi_kirchoff;
492  ++t_levi_kirchoff_dP;
493  ++t_levi_kirchoff_domega;
494  ++t_approx_P_adjont_dstretch;
495  ++t_approx_P_adjont_log_du;
496  ++t_approx_P_adjont_log_du_dP;
497  ++t_approx_P_adjont_log_du_domega;
498  ++t_approx_P;
499  ++t_R;
500  ++t_diff_R;
501  ++t_log_u;
502  ++t_u;
503  ++t_diff_u;
504  ++t_eigen_vals;
505  ++t_eigen_vecs;
506  ++t_omega;
507  ++t_grad_h1;
508  ++t_log_stretch_total;
509  }
510 
512 }
513 
514 MoFEMErrorCode OpSpatialEquilibrium::integrate(EntData &data) {
516  int nb_dofs = data.getIndices().size();
517  int nb_integration_pts = data.getN().size1();
518  auto v = getVolume();
519  auto t_w = getFTensor0IntegrationWeight();
520  auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
521  auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
522  if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
523  dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
524  dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
525  dataAtPts->wL2DotDotAtPts.clear();
526  }
527  auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
528 
529  int nb_base_functions = data.getN().size2();
530  auto t_row_base_fun = data.getFTensor0N();
531 
533  auto get_ftensor1 = [](auto &v) {
535  &v[2]);
536  };
537 
538  for (int gg = 0; gg != nb_integration_pts; ++gg) {
539  double a = v * t_w;
540  auto t_nf = get_ftensor1(nF);
541  int bb = 0;
542  for (; bb != nb_dofs / 3; ++bb) {
543  t_nf(i) += a * t_row_base_fun * t_div_P(i);
544  t_nf(i) -= a * t_row_base_fun * alphaW * t_s_dot_w(i);
545  t_nf(i) -= a * t_row_base_fun * alphaRho * t_s_dot_dot_w(i);
546  ++t_nf;
547  ++t_row_base_fun;
548  }
549  for (; bb != nb_base_functions; ++bb)
550  ++t_row_base_fun;
551  ++t_w;
552  ++t_div_P;
553  ++t_s_dot_w;
554  ++t_s_dot_dot_w;
555  }
556 
558 }
559 
560 MoFEMErrorCode OpSpatialRotation::integrate(EntData &data) {
562  int nb_dofs = data.getIndices().size();
563  int nb_integration_pts = getGaussPts().size2();
564  auto v = getVolume();
565  auto t_w = getFTensor0IntegrationWeight();
566  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
567  int nb_base_functions = data.getN().size2();
568  auto t_row_base_fun = data.getFTensor0N();
572  auto get_ftensor1 = [](auto &v) {
574  &v[2]);
575  };
576  for (int gg = 0; gg != nb_integration_pts; ++gg) {
577  double a = v * t_w;
578  auto t_nf = get_ftensor1(nF);
579  int bb = 0;
580  for (; bb != nb_dofs / 3; ++bb) {
581  t_nf(k) += (a * t_row_base_fun) * t_levi_kirchoff(k);
582  ++t_nf;
583  ++t_row_base_fun;
584  }
585  for (; bb != nb_base_functions; ++bb) {
586  ++t_row_base_fun;
587  }
588  ++t_w;
589  ++t_levi_kirchoff;
590  }
592 }
593 
594 MoFEMErrorCode OpSpatialPhysical::integrate(EntData &data) {
596 
597  if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
598  CHKERR integratePiola(data);
599  } else {
600  if (polyConvex) {
601  CHKERR integratePolyconvexHencky(data);
602  } else {
603  CHKERR integrateHencky(data);
604  }
605  }
606 
608 }
609 
610 MoFEMErrorCode OpSpatialPhysical::integrateHencky(EntData &data) {
612 
614  auto t_L = symm_L_tensor();
615 
616  int nb_dofs = data.getIndices().size();
617  int nb_integration_pts = data.getN().size1();
618  auto v = getVolume();
619  auto t_w = getFTensor0IntegrationWeight();
620  auto t_approx_P_adjont_log_du =
621  getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
622  auto t_log_streach_h1 =
623  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
624  auto t_dot_log_u =
625  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
626 
627  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
628 
633  auto get_ftensor2 = [](auto &v) {
635  &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
636  };
637 
638  int nb_base_functions = data.getN().size2();
639  auto t_row_base_fun = data.getFTensor0N();
640  for (int gg = 0; gg != nb_integration_pts; ++gg) {
641  double a = v * t_w;
642  auto t_nf = get_ftensor2(nF);
643 
645  t_T(i, j) =
646  t_D(i, j, k, l) * (t_log_streach_h1(k, l) + alphaU * t_dot_log_u(k, l));
648  t_residual(L) =
649  a * (t_approx_P_adjont_log_du(L) - t_L(i, j, L) * t_T(i, j));
650 
651  int bb = 0;
652  for (; bb != nb_dofs / 6; ++bb) {
653  t_nf(L) += t_row_base_fun * t_residual(L);
654  ++t_nf;
655  ++t_row_base_fun;
656  }
657  for (; bb != nb_base_functions; ++bb)
658  ++t_row_base_fun;
659 
660  ++t_w;
661  ++t_approx_P_adjont_log_du;
662  ++t_dot_log_u;
663  ++t_log_streach_h1;
664  }
666 }
667 
668 MoFEMErrorCode OpSpatialPhysical::integratePolyconvexHencky(EntData &data) {
670 
672  auto t_L = symm_L_tensor();
673 
674  int nb_dofs = data.getIndices().size();
675  int nb_integration_pts = data.getN().size1();
676  auto v = getVolume();
677  auto t_w = getFTensor0IntegrationWeight();
678  auto t_approx_P_adjont_log_du =
679  getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
680  auto t_log_streach_h1 =
681  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
682  auto t_dot_log_u =
683  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
684 
685  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
686 
691  auto get_ftensor2 = [](auto &v) {
693  &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
694  };
695 
696  constexpr double nohat_k = 1. / 4;
697  constexpr double hat_k = 1. / 8;
698  double mu = dataAtPts->mu;
699  double lambda = dataAtPts->lambda;
700 
701  constexpr double third = boost::math::constants::third<double>();
703  auto t_diff_deviator = diff_deviator(diff_tensor());
704 
705  int nb_base_functions = data.getN().size2();
706  auto t_row_base_fun = data.getFTensor0N();
707  for (int gg = 0; gg != nb_integration_pts; ++gg) {
708  double a = v * t_w;
709  auto t_nf = get_ftensor2(nF);
710 
711  double log_det = t_log_streach_h1(i, i);
712  double log_det2 = log_det * log_det;
714  t_dev(i, j) = t_log_streach_h1(i, j) - t_kd(i, j) * (third * log_det);
715  double dev_norm2 = t_dev(i, j) * t_dev(i, j);
716 
718  auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
719  auto B = lambda * std::exp(hat_k * log_det2) * log_det;
720  t_T(i, j) =
721 
722  A * (t_dev(k, l) * t_diff_deviator(k, l, i, j))
723 
724  +
725 
726  B * t_kd(i, j)
727 
728  +
729 
730  alphaU * t_D(i, j, k, l) * t_log_streach_h1(k, l);
731 
733  t_residual(L) =
734  a * (t_approx_P_adjont_log_du(L) - t_L(i, j, L) * t_T(i, j));
735 
736  int bb = 0;
737  for (; bb != nb_dofs / 6; ++bb) {
738  t_nf(L) += t_row_base_fun * t_residual(L);
739  ++t_nf;
740  ++t_row_base_fun;
741  }
742  for (; bb != nb_base_functions; ++bb)
743  ++t_row_base_fun;
744 
745  ++t_w;
746  ++t_approx_P_adjont_log_du;
747  ++t_dot_log_u;
748  ++t_log_streach_h1;
749  }
751 }
752 
753 MoFEMErrorCode OpSpatialPhysical::integratePiola(EntData &data) {
755 
757  auto t_L = symm_L_tensor();
758 
759  int nb_dofs = data.getIndices().size();
760  int nb_integration_pts = data.getN().size1();
761  auto v = getVolume();
762  auto t_w = getFTensor0IntegrationWeight();
763  auto t_approx_P_adjont_log_du =
764  getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
765  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
766  auto t_dot_log_u =
767  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
768  auto t_diff_u =
769  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
770 
775  auto get_ftensor2 = [](auto &v) {
777  &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
778  };
779 
780  int nb_base_functions = data.getN().size2();
781  auto t_row_base_fun = data.getFTensor0N();
782  for (int gg = 0; gg != nb_integration_pts; ++gg) {
783  double a = v * t_w;
784  auto t_nf = get_ftensor2(nF);
785 
787  t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
788 
790  t_viscous_P(i, j) =
791  alphaU *
792  t_dot_log_u(i,
793  j); // That is chit, should be split on axiator and deviator
794 
796  t_residual(L) =
797  a * (t_approx_P_adjont_log_du(L) - t_Ldiff_u(i, j, L) * t_P(i, j) -
798  t_L(i, j, L) * t_viscous_P(i, j));
799 
800  int bb = 0;
801  for (; bb != nb_dofs / 6; ++bb) {
802  t_nf(L) += t_row_base_fun * t_residual(L);
803  ++t_nf;
804  ++t_row_base_fun;
805  }
806  for (; bb != nb_base_functions; ++bb)
807  ++t_row_base_fun;
808 
809  ++t_w;
810  ++t_approx_P_adjont_log_du;
811  ++t_P;
812  ++t_dot_log_u;
813  ++t_diff_u;
814  }
816 }
817 
818 MoFEMErrorCode OpSpatialConsistencyP::integrate(EntData &data) {
820  int nb_dofs = data.getIndices().size();
821  int nb_integration_pts = data.getN().size1();
822  auto v = getVolume();
823  auto t_w = getFTensor0IntegrationWeight();
824  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
825  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
826 
827  int nb_base_functions = data.getN().size2() / 3;
828  auto t_row_base_fun = data.getFTensor1N<3>();
833 
834  auto get_ftensor1 = [](auto &v) {
836  &v[2]);
837  };
838 
839  for (int gg = 0; gg != nb_integration_pts; ++gg) {
840  double a = v * t_w;
841  auto t_nf = get_ftensor1(nF);
842 
843  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
845  t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
846 
847  int bb = 0;
848  for (; bb != nb_dofs / 3; ++bb) {
849  t_nf(i) += a * t_row_base_fun(j) * t_residuum(i, j);
850  ++t_nf;
851  ++t_row_base_fun;
852  }
853 
854  for (; bb != nb_base_functions; ++bb)
855  ++t_row_base_fun;
856 
857  ++t_w;
858  ++t_h;
859  ++t_omega_dot;
860  }
861 
863 }
864 
865 MoFEMErrorCode OpSpatialConsistencyBubble::integrate(EntData &data) {
867  int nb_dofs = data.getIndices().size();
868  int nb_integration_pts = data.getN().size1();
869  auto v = getVolume();
870  auto t_w = getFTensor0IntegrationWeight();
871  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
872  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
873 
874  int nb_base_functions = data.getN().size2() / 9;
875  auto t_row_base_fun = data.getFTensor2N<3, 3>();
880 
881  auto get_ftensor0 = [](auto &v) {
883  };
884 
885  for (int gg = 0; gg != nb_integration_pts; ++gg) {
886  double a = v * t_w;
887  auto t_nf = get_ftensor0(nF);
888 
889  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
891  t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
892 
893  int bb = 0;
894  for (; bb != nb_dofs; ++bb) {
895  t_nf += a * t_row_base_fun(i, j) * t_residuum(i, j);
896  ++t_nf;
897  ++t_row_base_fun;
898  }
899  for (; bb != nb_base_functions; ++bb) {
900  ++t_row_base_fun;
901  }
902  ++t_w;
903  ++t_h;
904  ++t_omega_dot;
905  }
906 
908 }
909 
910 MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
912  int nb_dofs = data.getIndices().size();
913  int nb_integration_pts = data.getN().size1();
914  auto v = getVolume();
915  auto t_w = getFTensor0IntegrationWeight();
916  auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
917  int nb_base_functions = data.getN().size2() / 3;
918  auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
920  auto get_ftensor1 = [](auto &v) {
922  &v[2]);
923  };
924 
925  for (int gg = 0; gg != nb_integration_pts; ++gg) {
926  double a = v * t_w;
927  auto t_nf = get_ftensor1(nF);
928  int bb = 0;
929  for (; bb != nb_dofs / 3; ++bb) {
930  double div_row_base = t_row_diff_base_fun(i, i);
931  t_nf(i) += a * div_row_base * t_w_l2(i);
932  ++t_nf;
933  ++t_row_diff_base_fun;
934  }
935  for (; bb != nb_base_functions; ++bb) {
936  ++t_row_diff_base_fun;
937  }
938  ++t_w;
939  ++t_w_l2;
940  }
941 
943 }
944 
945 MoFEMErrorCode OpDispBc::integrate(EntData &data) {
947  // get entity of face
948  EntityHandle fe_ent = getFEEntityHandle();
949  // interate over all boundary data
950  for (auto &bc : (*bcDispPtr)) {
951  // check if finite element entity is part of boundary condition
952  if (bc.faces.find(fe_ent) != bc.faces.end()) {
953  int nb_dofs = data.getIndices().size();
954  int nb_integration_pts = data.getN().size1();
955  auto t_normal = getFTensor1Normal();
956  auto t_w = getFTensor0IntegrationWeight();
957  int nb_base_functions = data.getN().size2() / 3;
958  auto t_row_base_fun = data.getFTensor1N<3>();
959 
962  auto get_ftensor1 = [](auto &v) {
964  &v[2]);
965  };
966 
967  double scale = 1;
968  for (auto &sm : scalingMethodsVec) {
969  scale *= sm->getScale(getFEMethod()->ts_t);
970  }
971 
972  // get bc data
973  FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
974  t_bc_disp(i) *= scale;
975 
976  for (int gg = 0; gg != nb_integration_pts; ++gg) {
978  t_bc_res(i) = t_bc_disp(i);
979  auto t_nf = get_ftensor1(nF);
980  int bb = 0;
981  for (; bb != nb_dofs / 3; ++bb) {
982  t_nf(i) -=
983  t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
984  ++t_nf;
985  ++t_row_base_fun;
986  }
987  for (; bb != nb_base_functions; ++bb)
988  ++t_row_base_fun;
989 
990  ++t_w;
991  }
992  }
993  }
995 }
996 
997 MoFEMErrorCode OpRotationBc::integrate(EntData &data) {
999  // get entity of face
1000  EntityHandle fe_ent = getFEEntityHandle();
1001  // interate over all boundary data
1002  for (auto &bc : (*bcRotPtr)) {
1003  // check if finite element entity is part of boundary condition
1004  if (bc.faces.find(fe_ent) != bc.faces.end()) {
1005  int nb_dofs = data.getIndices().size();
1006  int nb_integration_pts = data.getN().size1();
1007  auto t_normal = getFTensor1Normal();
1008  auto t_w = getFTensor0IntegrationWeight();
1009 
1010  int nb_base_functions = data.getN().size2() / 3;
1011  auto t_row_base_fun = data.getFTensor1N<3>();
1015  auto get_ftensor1 = [](auto &v) {
1017  &v[2]);
1018  };
1019 
1020  // get bc data
1021  FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
1022 
1023  double theta = bc.theta;
1024  theta *= getFEMethod()->ts_t;
1025 
1027  const double a = sqrt(t_normal(i) * t_normal(i));
1028  t_omega(i) = t_normal(i) * (theta / a);
1029  auto t_R = get_rotation_form_vector(t_omega, LARGE_ROT);
1030  auto t_coords = getFTensor1CoordsAtGaussPts();
1031 
1032  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1034  t_delta(i) = t_center(i) - t_coords(i);
1036  t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
1037 
1038  auto t_nf = get_ftensor1(nF);
1039  int bb = 0;
1040  for (; bb != nb_dofs / 3; ++bb) {
1041  t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
1042  ++t_nf;
1043  ++t_row_base_fun;
1044  }
1045  for (; bb != nb_base_functions; ++bb)
1046  ++t_row_base_fun;
1047 
1048  ++t_w;
1049  ++t_coords;
1050  }
1051  }
1052  }
1054 }
1055 
1056 MoFEMErrorCode FeTractionBc::preProcess() {
1058 
1059  struct FaceRule {
1060  int operator()(int p_row, int p_col, int p_data) const {
1061  return 2 * (p_data + 1);
1062  }
1063  };
1064 
1065  if (ts_ctx == CTX_TSSETIFUNCTION) {
1066 
1067  // Loop boundary elements with traction boundary conditions
1070  {HDIV});
1071  fe.getOpPtrVector().push_back(
1072  new OpTractionBc(std::string("P") /* + "_RT"*/, *this));
1073  fe.getRuleHook = FaceRule();
1074  fe.ts_t = ts_t;
1075  CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
1076  fe, nullptr, MF_ZERO,
1077  this->getCacheWeakPtr());
1078  }
1079 
1081 }
1082 
1083 MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type, EntData &data) {
1085 
1088 
1089  auto t_normal = getFTensor1Normal();
1090  const double nrm2 = sqrt(t_normal(i) * t_normal(i));
1091  FTensor::Tensor1<double, 3> t_unit_normal;
1092  t_unit_normal(i) = t_normal(i) / nrm2;
1093  int nb_dofs = data.getFieldData().size();
1094  int nb_integration_pts = data.getN().size1();
1095  int nb_base_functions = data.getN().size2() / 3;
1096  double ts_t = getFEMethod()->ts_t;
1097 
1098  auto integrate_matrix = [&]() {
1100 
1101  auto t_w = getFTensor0IntegrationWeight();
1102  matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
1103  matPP.clear();
1104 
1105  auto t_row_base_fun = data.getFTensor1N<3>();
1106  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1107 
1108  int rr = 0;
1109  for (; rr != nb_dofs / 3; ++rr) {
1110  const double a = t_row_base_fun(i) * t_unit_normal(i);
1111  auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
1112  for (int cc = 0; cc != nb_dofs / 3; ++cc) {
1113  const double b = t_col_base_fun(i) * t_unit_normal(i);
1114  matPP(rr, cc) += t_w * a * b;
1115  ++t_col_base_fun;
1116  }
1117  ++t_row_base_fun;
1118  }
1119 
1120  for (; rr != nb_base_functions; ++rr)
1121  ++t_row_base_fun;
1122 
1123  ++t_w;
1124  }
1125 
1127  };
1128 
1129  auto integrate_rhs = [&](auto &bc) {
1131 
1132  auto t_w = getFTensor0IntegrationWeight();
1133  vecPv.resize(3, nb_dofs / 3, false);
1134  vecPv.clear();
1135 
1136  auto t_row_base_fun = data.getFTensor1N<3>();
1137  double ts_t = getFEMethod()->ts_t;
1138 
1139  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1140  int rr = 0;
1141  for (; rr != nb_dofs / 3; ++rr) {
1142  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
1143  for (int dd = 0; dd != 3; ++dd)
1144  if (bc.flags[dd])
1145  vecPv(dd, rr) += t * bc.vals[dd];
1146  ++t_row_base_fun;
1147  }
1148  for (; rr != nb_base_functions; ++rr)
1149  ++t_row_base_fun;
1150  ++t_w;
1151  }
1153  };
1154 
1155  auto integrate_rhs_cook = [&](auto &bc) {
1157 
1158  vecPv.resize(3, nb_dofs / 3, false);
1159  vecPv.clear();
1160 
1161  auto t_w = getFTensor0IntegrationWeight();
1162  auto t_row_base_fun = data.getFTensor1N<3>();
1163  auto t_coords = getFTensor1CoordsAtGaussPts();
1164 
1165  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1166 
1167  auto calc_tau = [](double y) {
1168  y -= 44;
1169  y /= (60 - 44);
1170  return -y * (y - 1) / 0.25;
1171  };
1172 
1173  const double tau = calc_tau(t_coords(1));
1174 
1175  int rr = 0;
1176  for (; rr != nb_dofs / 3; ++rr) {
1177  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
1178  for (int dd = 0; dd != 3; ++dd)
1179  if (bc.flags[dd])
1180  vecPv(dd, rr) += t * tau * bc.vals[dd];
1181  ++t_row_base_fun;
1182  }
1183 
1184  for (; rr != nb_base_functions; ++rr)
1185  ++t_row_base_fun;
1186  ++t_w;
1187  ++t_coords;
1188  }
1190  };
1191 
1192  // get entity of face
1193  EntityHandle fe_ent = getFEEntityHandle();
1194  for (auto &bc : *(bcFe.bcData)) {
1195  if (bc.faces.find(fe_ent) != bc.faces.end()) {
1196 
1197  int nb_dofs = data.getFieldData().size();
1198  if (nb_dofs) {
1199 
1200  CHKERR integrate_matrix();
1201  if (std::regex_match(bc.blockName, std::regex(".*COOK.*")))
1202  CHKERR integrate_rhs_cook(bc);
1203  else
1204  CHKERR integrate_rhs(bc);
1205 
1206  const auto info =
1207  lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
1208  nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
1209  if (info > 0)
1210  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1211  "The leading minor of order %i is "
1212  "not positive; definite;\nthe "
1213  "solution could not be computed",
1214  info);
1215 
1216  for (int dd = 0; dd != 3; ++dd)
1217  if (bc.flags[dd])
1218  for (int rr = 0; rr != nb_dofs / 3; ++rr)
1219  data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
1220  }
1221  }
1222  }
1223 
1225 }
1226 
1227 MoFEMErrorCode OpSpatialEquilibrium_dw_dP::integrate(EntData &row_data,
1228  EntData &col_data) {
1230  int nb_integration_pts = row_data.getN().size1();
1231  int row_nb_dofs = row_data.getIndices().size();
1232  int col_nb_dofs = col_data.getIndices().size();
1233  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1235  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
1236  };
1238  auto v = getVolume();
1239  auto t_w = getFTensor0IntegrationWeight();
1240  int row_nb_base_functions = row_data.getN().size2();
1241  auto t_row_base_fun = row_data.getFTensor0N();
1242  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1243  double a = v * t_w;
1244  int rr = 0;
1245  for (; rr != row_nb_dofs / 3; ++rr) {
1246  auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
1247  auto t_m = get_ftensor1(K, 3 * rr, 0);
1248  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1249  double div_col_base = t_col_diff_base_fun(i, i);
1250  t_m(i) += a * t_row_base_fun * div_col_base;
1251  ++t_m;
1252  ++t_col_diff_base_fun;
1253  }
1254  ++t_row_base_fun;
1255  }
1256  for (; rr != row_nb_base_functions; ++rr)
1257  ++t_row_base_fun;
1258  ++t_w;
1259  }
1261 }
1262 
1263 MoFEMErrorCode OpSpatialEquilibrium_dw_dw::integrate(EntData &row_data,
1264  EntData &col_data) {
1266 
1267  if (alphaW < std::numeric_limits<double>::epsilon() &&
1268  alphaRho < std::numeric_limits<double>::epsilon())
1270 
1271  const int nb_integration_pts = row_data.getN().size1();
1272  const int row_nb_dofs = row_data.getIndices().size();
1273  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1275  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
1276 
1277  );
1278  };
1280 
1281  auto v = getVolume();
1282  auto t_w = getFTensor0IntegrationWeight();
1283 
1284  int row_nb_base_functions = row_data.getN().size2();
1285  auto t_row_base_fun = row_data.getFTensor0N();
1286 
1287  double ts_scale = alphaW * getTSa();
1288  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
1289  ts_scale += alphaRho * getTSaa();
1290 
1291  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1292  double a = v * t_w * ts_scale;
1293 
1294  int rr = 0;
1295  for (; rr != row_nb_dofs / 3; ++rr) {
1296 
1297  auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
1298  auto t_m = get_ftensor1(K, 3 * rr, 0);
1299  for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1300  const double b = a * t_row_base_fun * t_col_base_fun;
1301  t_m(i) -= b;
1302  ++t_m;
1303  ++t_col_base_fun;
1304  }
1305 
1306  ++t_row_base_fun;
1307  }
1308 
1309  for (; rr != row_nb_base_functions; ++rr)
1310  ++t_row_base_fun;
1311 
1312  ++t_w;
1313  }
1314 
1316 }
1317 
1318 MoFEMErrorCode OpSpatialPhysical_du_dP::integrate(EntData &row_data,
1319  EntData &col_data) {
1321 
1323 
1324  int nb_integration_pts = row_data.getN().size1();
1325  int row_nb_dofs = row_data.getIndices().size();
1326  int col_nb_dofs = col_data.getIndices().size();
1327  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1329 
1330  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1331 
1332  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1333 
1334  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1335 
1336  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1337 
1338  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1339 
1340  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
1341  };
1342 
1345 
1346  auto v = getVolume();
1347  auto t_w = getFTensor0IntegrationWeight();
1348  auto t_approx_P_adjont_log_du_dP =
1349  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1350  int row_nb_base_functions = row_data.getN().size2();
1351  auto t_row_base_fun = row_data.getFTensor0N();
1352  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1353  double a = v * t_w;
1354  int rr = 0;
1355  for (; rr != row_nb_dofs / 6; ++rr) {
1356 
1357  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1358  auto t_m = get_ftensor3(K, 6 * rr, 0);
1359 
1360  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1361  t_m(L, i) +=
1362  a * (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(j)) *
1363  t_row_base_fun;
1364  ++t_col_base_fun;
1365  ++t_m;
1366  }
1367 
1368  ++t_row_base_fun;
1369  }
1370  for (; rr != row_nb_base_functions; ++rr)
1371  ++t_row_base_fun;
1372  ++t_w;
1373  ++t_approx_P_adjont_log_du_dP;
1374  }
1376 }
1377 
1378 MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
1379  EntData &col_data) {
1381 
1383 
1384  int nb_integration_pts = row_data.getN().size1();
1385  int row_nb_dofs = row_data.getIndices().size();
1386  int col_nb_dofs = col_data.getIndices().size();
1387  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1389  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
1390  &m(r + 5, c));
1391  };
1392 
1395 
1396  auto v = getVolume();
1397  auto t_w = getFTensor0IntegrationWeight();
1398  auto t_approx_P_adjont_log_du_dP =
1399  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1400 
1401  int row_nb_base_functions = row_data.getN().size2();
1402  auto t_row_base_fun = row_data.getFTensor0N();
1403  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1404  double a = v * t_w;
1405  int rr = 0;
1406  for (; rr != row_nb_dofs / 6; ++rr) {
1407  auto t_m = get_ftensor2(K, 6 * rr, 0);
1408  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1409  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1410  t_m(L) +=
1411  a * (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(i, j)) *
1412  t_row_base_fun;
1413  ++t_m;
1414  ++t_col_base_fun;
1415  }
1416  ++t_row_base_fun;
1417  }
1418  for (; rr != row_nb_base_functions; ++rr)
1419  ++t_row_base_fun;
1420  ++t_w;
1421  ++t_approx_P_adjont_log_du_dP;
1422  }
1424 }
1425 
1426 MoFEMErrorCode OpSpatialPhysical_du_du::integrate(EntData &row_data,
1427  EntData &col_data) {
1429  if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
1430  CHKERR integratePiola(row_data, col_data);
1431  } else {
1432  if (polyConvex) {
1433  CHKERR integratePolyconvexHencky(row_data, col_data);
1434  } else {
1435  CHKERR integrateHencky(row_data, col_data);
1436  }
1437  }
1439 }
1440 
1441 MoFEMErrorCode OpSpatialPhysical_du_du::integrateHencky(EntData &row_data,
1442  EntData &col_data) {
1444 
1447  auto t_L = symm_L_tensor();
1448  auto t_diff = diff_tensor();
1449 
1450  int nb_integration_pts = row_data.getN().size1();
1451  int row_nb_dofs = row_data.getIndices().size();
1452  int col_nb_dofs = col_data.getIndices().size();
1453 
1454  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1456  size_symm>(
1457 
1458  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1459  &m(r + 0, c + 4), &m(r + 0, c + 5),
1460 
1461  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1462  &m(r + 1, c + 4), &m(r + 1, c + 5),
1463 
1464  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1465  &m(r + 2, c + 4), &m(r + 2, c + 5),
1466 
1467  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1468  &m(r + 3, c + 4), &m(r + 3, c + 5),
1469 
1470  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1471  &m(r + 4, c + 4), &m(r + 4, c + 5),
1472 
1473  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1474  &m(r + 5, c + 4), &m(r + 5, c + 5)
1475 
1476  );
1477  };
1484 
1485  auto v = getVolume();
1486  auto t_w = getFTensor0IntegrationWeight();
1487 
1488  auto t_approx_P_adjont_dstretch =
1489  getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1490  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1491  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1492  auto &nbUniq = dataAtPts->nbUniq;
1493 
1494  int row_nb_base_functions = row_data.getN().size2();
1495  auto t_row_base_fun = row_data.getFTensor0N();
1496 
1497  auto get_dP = [&]() {
1498  dP.resize(size_symm * size_symm, nb_integration_pts, false);
1499  auto ts_a = getTSa();
1500 
1501  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
1503  t_dP_tmp(L, J) = -(1 + alphaU * ts_a) *
1504  (t_L(i, j, L) *
1505  ((t_D(i, j, m, n) * t_diff(m, n, k, l)) * t_L(k, l, J)));
1506 
1507 
1508  if (EshelbianCore::stretchSelector > LINEAR) {
1509  auto t_approx_P_adjont_dstretch =
1510  getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1511  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1512  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1513  auto &nbUniq = dataAtPts->nbUniq;
1514 
1515  auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1516  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1517 
1518  // Work of symmetric tensor on undefined tensor is equal to the work
1519  // of the symmetric part of it
1521  t_sym(i, j) = (t_approx_P_adjont_dstretch(i, j) ||
1522  t_approx_P_adjont_dstretch(j, i));
1523  t_sym(i, j) /= 2.0;
1524  auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1525  t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
1526  EshelbianCore::dd_f, t_sym, nbUniq[gg]);
1527  t_dP(L, J) = t_L(i, j, L) *
1528  ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
1529  t_L(k, l, J)) /
1530  2. +
1531  t_dP_tmp(L, J);
1532 
1533  ++t_dP;
1534  ++t_approx_P_adjont_dstretch;
1535  ++t_eigen_vals;
1536  ++t_eigen_vecs;
1537  }
1538  } else {
1539  auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1540  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1541  t_dP(L, J) = t_dP_tmp(L, J);
1542  ++t_dP;
1543  }
1544  }
1545 
1546  return getFTensor2FromMat<size_symm, size_symm>(dP);
1547  };
1548 
1549  auto t_dP = get_dP();
1550  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1551  double a = v * t_w;
1552 
1553  int rr = 0;
1554  for (; rr != row_nb_dofs / 6; ++rr) {
1555  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1556  auto t_m = get_ftensor2(K, 6 * rr, 0);
1557  for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1558  const double b = a * t_row_base_fun * t_col_base_fun;
1559  t_m(L, J) += b * t_dP(L, J);
1560  ++t_m;
1561  ++t_col_base_fun;
1562  }
1563  ++t_row_base_fun;
1564  }
1565 
1566  for (; rr != row_nb_base_functions; ++rr) {
1567  ++t_row_base_fun;
1568  }
1569 
1570  ++t_w;
1571  ++t_dP;
1572  }
1574 }
1575 
1577 OpSpatialPhysical_du_du::integratePolyconvexHencky(EntData &row_data,
1578  EntData &col_data) {
1580 
1583  auto t_L = symm_L_tensor();
1584  auto t_diff = diff_tensor();
1585 
1586  int nb_integration_pts = row_data.getN().size1();
1587  int row_nb_dofs = row_data.getIndices().size();
1588  int col_nb_dofs = col_data.getIndices().size();
1589 
1590  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1592  size_symm>(
1593 
1594  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1595  &m(r + 0, c + 4), &m(r + 0, c + 5),
1596 
1597  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1598  &m(r + 1, c + 4), &m(r + 1, c + 5),
1599 
1600  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1601  &m(r + 2, c + 4), &m(r + 2, c + 5),
1602 
1603  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1604  &m(r + 3, c + 4), &m(r + 3, c + 5),
1605 
1606  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1607  &m(r + 4, c + 4), &m(r + 4, c + 5),
1608 
1609  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1610  &m(r + 5, c + 4), &m(r + 5, c + 5)
1611 
1612  );
1613  };
1620 
1621  auto v = getVolume();
1622  auto t_w = getFTensor0IntegrationWeight();
1623 
1624  int row_nb_base_functions = row_data.getN().size2();
1625  auto t_row_base_fun = row_data.getFTensor0N();
1626 
1627  auto get_dP = [&]() {
1628  dP.resize(size_symm * size_symm, nb_integration_pts, false);
1629  auto ts_a = getTSa();
1630 
1631  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
1632 
1633  constexpr double nohat_k = 1. / 4;
1634  constexpr double hat_k = 1. / 8;
1635  double mu = dataAtPts->mu;
1636  double lambda = dataAtPts->lambda;
1637 
1638  constexpr double third = boost::math::constants::third<double>();
1640  auto t_diff_deviator = diff_deviator(diff_tensor());
1641 
1642  auto t_approx_P_adjont_dstretch =
1643  getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1644  auto t_log_streach_h1 =
1645  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
1646  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1647  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1648  auto &nbUniq = dataAtPts->nbUniq;
1649 
1650  auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1651  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1652 
1653  double log_det = t_log_streach_h1(i, i);
1654  double log_det2 = log_det * log_det;
1656  t_dev(i, j) = t_log_streach_h1(i, j) - t_kd(i, j) * (third * log_det);
1657  double dev_norm2 = t_dev(i, j) * t_dev(i, j);
1658 
1659  auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
1660  auto B = lambda * std::exp(hat_k * log_det2) * log_det;
1661 
1662  FTensor::Tensor2_symmetric<double, 3> t_A_diff, t_B_diff;
1663  t_A_diff(i, j) =
1664  (A * 2 * nohat_k) * (t_dev(k, l) * t_diff_deviator(k, l, i, j));
1665  t_B_diff(i, j) = (B * 2 * hat_k) * log_det * t_kd(i, j) +
1666  lambda * std::exp(hat_k * log_det2) * t_kd(i, j);
1668  t_dT(i, j, k, l) =
1669  t_A_diff(i, j) * (t_dev(m, n) * t_diff_deviator(m, n, k, l))
1670 
1671  +
1672 
1673  A * t_diff_deviator(m, n, i, j) * t_diff_deviator(m, n, k, l)
1674 
1675  +
1676 
1677  t_B_diff(i, j) * t_kd(k, l);
1678 
1679  t_dP(L, J) = -t_L(i, j, L) *
1680  ((
1681 
1682  t_dT(i, j, k, l)
1683 
1684  +
1685 
1686  (alphaU * ts_a) * (t_D(i, j, m, n) * t_diff(m, n, k, l)
1687 
1688  )) *
1689  t_L(k, l, J));
1690 
1691  // Work of symmetric tensor on undefined tensor is equal to the work
1692  // of the symmetric part of it
1693  if (EshelbianCore::stretchSelector > LINEAR) {
1695  t_sym(i, j) = (t_approx_P_adjont_dstretch(i, j) ||
1696  t_approx_P_adjont_dstretch(j, i));
1697  t_sym(i, j) /= 2.0;
1698  auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1699  t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
1700  EshelbianCore::dd_f, t_sym, nbUniq[gg]);
1701  t_dP(L, J) += t_L(i, j, L) *
1702  ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
1703  t_L(k, l, J)) /
1704  2.;
1705  }
1706 
1707  ++t_dP;
1708  ++t_approx_P_adjont_dstretch;
1709  ++t_log_streach_h1;
1710  ++t_eigen_vals;
1711  ++t_eigen_vecs;
1712  }
1713 
1714  return getFTensor2FromMat<size_symm, size_symm>(dP);
1715  };
1716 
1717  auto t_dP = get_dP();
1718  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1719  double a = v * t_w;
1720 
1721  int rr = 0;
1722  for (; rr != row_nb_dofs / 6; ++rr) {
1723  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1724  auto t_m = get_ftensor2(K, 6 * rr, 0);
1725  for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1726  const double b = a * t_row_base_fun * t_col_base_fun;
1727  t_m(L, J) += b * t_dP(L, J);
1728  ++t_m;
1729  ++t_col_base_fun;
1730  }
1731  ++t_row_base_fun;
1732  }
1733 
1734  for (; rr != row_nb_base_functions; ++rr) {
1735  ++t_row_base_fun;
1736  }
1737 
1738  ++t_w;
1739  ++t_dP;
1740  }
1742 }
1743 
1744 MoFEMErrorCode OpSpatialPhysical_du_du::integratePiola(EntData &row_data,
1745  EntData &col_data) {
1747 
1750  auto t_L = symm_L_tensor();
1751  auto t_diff = diff_tensor();
1752 
1753  int nb_integration_pts = row_data.getN().size1();
1754  int row_nb_dofs = row_data.getIndices().size();
1755  int col_nb_dofs = col_data.getIndices().size();
1756 
1757  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1759  size_symm>(
1760 
1761  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1762  &m(r + 0, c + 4), &m(r + 0, c + 5),
1763 
1764  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1765  &m(r + 1, c + 4), &m(r + 1, c + 5),
1766 
1767  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1768  &m(r + 2, c + 4), &m(r + 2, c + 5),
1769 
1770  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1771  &m(r + 3, c + 4), &m(r + 3, c + 5),
1772 
1773  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1774  &m(r + 4, c + 4), &m(r + 4, c + 5),
1775 
1776  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1777  &m(r + 5, c + 4), &m(r + 5, c + 5)
1778 
1779  );
1780  };
1787 
1788  auto v = getVolume();
1789  auto t_w = getFTensor0IntegrationWeight();
1790 
1791  auto get_dP = [&]() {
1792  dP.resize(size_symm * size_symm, nb_integration_pts, false);
1793  auto ts_a = getTSa();
1794 
1795  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
1796  auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
1797  auto t_approx_P_adjont_dstretch =
1798  getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1799  auto t_dot_log_u =
1800  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
1801  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1802  auto t_diff_u =
1803  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
1804  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1805  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1806  auto &nbUniq = dataAtPts->nbUniq;
1807 
1808  auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1809  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1810 
1812  t_deltaP(i, j) = t_approx_P_adjont_dstretch(i, j) - t_P(i, j);
1813 
1815  t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
1816  t_dP(L, J) =
1817  -t_Ldiff_u(i, j, L) * (r_P_du(i, j, k, l) * t_Ldiff_u(k, l, J));
1818  // viscous stress derivative
1819  t_dP(L, J) -= (alphaU * ts_a) *
1820  (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
1821 
1822  if (EshelbianCore::stretchSelector > LINEAR) {
1824  t_deltaP_sym(i, j) = (t_deltaP(i, j) || t_deltaP(j, i));
1825  t_deltaP_sym(i, j) /= 2.0;
1826  auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1827  t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
1828  EshelbianCore::dd_f, t_deltaP_sym, nbUniq[gg]);
1829  t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP2(i, j, k, l) * t_L(k, l, J));
1830  }
1831 
1832  ++t_P;
1833  ++t_dP;
1834  ++r_P_du;
1835  ++t_approx_P_adjont_dstretch;
1836  ++t_dot_log_u;
1837  ++t_u;
1838  ++t_diff_u;
1839  ++t_eigen_vals;
1840  ++t_eigen_vecs;
1841  }
1842 
1843  return getFTensor2FromMat<size_symm, size_symm>(dP);
1844  };
1845 
1846  int row_nb_base_functions = row_data.getN().size2();
1847  auto t_row_base_fun = row_data.getFTensor0N();
1848 
1849  auto t_dP = get_dP();
1850  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1851  double a = v * t_w;
1852 
1853  int rr = 0;
1854  for (; rr != row_nb_dofs / 6; ++rr) {
1855  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1856  auto t_m = get_ftensor2(K, 6 * rr, 0);
1857  for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1858  const double b = a * t_row_base_fun * t_col_base_fun;
1859  t_m(L, J) += b * t_dP(L, J);
1860  ++t_m;
1861  ++t_col_base_fun;
1862  }
1863  ++t_row_base_fun;
1864  }
1865 
1866  for (; rr != row_nb_base_functions; ++rr) {
1867  ++t_row_base_fun;
1868  }
1869 
1870  ++t_w;
1871  ++t_dP;
1872  }
1874 }
1875 
1876 MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
1877  EntData &col_data) {
1879 
1881  auto t_L = symm_L_tensor();
1882 
1883  int row_nb_dofs = row_data.getIndices().size();
1884  int col_nb_dofs = col_data.getIndices().size();
1885  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1887 
1888  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1889 
1890  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1891 
1892  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1893 
1894  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1895 
1896  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1897 
1898  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1899 
1900  );
1901  };
1907 
1908  auto v = getVolume();
1909  auto t_w = getFTensor0IntegrationWeight();
1910  auto t_approx_P_adjont_log_du_domega =
1911  getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
1912 
1913  int row_nb_base_functions = row_data.getN().size2();
1914  auto t_row_base_fun = row_data.getFTensor0N();
1915 
1916  int nb_integration_pts = row_data.getN().size1();
1917 
1918  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1919  double a = v * t_w;
1920 
1921  int rr = 0;
1922  for (; rr != row_nb_dofs / 6; ++rr) {
1923  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1924  auto t_m = get_ftensor3(K, 6 * rr, 0);
1925  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1926  double v = a * t_row_base_fun * t_col_base_fun;
1927  t_m(L, k) += v * t_approx_P_adjont_log_du_domega(k, L);
1928  ++t_m;
1929  ++t_col_base_fun;
1930  }
1931  ++t_row_base_fun;
1932  }
1933 
1934  for (; rr != row_nb_base_functions; ++rr)
1935  ++t_row_base_fun;
1936 
1937  ++t_w;
1938  ++t_approx_P_adjont_log_du_domega;
1939  }
1940 
1942 }
1943 
1944 MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
1945  EntData &col_data) {
1947  int nb_integration_pts = getGaussPts().size2();
1948  int row_nb_dofs = row_data.getIndices().size();
1949  int col_nb_dofs = col_data.getIndices().size();
1950  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1952 
1953  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1954 
1955  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1956 
1957  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1958 
1959  );
1960  };
1967  auto v = getVolume();
1968  auto t_w = getFTensor0IntegrationWeight();
1969  auto t_levi_kirchoff_dP =
1970  getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
1971 
1972  int row_nb_base_functions = row_data.getN().size2();
1973  auto t_row_base_fun = row_data.getFTensor0N();
1974 
1975  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1976  double a = v * t_w;
1977  int rr = 0;
1978  for (; rr != row_nb_dofs / 3; ++rr) {
1979  double b = a * t_row_base_fun;
1980  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1981  auto t_m = get_ftensor2(K, 3 * rr, 0);
1982  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1983  t_m(k, i) += b * (t_levi_kirchoff_dP(i, l, k) * t_col_base_fun(l));
1984  ++t_m;
1985  ++t_col_base_fun;
1986  }
1987  ++t_row_base_fun;
1988  }
1989  for (; rr != row_nb_base_functions; ++rr) {
1990  ++t_row_base_fun;
1991  }
1992 
1993  ++t_w;
1994  ++t_levi_kirchoff_dP;
1995  }
1997 }
1998 
1999 MoFEMErrorCode OpSpatialRotation_domega_dBubble::integrate(EntData &row_data,
2000  EntData &col_data) {
2002  int nb_integration_pts = getGaussPts().size2();
2003  int row_nb_dofs = row_data.getIndices().size();
2004  int col_nb_dofs = col_data.getIndices().size();
2005 
2006  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
2008  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
2009  };
2010 
2016 
2017  auto v = getVolume();
2018  auto t_w = getFTensor0IntegrationWeight();
2019  auto t_levi_kirchoff_dP =
2020  getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2021 
2022  int row_nb_base_functions = row_data.getN().size2();
2023  auto t_row_base_fun = row_data.getFTensor0N();
2024 
2025  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2026  double a = v * t_w;
2027  int rr = 0;
2028  for (; rr != row_nb_dofs / 3; ++rr) {
2029  double b = a * t_row_base_fun;
2030  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
2031  auto t_m = get_ftensor1(K, 3 * rr, 0);
2032  for (int cc = 0; cc != col_nb_dofs; ++cc) {
2033  t_m(k) += b * (t_levi_kirchoff_dP(i, j, k) * t_col_base_fun(i, j));
2034  ++t_m;
2035  ++t_col_base_fun;
2036  }
2037  ++t_row_base_fun;
2038  }
2039 
2040  for (; rr != row_nb_base_functions; ++rr) {
2041  ++t_row_base_fun;
2042  }
2043  ++t_w;
2044  ++t_levi_kirchoff_dP;
2045  }
2047 }
2048 
2049 MoFEMErrorCode OpSpatialRotation_domega_domega::integrate(EntData &row_data,
2050  EntData &col_data) {
2052  int nb_integration_pts = getGaussPts().size2();
2053  int row_nb_dofs = row_data.getIndices().size();
2054  int col_nb_dofs = col_data.getIndices().size();
2055  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2057 
2058  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2059 
2060  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2061 
2062  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
2063 
2064  );
2065  };
2072 
2073  auto v = getVolume();
2074  auto t_w = getFTensor0IntegrationWeight();
2075  auto t_levi_kirchoff_domega =
2076  getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
2077  int row_nb_base_functions = row_data.getN().size2();
2078  auto t_row_base_fun = row_data.getFTensor0N();
2079  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2080  double a = v * t_w;
2081  int rr = 0;
2082  for (; rr != row_nb_dofs / 3; ++rr) {
2083  auto t_m = get_ftensor2(K, 3 * rr, 0);
2084  const double b = a * t_row_base_fun;
2085  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2086  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2087  t_m(k, l) += (b * t_col_base_fun) * t_levi_kirchoff_domega(k, l);
2088  ++t_m;
2089  ++t_col_base_fun;
2090  }
2091  ++t_row_base_fun;
2092  }
2093  for (; rr != row_nb_base_functions; ++rr) {
2094  ++t_row_base_fun;
2095  }
2096  ++t_w;
2097  ++t_levi_kirchoff_domega;
2098  }
2100 }
2101 
2102 MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
2103  EntData &col_data) {
2105 
2106  int nb_integration_pts = row_data.getN().size1();
2107  int row_nb_dofs = row_data.getIndices().size();
2108  int col_nb_dofs = col_data.getIndices().size();
2109 
2110  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2112 
2113  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
2114 
2115  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
2116 
2117  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
2118 
2119  );
2120  };
2121 
2128 
2129  auto v = getVolume();
2130  auto t_w = getFTensor0IntegrationWeight();
2131  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2132  int row_nb_base_functions = row_data.getN().size2() / 3;
2133  auto t_row_base_fun = row_data.getFTensor1N<3>();
2134 
2135  const double ts_a = getTSa();
2136 
2137  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2138  double a = v * t_w;
2139 
2140  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
2141 
2142  int rr = 0;
2143  for (; rr != row_nb_dofs / 3; ++rr) {
2144 
2146  t_PRT(i, k) = t_row_base_fun(j) * t_h_domega(i, j, k);
2147 
2148  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2149  auto t_m = get_ftensor2(K, 3 * rr, 0);
2150  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2151  t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
2152  ++t_m;
2153  ++t_col_base_fun;
2154  }
2155 
2156  ++t_row_base_fun;
2157  }
2158 
2159  for (; rr != row_nb_base_functions; ++rr)
2160  ++t_row_base_fun;
2161  ++t_w;
2162  ++t_h_domega;
2163  }
2165 }
2166 
2168 OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
2169  EntData &col_data) {
2171 
2172  int nb_integration_pts = row_data.getN().size1();
2173  int row_nb_dofs = row_data.getIndices().size();
2174  int col_nb_dofs = col_data.getIndices().size();
2175 
2176  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
2178  &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
2179  };
2180 
2187 
2188  auto v = getVolume();
2189  auto t_w = getFTensor0IntegrationWeight();
2190  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2191  int row_nb_base_functions = row_data.getN().size2() / 9;
2192  auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
2193 
2194  const double ts_a = getTSa();
2195 
2196  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2197  double a = v * t_w;
2198 
2199  int rr = 0;
2200  for (; rr != row_nb_dofs; ++rr) {
2201 
2203  t_PRT(k) = t_row_base_fun(i, j) * t_h_domega(i, j, k);
2204 
2205  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2206  auto t_m = get_ftensor2(K, rr, 0);
2207  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2208  t_m(j) += (a * t_col_base_fun) * t_PRT(j);
2209  ++t_m;
2210  ++t_col_base_fun;
2211  }
2212 
2213  ++t_row_base_fun;
2214  }
2215 
2216  for (; rr != row_nb_base_functions; ++rr)
2217  ++t_row_base_fun;
2218 
2219  ++t_w;
2220  ++t_h_domega;
2221  }
2223 }
2224 
2225 MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
2226  EntData &data) {
2228 
2229  auto create_tag = [this](const std::string tag_name, const int size) {
2230  double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2231  Tag th;
2232  CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
2233  th, MB_TAG_CREAT | MB_TAG_SPARSE,
2234  def_VAL);
2235  return th;
2236  };
2237 
2238  Tag th_w = create_tag("SpatialDisplacement", 3);
2239  Tag th_omega = create_tag("Omega", 3);
2240  Tag th_approxP = create_tag("Piola1Stress", 9);
2241  Tag th_sigma = create_tag("CauchyStress", 9);
2242  Tag th_log_u = create_tag("LogSpatialStretch", 9);
2243  Tag th_u = create_tag("SpatialStretch", 9);
2244  Tag th_h = create_tag("h", 9);
2245  Tag th_X = create_tag("X", 3);
2246  Tag th_detF = create_tag("detF", 1);
2247  Tag th_angular_momentum = create_tag("AngularMomentum", 3);
2248 
2249  Tag th_u_eig_vec = create_tag("SpatialStretchEigenVec", 9);
2250  Tag th_u_eig_vals = create_tag("SpatialStretchEigenVals", 3);
2251  Tag th_traction = create_tag("traction", 3);
2252 
2253  Tag th_disp = create_tag("H1Displacement", 3);
2254  Tag th_disp_error = create_tag("DisplacementError", 1);
2255  Tag th_lambda_disp = create_tag("ContactDisplacement", 3);
2256 
2257  auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
2258  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
2259  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2260  auto t_log_u =
2261  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
2262  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
2263  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
2264  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2265  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
2266  auto t_coords = getFTensor1CoordsAtGaussPts();
2267  auto t_normal = getFTensor1NormalsAtGaussPts();
2268  auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
2269  auto t_lambda_disp = getFTensor1FromMat<3>(dataAtPts->contactL2AtPts);
2270 
2275 
2276  auto set_float_precision = [](const double x) {
2277  if (std::abs(x) < std::numeric_limits<float>::epsilon())
2278  return 0.;
2279  else
2280  return x;
2281  };
2282 
2283  // scalars
2284  auto save_scal_tag = [&](auto &th, auto v, const int gg) {
2286  v = set_float_precision(v);
2287  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, &v);
2289  };
2290 
2291  // vectors
2292  VectorDouble3 v(3);
2293  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
2294  auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
2296  t_v(i) = t_d(i);
2297  for (auto &a : v.data())
2298  a = set_float_precision(a);
2299  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2300  &*v.data().begin());
2302  };
2303 
2304  // tensors
2305 
2306  MatrixDouble3by3 m(3, 3);
2308  &m(0, 0), &m(0, 1), &m(0, 2),
2309 
2310  &m(1, 0), &m(1, 1), &m(1, 2),
2311 
2312  &m(2, 0), &m(2, 1), &m(2, 2));
2313 
2314  auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
2316  t_m(i, j) = t_d(i, j);
2317  for (auto &v : m.data())
2318  v = set_float_precision(v);
2319  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2320  &*m.data().begin());
2322  };
2323 
2324  const auto nb_gauss_pts = getGaussPts().size2();
2325  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
2326 
2327  FTensor::Tensor1<double, 3> t_traction;
2328  t_traction(i) = t_approx_P(i, j) * t_normal(j) / t_normal.l2();
2329 
2330  // vectors
2331  CHKERR save_vec_tag(th_w, t_w, gg);
2332  CHKERR save_vec_tag(th_X, t_coords, gg);
2333  CHKERR save_vec_tag(th_omega, t_omega, gg);
2334  CHKERR save_vec_tag(th_traction, t_traction, gg);
2335 
2336  // tensors
2337  CHKERR save_mat_tag(th_h, t_h, gg);
2338 
2339  FTensor::Tensor2<double, 3, 3> t_log_u_tmp;
2340  for (int ii = 0; ii != 3; ++ii)
2341  for (int jj = 0; jj != 3; ++jj)
2342  t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
2343 
2344  CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
2345 
2347  for (int ii = 0; ii != 3; ++ii)
2348  for (int jj = 0; jj != 3; ++jj)
2349  t_u_tmp(ii, jj) = t_u(ii, jj);
2350 
2351  CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2352  CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2353  CHKERR save_vec_tag(th_disp, t_disp, gg);
2354  CHKERR save_vec_tag(th_lambda_disp, t_lambda_disp, gg);
2355 
2356  double u_error = sqrt((t_disp(i) - t_w(i)) * (t_disp(i) - t_w(i)));
2357  CHKERR save_scal_tag(th_disp_error, u_error, gg);
2358 
2359  const double jac = determinantTensor3by3(t_h);
2361  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
2362  CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2363  CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
2364 
2366  t_levi(k) = t_levi_kirchoff(k);
2367  CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
2368  &t_levi(0));
2369 
2370  auto get_eiegn_vector_symmetric = [&](auto &t_u) {
2372 
2373  for (int ii = 0; ii != 3; ++ii)
2374  for (int jj = 0; jj != 3; ++jj)
2375  t_m(ii, jj) = t_u(ii, jj);
2376 
2377  VectorDouble3 eigen_values(3);
2378  auto t_eigen_values = getFTensor1FromArray<3>(eigen_values);
2379  CHKERR computeEigenValuesSymmetric(t_m, t_eigen_values);
2380 
2381  CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
2382  &*m.data().begin());
2383  CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
2384  &*eigen_values.data().begin());
2385 
2387  };
2388 
2389  CHKERR get_eiegn_vector_symmetric(t_u);
2390 
2391  ++t_w;
2392  ++t_h;
2393  ++t_log_u;
2394  ++t_u;
2395  ++t_omega;
2396  ++t_R;
2397  ++t_approx_P;
2398  ++t_levi_kirchoff;
2399  ++t_coords;
2400  ++t_normal;
2401  ++t_disp;
2402  ++t_lambda_disp;
2403  }
2404 
2406 }
2407 
2409 OpCalculateStrainEnergy::doWork(int side, EntityType type,
2412  if (type == MBTET) {
2413  int nb_integration_pts = data.getN().size1();
2414  auto v = getVolume();
2415  auto t_w = getFTensor0IntegrationWeight();
2416  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2417  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2418 
2421 
2422  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2423  const double a = t_w * v;
2424  // FIXME: this is wrong, energy should be calculated in material
2425  (*energy) += a * t_P(i, J) * t_h(i, J);
2426  ++t_w;
2427  ++t_P;
2428  ++t_h;
2429  }
2430  }
2432 }
2433 
2434 } // namespace EshelbianPlasticity
EshelbianPlasticity::sort_eigen_vals
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
Definition: EshelbianOperators.cpp:218
EigenMatrix::getMat
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
Definition: MatrixFunction.cpp:53
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
EshelbianPlasticity::is_eq
auto is_eq(const double &a, const double &b)
Definition: EshelbianOperators.cpp:203
FTensor
JSON compatible output.
Definition: Christof_constructor.hpp:6
EshelbianPlasticity::get_rotation_form_vector
auto get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:83
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:21
EshelbianPlasticity::size_symm
constexpr static auto size_symm
Definition: EshelbianOperators.cpp:47
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::ForcesAndSourcesCore::getRuleHook
RuleHookFun getRuleHook
Hook to get rule.
Definition: ForcesAndSourcesCore.hpp:42
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FTensor::Tensor1::l2
T l2() const
Definition: Tensor1_value.hpp:84
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
EshelbianPlasticity::isEq
Definition: EshelbianOperators.cpp:193
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1241
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
EshelbianPlasticity::diff_tensor
auto diff_tensor()
Definition: EshelbianOperators.cpp:49
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
EshelbianPlasticity::isEq::absMax
static double absMax
Definition: EshelbianOperators.cpp:198
EshelbianPlasticity::get_uniq_nb
auto get_uniq_nb(double *ptr)
Definition: EshelbianOperators.cpp:207
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::EntitiesFieldData::EntData::getFTensor2DiffN
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Definition: EntitiesFieldData.cpp:660
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
convert.type
type
Definition: convert.py:64
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
lapack_wrap.h
MatrixFunction.hpp
MoFEM::EntitiesFieldData::EntData::getFTensor2N
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:762
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
EshelbianPlasticity::FaceRule
Definition: EshelbianPlasticity.cpp:748
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1256
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
EshelbianPlasticity::TensorTypeExtractor
Definition: EshelbianOperators.cpp:75
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
EshelbianPlasticity::diff_deviator
auto diff_deviator(FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
Definition: EshelbianOperators.cpp:20
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
EshelbianPlasticity::symm_L_tensor
auto symm_L_tensor()
Definition: EshelbianOperators.cpp:60
FTensor::Index< 'i', 3 >
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
convert.n
n
Definition: convert.py:82
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1518
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:20
FTensor::dd
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MoFEM::TSMethod::ts_t
PetscReal ts_t
time
Definition: LoopMethods.hpp:162
FTensor::Dg
Definition: Dg_value.hpp:9
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:20
FTensor::Tensor0
Definition: Tensor0.hpp:16
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:21
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1305
EshelbianPlasticity::get_diff2_rotation_form_vector
auto get_diff2_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:167
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::EntitiesFieldData::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:640
EshelbianPlasticity::RotSelector
RotSelector
Definition: EshelbianPlasticity.hpp:20
FTensor::Ddg< double, 3, 3 >
lapack_dposv
static __CLPK_integer lapack_dposv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:211
EshelbianPlasticity::OpTractionBc
Definition: EshelbianPlasticity.hpp:612
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
EigenMatrix::getDiffDiffMat
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
Definition: MatrixFunction.cpp:78
EshelbianPlasticity::get_diff_rotation_form_vector
auto get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:119
MoFEM::Types::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
third
constexpr double third
Definition: EshelbianADOL-C.cpp:14
EshelbianPlasticity::isEq::check
static auto check(const double &a, const double &b)
Definition: EshelbianOperators.cpp:194
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:20
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EshelbianPlasticity::TensorTypeExtractor::Type
std::remove_pointer< T >::type Type
Definition: EshelbianOperators.cpp:76
EigenMatrix::getDiffMat
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
Definition: MatrixFunction.cpp:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
EshelbianPlasticity::TensorTypeExtractor< FTensor::PackPtr< T, I > >::Type
std::remove_pointer< T >::type Type
Definition: EshelbianOperators.cpp:79
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1430