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
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
282  dataAtPts->adjointPdUdPAtPts.resize(9 * size_symm, nb_integration_pts, false);
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);
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;
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
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);
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
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);
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
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);
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);
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);
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
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;
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;
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();
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;
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();
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;
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();
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;
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();
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;
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();
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;
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
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) {
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) ||
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;
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
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) ||
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;
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);
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;
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();
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;
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
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::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
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
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 >
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
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
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
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