v0.10.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 <BasicFiniteElements.hpp>
12 
13 #include <EshelbianPlasticity.hpp>
14 #include <boost/math/constants/constants.hpp>
15 
16 #include <lapack_wrap.h>
17 
18 namespace EshelbianPlasticity {
19 
20 template <class T> struct TensorTypeExtractor {
22 };
23 template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
25 };
26 
27 template <typename T>
30  &t_omega) {
35  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
36  t_R(i, j) = t_kd(i, j);
37 
38  const typename TensorTypeExtractor<T>::Type angle =
39  sqrt(t_omega(i) * t_omega(i));
40 
41  constexpr double tol = 1e-18;
42  if (std::abs(angle) < tol) {
43  return t_R;
44  }
45 
47  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
48  const typename TensorTypeExtractor<T>::Type a = sin(angle) / angle;
49  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
50  const typename TensorTypeExtractor<T>::Type b =
51  2. * ss_2 * ss_2 / (angle * angle);
52  t_R(i, j) += a * t_Omega(i, j);
53  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
54 
55  return t_R;
56 }
57 
58 template <typename T>
61  &t_omega) {
62 
67 
68  const typename TensorTypeExtractor<T>::Type angle =
69  sqrt(t_omega(i) * t_omega(i));
70 
71  constexpr double tol = 1e-18;
72  if (std::abs(angle) < tol) {
74  auto t_R = get_rotation_form_vector(t_omega);
75  t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, l, k) * t_R(l, j);
76  return t_diff_R;
77  }
78 
81  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
82 
83  const typename TensorTypeExtractor<T>::Type angle2 = angle * angle;
84  const typename TensorTypeExtractor<T>::Type ss = sin(angle);
85  const typename TensorTypeExtractor<T>::Type a = ss / angle;
86  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
87  const typename TensorTypeExtractor<T>::Type b = 2. * ss_2 * ss_2 / angle2;
88  t_diff_R(i, j, k) = 0;
89  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
90  t_diff_R(i, j, k) +=
91  b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
92  FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
93 
94  const typename TensorTypeExtractor<T>::Type cc = cos(angle);
95  const typename TensorTypeExtractor<T>::Type cc_2 = cos(angle / 2.);
96  const typename TensorTypeExtractor<T>::Type diff_a =
97  (angle * cc - ss) / angle2;
98 
99  const typename TensorTypeExtractor<T>::Type diff_b =
100  (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
101  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
102  t_diff_R(i, j, k) +=
103  diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
104 
105  return t_diff_R;
106 }
107 
108 template <typename T>
109 static MoFEMErrorCode
112  FTensor::Tensor4<T, 3, 3, 3, 3> &t_diff_s_u) {
114 
115  typedef typename TensorTypeExtractor<T>::Type A;
116 
117  auto calc_exponent = [](FTensor::Tensor2<A, 3, 3> &t_log_u,
121 
122  constexpr int N = 800;
123  constexpr double eps = 1e-12;
124 
130 
135 
137  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
138  t_one4(i, j, k, l) = t_kd(j, l) * t_kd(i, k);
139 
140  t_p(i, j) = t_log_u(i, j);
141  t_diff_p(i, j, k, l) = t_one4(i, j, k, l);
142 
143  t_u(i, j) = t_p(i, j);
144  for (auto ii : {0, 1, 2})
145  t_u(ii, ii) += 1;
146  t_diff_u(i, j, k, l) = t_diff_p(i, j, k, l);
147 
148  for (int nn = 2; nn <= N; ++nn) {
149 
150  t_diff_tmp(i, j, k, l) =
151  (t_p(i, m) / static_cast<double>(nn)) * t_one4(m, j, k, l) +
152  t_diff_p(i, m, k, l) * (t_log_u(m, j) / static_cast<double>(nn));
153 
154  t_tmp(i, j) = t_p(i, k) * (t_log_u(k, j) / static_cast<double>(nn));
155 
156  t_p(i, j) = t_tmp(i, j);
157  t_diff_p(i, j, k, l) = t_diff_tmp(i, j, k, l);
158 
159  t_u(i, j) += t_p(i, j);
160  t_diff_u(i, j, k, l) += t_diff_p(i, j, k, l);
161 
162  double e = std::abs(sqrt(t_p(i, j) * t_p(i, j)));
163  if (e < eps)
165  }
167  };
168 
171 
173  for (auto ii : {0, 1, 2})
174  for (auto jj : {0, 1, 2})
175  t_log_u(ii, jj) = t_log_s_u(ii, jj);
176 
177  CHKERR calc_exponent(t_log_u, t_u, t_diff_u);
178 
179  for (auto ii : {0, 1, 2})
180  for (auto jj = 0; jj != 3; ++jj) {
181  t_s_u(ii, jj) = (t_u(ii, jj) + t_u(jj, ii)) / 2.;
182  for (auto kk : {0, 1, 2})
183  for (auto ll = 0; ll <= kk; ++ll) {
184 
185  if (ll != kk)
186  t_diff_s_u(ii, jj, kk, ll) =
187 
188  t_diff_u(ii, jj, kk, ll) + t_diff_u(ii, jj, ll, kk);
189 
190  else
191  t_diff_s_u(ii, jj, kk, ll) = t_diff_u(ii, jj, kk, ll);
192  }
193  }
194 
196 }
197 
198 static MoFEMErrorCode
200  std::array<MatrixDouble, 6> &diff2_u) {
206 
207  for (auto ll : {0, 1, 2, 3, 4, 5})
208  diff2_u[ll].resize(81, log_u.size2(), false);
209 
210  std::array<FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>, 6>
211  t_diff2_u = {
212 
214  &diff2_u[0](0, 0), diff2_u[0].size2()),
216  &diff2_u[1](0, 0), diff2_u[1].size2()),
218  &diff2_u[2](0, 0), diff2_u[2].size2()),
220  &diff2_u[3](0, 0), diff2_u[3].size2()),
222  &diff2_u[4](0, 0), diff2_u[4].size2()),
224  &diff2_u[5](0, 0), diff2_u[5].size2())
225 
226  };
227 
228  constexpr double h = 1e-18;
229  const std::complex<double> ih = h * 1i;
230 
231  int ss = 0;
232  for (auto ii : {0, 1, 2})
233  for (auto jj = 0; jj <= ii; ++jj, ++ss) {
234 
235  auto t_log_u = getFTensor2SymmetricFromMat<3>(log_u);
236  auto t_diff_u =
238  &diff_u(0, 0), diff_u.size2());
239 
241  tc_h(i, j) = 0;
242  tc_h(ii, jj) = ih;
243 
244  for (int gg = 0; gg != log_u.size2(); ++gg) {
245 
248  FTensor::Tensor4<std::complex<double>, 3, 3, 3, 3> tc_diff_u;
249 
250  tc_log_u(i, j) = t_log_u(i, j) + tc_h(i, j);
251  CHKERR tensor_exponent(tc_log_u, tc_u, tc_diff_u);
252 
253  (t_diff2_u[ss])(i, j, k, l) = 0;
254  for (auto kk : {0, 1, 2})
255  for (auto ll : {0, 1, 2})
256  for (auto mm : {0, 1, 2})
257  for (auto nn : {0, 1, 2}) {
258  (t_diff2_u[ss])(kk, ll, mm, nn) =
259  std::imag(tc_diff_u(kk, ll, mm, nn)) / h;
260  }
261 
262  // for (auto kk : {0, 1, 2})
263  // for (auto ll = 0; ll <= 2; ++ll) {
264 
265  // constexpr double eps = 1e-6;
266  // if (std::abs(t_diff_u(kk, ll, ii, jj) -
267  // std::imag(tc_u(kk, ll)) / h) > eps) {
268  // std::ostringstream ss;
269  // ss << "Error " << kk << " " << ll << " : "
270  // << t_diff_u(kk, ll, ii, jj) << " "
271  // << std::imag(tc_u(kk, ll)) / h << " error "
272  // << t_diff_u(kk, ll, ii, jj) - std::imag(tc_u(kk, ll)) / h
273  // << " "
274  // << t_diff_u(kk, ll, ii, jj) / (std::imag(tc_u(kk, ll)) / h)
275  // << endl;
276  // SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "%s",
277  // ss.str().c_str());
278  // }
279  // }
280 
281  ++t_log_u;
282  ++t_diff_u;
283  ++(t_diff2_u[ss]);
284  }
285  }
286 
288 }
289 
290 static MoFEMErrorCode
292  std::array<MatrixDouble, 3> &diff2_omega) {
297 
298  for (auto ll : {0, 1, 2})
299  diff2_omega[ll].resize(27, omega.size2(), false);
300 
301  std::array<FTensor::Tensor3<FTensor::PackPtr<double *, 1>, 3, 3, 3>, 3>
302  t_diff2_omega = {
303 
304  getFTensor3FromMat(diff2_omega[0]),
305 
306  getFTensor3FromMat(diff2_omega[1]),
307 
308  getFTensor3FromMat(diff2_omega[2])
309 
310  };
311 
312  constexpr double h = 1e-18;
313  const std::complex<double> ih = h * 1i;
314 
315  for (auto ii : {0, 1, 2}) {
316 
317  auto t_omega = getFTensor1FromMat<3>(omega);
318 
320  tc_h(i) = 0;
321  tc_h(ii) = ih;
322 
323  for (int gg = 0; gg != omega.size2(); ++gg) {
324 
326  tc_omega(i) = t_omega(i) + tc_h(i);
327 
328  auto tc_diff_R = get_diff_rotation_form_vector(tc_omega);
329 
330  (t_diff2_omega[ii])(i, j, k) = 0;
331  for (auto kk : {0, 1, 2})
332  for (auto ll : {0, 1, 2})
333  for (auto mm : {0, 1, 2})
334  (t_diff2_omega[ii])(kk, ll, mm) =
335  std::imag(tc_diff_R(kk, ll, mm)) / h;
336 
337  ++t_omega;
338  ++(t_diff2_omega[ii]);
339  }
340  }
342 }
343 
344 template <typename T>
347 
351 
352  auto get_dot_fun = [i, j, k](auto &t) {
354  t_c(i, j) = t(i, j);
355  return t_c;
356  };
357 
358  auto get_trace = [](auto &t) {
359  typename TensorTypeExtractor<T>::Type tr = 0;
360  for (auto ii = 0; ii != 3; ++ii)
361  tr += t(ii, ii) / 3.;
362  return tr;
363  };
364 
365  auto get_deviator = [i, j, get_trace](auto &&t) {
366  const auto tr = get_trace(t);
368  t_dev(i, j) = t(i, j);
369  for (auto ii = 0; ii != 3; ++ii)
370  t_dev(ii, ii) -= tr;
371  return t_dev;
372  };
373 
374  return get_deviator(get_dot_fun(t));
375 }
376 
377 template <typename T>
381 
387 
388  auto get_dot_fun = [i, j, k, m, n](auto &t, auto &t_diff) {
390  t_diff_c;
391  t_diff_c(i, j, m, n) = t_diff(i, j, m, n);
392  return t_diff_c;
393  };
394 
395  auto get_trace = [i, j](auto &t_diff) {
397  t_diff_tr(i, j) = 0;
398  for (auto ii = 0; ii != 3; ++ii)
399  for (auto mm = 0; mm != 3; ++mm)
400  for (auto nn = 0; nn != 3; ++nn)
401  t_diff_tr(mm, nn) += t_diff(ii, ii, mm, nn) / 3.;
402  return t_diff_tr;
403  };
404 
405  auto get_deviator = [i, j, m, n, get_trace](auto &&t_diff) {
406  auto t_diff_tr = get_trace(t_diff);
408  t_diff_dev;
409  t_diff_dev(i, j, m, n) = t_diff(i, j, m, n);
410  for (auto mm = 0; mm != 3; ++mm)
411  for (auto nn = 0; nn != 3; ++nn)
412  for (auto ii = 0; ii != 3; ++ii)
413  t_diff_dev(ii, ii, mm, nn) -= t_diff_tr(mm, nn);
414  return t_diff_dev;
415  };
416 
417  return get_deviator(get_dot_fun(t, t_diff));
418 }
419 
420 MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
421  EntityType type,
422  EntData &data) {
424  if (type != MBTET)
426  int nb_integration_pts = data.getN().size1();
431 
432  dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
433  dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
434  dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
435  dataAtPts->streachTensorAtPts.resize(6, nb_integration_pts, false);
436  dataAtPts->diffStreachTensorAtPts.resize(81, nb_integration_pts, false);
437 
438  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
439  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
440  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
441  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
442  auto t_log_u =
443  getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachTensorAtPts);
444  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
445  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
446  &dataAtPts->diffStreachTensorAtPts(0, 0), nb_integration_pts);
447 
448  for (int gg = 0; gg != nb_integration_pts; ++gg) {
449 
450  auto t0_diff = get_diff_rotation_form_vector(t_omega);
451  auto t0 = get_rotation_form_vector(t_omega);
452 
453  t_diff_R(i, j, k) = t0_diff(i, j, k);
454  t_R(i, j) = t0(i, j);
455 
456  CHKERR tensor_exponent(t_log_u, t_u, t_diff_u);
457 
458  t_h(i, j) = t_R(i, k) * t_u(k, j);
459 
460  ++t_h;
461  ++t_R;
462  ++t_diff_R;
463  ++t_log_u;
464  ++t_u;
465  ++t_diff_u;
466  ++t_omega;
467  }
468 
469  CHKERR matrix_exponent_derivative(dataAtPts->logStreachTensorAtPts,
470  dataAtPts->diffStreachTensorAtPts,
471  dataAtPts->expLogUHessian);
472 
474 }
475 
476 MoFEMErrorCode OpSpatialEquilibrium::integrate(EntData &data) {
478  int nb_dofs = data.getIndices().size();
479  int nb_integration_pts = data.getN().size1();
480  auto v = getVolume();
481  auto t_w = getFTensor0IntegrationWeight();
482  auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
483  auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wDotAtPts);
484  int nb_base_functions = data.getN().size2();
485  auto t_row_base_fun = data.getFTensor0N();
486 
488  auto get_ftensor1 = [](auto &v) {
489  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
490  &v[2]);
491  };
492 
493  for (int gg = 0; gg != nb_integration_pts; ++gg) {
494  double a = v * t_w;
495  auto t_nf = get_ftensor1(nF);
496  int bb = 0;
497  for (; bb != nb_dofs / 3; ++bb) {
498  t_nf(i) += a * t_row_base_fun * t_div_P(i);
499  t_nf(i) -= a * t_row_base_fun * alpha_w * t_s_dot_w(i);
500  ++t_nf;
501  ++t_row_base_fun;
502  }
503  for (; bb != nb_base_functions; ++bb)
504  ++t_row_base_fun;
505  ++t_w;
506  ++t_div_P;
507  ++t_s_dot_w;
508  }
509 
511 }
512 
513 MoFEMErrorCode OpSpatialRotation::integrate(EntData &data) {
515  int nb_dofs = data.getIndices().size();
516  int nb_integration_pts = data.getN().size1();
517  auto v = getVolume();
518  auto t_w = getFTensor0IntegrationWeight();
519  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
520  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
521 
522  int nb_base_functions = data.getN().size2();
523  auto t_row_base_fun = data.getFTensor0N();
529  auto get_ftensor1 = [](auto &v) {
530  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
531  &v[2]);
532  };
533 
534  for (int gg = 0; gg != nb_integration_pts; ++gg) {
535  double a = v * t_w;
536  auto t_nf = get_ftensor1(nF);
538  t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
539  FTensor::Tensor1<double, 3> t_leviPRT;
540  t_leviPRT(k) = levi_civita(i, m, k) * t_PRT(i, m);
541  int bb = 0;
542  for (; bb != nb_dofs / 3; ++bb) {
543  t_nf(k) += a * t_row_base_fun * t_leviPRT(k);
544  ++t_nf;
545  ++t_row_base_fun;
546  }
547  for (; bb != nb_base_functions; ++bb)
548  ++t_row_base_fun;
549  ++t_w;
550  ++t_approx_P;
551  ++t_R;
552  }
554 }
555 
556 MoFEMErrorCode OpSpatialPhysical::integrate(EntData &data) {
558  int nb_dofs = data.getIndices().size();
559  int nb_integration_pts = data.getN().size1();
560  auto v = getVolume();
561  auto t_w = getFTensor0IntegrationWeight();
562  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
563  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
564  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
565  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
566  auto t_dot_log_u =
567  getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachDotTensorAtPts);
568  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
569  &dataAtPts->diffStreachTensorAtPts(0, 0),
570  dataAtPts->diffStreachTensorAtPts.size2());
571 
576  auto get_ftensor2 = [](auto &v) {
578  &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
579  };
580 
581  int nb_base_functions = data.getN().size2();
582  auto t_row_base_fun = data.getFTensor0N();
583  for (int gg = 0; gg != nb_integration_pts; ++gg) {
584  double a = v * t_w;
585  auto t_nf = get_ftensor2(nF);
586 
588  t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
590  t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
591 
592  // FTensor::Tensor2<double, 3, 3> t_dot_u;
593  // t_dot_u(i, j) = t_u(i, k) * t_dot_log_u(k, j);
594 
595 
596  // auto t_viscous_stress = calculate_viscous_stress(t_dot_u);
597  // t_viscous_stress(i, j) = t_dot_u(i, j);
598  FTensor::Tensor2_symmetric<double, 3> t_viscous_stress;
599  t_viscous_stress(i, j) = alpha_u * t_dot_log_u(i, j);
600 
602  t1(i, j) = 0;
603  for (int ii : {0, 1, 2})
604  for (int jj = 0; jj <= ii; ++jj)
605  for (int mm : {0, 1, 2})
606  for (int nn : {0, 1, 2})
607  t1(ii, jj) += a * t_diff_u(mm, nn, ii, jj) *
608  (t_deltaP(mm, nn) - t_viscous_stress(mm, nn));
609 
610  int bb = 0;
611  for (; bb != nb_dofs / 6; ++bb) {
612 
613  for (int ii : {0, 1, 2})
614  for (int jj = 0; jj <= ii; ++jj) {
615  t_nf(ii, jj) += t_row_base_fun * t1(ii, jj);
616  }
617 
618  ++t_nf;
619  ++t_row_base_fun;
620  }
621  for (; bb != nb_base_functions; ++bb)
622  ++t_row_base_fun;
623 
624  ++t_w;
625  ++t_approx_P;
626  ++t_P;
627  ++t_R;
628  ++t_dot_log_u;
629  ++t_u;
630  ++t_diff_u;
631  }
633 }
634 
635 MoFEMErrorCode OpSpatialConsistencyP::integrate(EntData &data) {
637  int nb_dofs = data.getIndices().size();
638  int nb_integration_pts = data.getN().size1();
639  auto v = getVolume();
640  auto t_w = getFTensor0IntegrationWeight();
641  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
642  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
643 
644  int nb_base_functions = data.getN().size2() / 3;
645  auto t_row_base_fun = data.getFTensor1N<3>();
650 
651  auto get_ftensor1 = [](auto &v) {
652  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
653  &v[2]);
654  };
655 
656  for (int gg = 0; gg != nb_integration_pts; ++gg) {
657  double a = v * t_w;
658  auto t_nf = get_ftensor1(nF);
659 
660  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
662  t_residuum(i, j) = t_R(i, k) * t_u(k, j) - t_kd(i, j);
663 
664  int bb = 0;
665  for (; bb != nb_dofs / 3; ++bb) {
666  t_nf(i) += a * t_row_base_fun(j) * t_residuum(i, j);
667 
668  ++t_nf;
669  ++t_row_base_fun;
670  }
671 
672  for (; bb != nb_base_functions; ++bb)
673  ++t_row_base_fun;
674 
675  ++t_w;
676  ++t_u;
677  ++t_R;
678  }
679 
681 }
682 
683 MoFEMErrorCode OpSpatialConsistencyBubble::integrate(EntData &data) {
685  int nb_dofs = data.getIndices().size();
686  int nb_integration_pts = data.getN().size1();
687  auto v = getVolume();
688  auto t_w = getFTensor0IntegrationWeight();
689  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
690  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
691 
692  int nb_base_functions = data.getN().size2() / 9;
693  auto t_row_base_fun = data.getFTensor2N<3, 3>();
698 
699  auto get_ftensor0 = [](auto &v) {
701  };
702 
703  for (int gg = 0; gg != nb_integration_pts; ++gg) {
704  double a = v * t_w;
705  auto t_nf = get_ftensor0(nF);
706 
707  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
709  t_residuum(i, j) = t_R(i, k) * t_u(k, j) - t_kd(i, j);
710 
711  int bb = 0;
712  for (; bb != nb_dofs; ++bb) {
713  t_nf += a * t_row_base_fun(i, j) * t_residuum(i, j);
714  ++t_nf;
715  ++t_row_base_fun;
716  }
717  for (; bb != nb_base_functions; ++bb) {
718  ++t_row_base_fun;
719  }
720  ++t_w;
721  ++t_R;
722  ++t_u;
723  }
724 
726 }
727 
728 MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
730  int nb_dofs = data.getIndices().size();
731  int nb_integration_pts = data.getN().size1();
732  auto v = getVolume();
733  auto t_w = getFTensor0IntegrationWeight();
734  auto t_s_w = getFTensor1FromMat<3>(dataAtPts->wAtPts);
735  int nb_base_functions = data.getN().size2() / 3;
736  auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
738  auto get_ftensor1 = [](auto &v) {
739  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
740  &v[2]);
741  };
742 
743  for (int gg = 0; gg != nb_integration_pts; ++gg) {
744  double a = v * t_w;
745  auto t_nf = get_ftensor1(nF);
746  int bb = 0;
747  for (; bb != nb_dofs / 3; ++bb) {
748  double div_row_base = t_row_diff_base_fun(i, i);
749  t_nf(i) += a * div_row_base * t_s_w(i);
750  ++t_nf;
751  ++t_row_diff_base_fun;
752  }
753  for (; bb != nb_base_functions; ++bb) {
754  ++t_row_diff_base_fun;
755  }
756  ++t_w;
757  ++t_s_w;
758  }
759 
761 }
762 
763 MoFEMErrorCode OpDispBc::integrate(EntData &data) {
765  // get entity of face
766  EntityHandle fe_ent = getFEEntityHandle();
767  // interate over all boundary data
768  for (auto &bc : (*bcDispPtr)) {
769  // check if finite element entity is part of boundary condition
770  if (bc.faces.find(fe_ent) != bc.faces.end()) {
771  int nb_dofs = data.getIndices().size();
772  int nb_integration_pts = data.getN().size1();
773  auto t_normal = getFTensor1Normal();
774  auto t_w = getFTensor0IntegrationWeight();
775  int nb_base_functions = data.getN().size2() / 3;
776  auto t_row_base_fun = data.getFTensor1N<3>();
777 
780  auto get_ftensor1 = [](auto &v) {
781  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
782  &v[2]);
783  };
784 
785  // get bc data
786  FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
787  t_bc_disp(i) *= getFEMethod()->ts_t;
788 
789  for (int gg = 0; gg != nb_integration_pts; ++gg) {
791  t_bc_res(i) = t_bc_disp(i);
792  auto t_nf = get_ftensor1(nF);
793  int bb = 0;
794  for (; bb != nb_dofs / 3; ++bb) {
795  t_nf(i) -=
796  t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
797  ++t_nf;
798  ++t_row_base_fun;
799  }
800  for (; bb != nb_base_functions; ++bb)
801  ++t_row_base_fun;
802 
803  ++t_w;
804  }
805  }
806  }
808 }
809 
810 MoFEMErrorCode OpRotationBc::integrate(EntData &data) {
812  // get entity of face
813  EntityHandle fe_ent = getFEEntityHandle();
814  // interate over all boundary data
815  for (auto &bc : (*bcRotPtr)) {
816  // check if finite element entity is part of boundary condition
817  if (bc.faces.find(fe_ent) != bc.faces.end()) {
818  int nb_dofs = data.getIndices().size();
819  int nb_integration_pts = data.getN().size1();
820  auto t_normal = getFTensor1Normal();
821  auto t_w = getFTensor0IntegrationWeight();
822 
823  int nb_base_functions = data.getN().size2() / 3;
824  auto t_row_base_fun = data.getFTensor1N<3>();
828  auto get_ftensor1 = [](auto &v) {
829  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
830  &v[2]);
831  };
832 
833  // get bc data
834  FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
835 
836  double theta = bc.theta;
837  theta *= getFEMethod()->ts_t;
838 
840  const double a = sqrt(t_normal(i) * t_normal(i));
841  t_omega(i) = t_normal(i) * (theta / a);
842  auto t_R = get_rotation_form_vector(t_omega);
843  auto t_coords = getFTensor1CoordsAtGaussPts();
844 
845  for (int gg = 0; gg != nb_integration_pts; ++gg) {
847  t_delta(i) = t_center(i) - t_coords(i);
849  t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
850 
851  auto t_nf = get_ftensor1(nF);
852  int bb = 0;
853  for (; bb != nb_dofs / 3; ++bb) {
854  t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
855  ++t_nf;
856  ++t_row_base_fun;
857  }
858  for (; bb != nb_base_functions; ++bb)
859  ++t_row_base_fun;
860 
861  ++t_w;
862  ++t_coords;
863  }
864  }
865  }
867 }
868 
869 MoFEMErrorCode FeTractionBc::preProcess() {
871 
872  struct FaceRule {
873  int operator()(int p_row, int p_col, int p_data) const {
874  return 2 * (p_data + 1);
875  }
876  };
877 
878  if (ts_ctx == CTX_TSSETIFUNCTION) {
879 
880  // Loop boundary elements with traction boundary conditions
882  fe.getOpPtrVector().push_back(
883  new OpTractionBc(std::string("P") /* + "_RT"*/, *this));
884  fe.getRuleHook = FaceRule();
885  fe.ts_t = ts_t;
886  CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
887  fe);
888  }
889 
891 }
892 
893 MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type, EntData &data) {
895 
898 
899  auto t_normal = getFTensor1Normal();
900  const double nrm2 = sqrt(t_normal(i) * t_normal(i));
901  FTensor::Tensor1<double, 3> t_unit_normal;
902  t_unit_normal(i) = t_normal(i) / nrm2;
903  int nb_dofs = data.getFieldData().size();
904  int nb_integration_pts = data.getN().size1();
905  int nb_base_functions = data.getN().size2() / 3;
906  double ts_t = getFEMethod()->ts_t;
907 
908  auto integrate_matrix = [&]() {
910 
911  auto t_w = getFTensor0IntegrationWeight();
912  matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
913  matPP.clear();
914 
915  auto t_row_base_fun = data.getFTensor1N<3>();
916  for (int gg = 0; gg != nb_integration_pts; ++gg) {
917 
918  int rr = 0;
919  for (; rr != nb_dofs / 3; ++rr) {
920  const double a = t_row_base_fun(i) * t_unit_normal(i);
921  auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
922  for (int cc = 0; cc != nb_dofs / 3; ++cc) {
923  const double b = t_col_base_fun(i) * t_unit_normal(i);
924  matPP(rr, cc) += t_w * a * b;
925  ++t_col_base_fun;
926  }
927  ++t_row_base_fun;
928  }
929 
930  for (; rr != nb_base_functions; ++rr)
931  ++t_row_base_fun;
932 
933  ++t_w;
934  }
935 
937  };
938 
939  auto integrate_rhs = [&](auto &bc) {
941 
942  auto t_w = getFTensor0IntegrationWeight();
943  vecPv.resize(3, nb_dofs / 3, false);
944  vecPv.clear();
945 
946  auto t_row_base_fun = data.getFTensor1N<3>();
947  double ts_t = getFEMethod()->ts_t;
948 
949  for (int gg = 0; gg != nb_integration_pts; ++gg) {
950  int rr = 0;
951  for (; rr != nb_dofs / 3; ++rr) {
952  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
953  for (int dd = 0; dd != 3; ++dd)
954  if (bc.flags[dd])
955  vecPv(dd, rr) += t * bc.vals[dd];
956  ++t_row_base_fun;
957  }
958  for (; rr != nb_base_functions; ++rr)
959  ++t_row_base_fun;
960  ++t_w;
961  }
963  };
964 
965  auto integrate_rhs_cook = [&](auto &bc) {
967 
968  vecPv.resize(3, nb_dofs / 3, false);
969  vecPv.clear();
970 
971  auto t_w = getFTensor0IntegrationWeight();
972  auto t_row_base_fun = data.getFTensor1N<3>();
973  auto t_coords = getFTensor1CoordsAtGaussPts();
974 
975  for (int gg = 0; gg != nb_integration_pts; ++gg) {
976 
977  auto calc_tau = [](double y) {
978  y -= 44;
979  y /= (60 - 44);
980  return -y * (y - 1) / 0.25;
981  };
982 
983  const double tau = calc_tau(t_coords(1));
984 
985  int rr = 0;
986  for (; rr != nb_dofs / 3; ++rr) {
987  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
988  for (int dd = 0; dd != 3; ++dd)
989  if (bc.flags[dd])
990  vecPv(dd, rr) += t * tau * bc.vals[dd];
991  ++t_row_base_fun;
992  }
993 
994  for (; rr != nb_base_functions; ++rr)
995  ++t_row_base_fun;
996  ++t_w;
997  ++t_coords;
998  }
1000  };
1001 
1002  // get entity of face
1003  EntityHandle fe_ent = getFEEntityHandle();
1004  for (auto &bc : *(bcFe.bcData)) {
1005  if (bc.faces.find(fe_ent) != bc.faces.end()) {
1006 
1007  int nb_dofs = data.getFieldData().size();
1008  if (nb_dofs) {
1009 
1010  CHKERR integrate_matrix();
1011  if (bc.blockName.compare(20, 4, "COOK") == 0)
1012  CHKERR integrate_rhs_cook(bc);
1013  else
1014  CHKERR integrate_rhs(bc);
1015 
1016  const auto info =
1017  lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
1018  nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
1019  if (info > 0)
1020  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1021  "The leading minor of order %i is "
1022  "not positive; definite;\nthe "
1023  "solution could not be computed",
1024  info);
1025 
1026  for (int dd = 0; dd != 3; ++dd)
1027  if (bc.flags[dd])
1028  for (int rr = 0; rr != nb_dofs / 3; ++rr)
1029  data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
1030 
1031  // int nb_integration_pts = data.getN().size1();
1032  // auto t_row_base_fun = data.getFTensor1N<3>();
1033  // auto t_coords = getFTensor1CoordsAtGaussPts();
1034  // for (int gg = 0; gg != nb_integration_pts; ++gg) {
1035 
1036  // FTensor::Tensor1<double, 3> t_v;
1037  // t_v(i) = 0;
1038 
1039  // int rr = 0;
1040  // for (; rr != nb_dofs / 3; ++rr) {
1041  // for (int dd = 0; dd != 3; ++dd) {
1042  // const double v = data.getFieldDofs()[3 * rr +
1043  // dd]->getFieldData(); t_v(dd) += t_row_base_fun(i) *
1044  // t_unit_normal(i) * v;
1045  // }
1046  // ++t_row_base_fun;
1047  // }
1048  // for (; rr != nb_base_functions; ++rr)
1049  // ++t_row_base_fun;
1050  // ++t_coords;
1051 
1052  // cerr << t_v << endl;
1053  // }
1054  }
1055  }
1056  }
1057 
1059 }
1060 
1061 MoFEMErrorCode OpSpatialEquilibrium_dw_dP::integrate(EntData &row_data,
1062  EntData &col_data) {
1064  int nb_integration_pts = row_data.getN().size1();
1065  int row_nb_dofs = row_data.getIndices().size();
1066  int col_nb_dofs = col_data.getIndices().size();
1067  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1069  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
1070  };
1072  auto v = getVolume();
1073  auto t_w = getFTensor0IntegrationWeight();
1074  int row_nb_base_functions = row_data.getN().size2();
1075  auto t_row_base_fun = row_data.getFTensor0N();
1076  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1077  double a = v * t_w;
1078  int rr = 0;
1079  for (; rr != row_nb_dofs / 3; ++rr) {
1080  auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
1081  auto t_m = get_ftensor1(K, 3 * rr, 0);
1082  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1083  double div_col_base = t_col_diff_base_fun(i, i);
1084  t_m(i) += a * t_row_base_fun * div_col_base;
1085  ++t_m;
1086  ++t_col_diff_base_fun;
1087  }
1088  ++t_row_base_fun;
1089  }
1090  for (; rr != row_nb_base_functions; ++rr)
1091  ++t_row_base_fun;
1092  ++t_w;
1093  }
1095 }
1096 
1097 MoFEMErrorCode OpSpatialEquilibrium_dw_dw::integrate(EntData &row_data,
1098  EntData &col_data) {
1100 
1101  if (alpha_w < std::numeric_limits<double>::epsilon())
1103 
1104  const int nb_integration_pts = row_data.getN().size1();
1105  const int row_nb_dofs = row_data.getIndices().size();
1106  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1108  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
1109 
1110  );
1111  };
1113 
1114  auto v = getVolume();
1115  auto t_w = getFTensor0IntegrationWeight();
1116 
1117  int row_nb_base_functions = row_data.getN().size2();
1118  auto t_row_base_fun = row_data.getFTensor0N();
1119 
1120  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1121  double a = v * t_w * alpha_w * getTSa();
1122 
1123  int rr = 0;
1124  for (; rr != row_nb_dofs / 3; ++rr) {
1125 
1126  auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
1127  auto t_m = get_ftensor1(K, 3 * rr, 0);
1128  for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1129  const double b = a * t_row_base_fun * t_col_base_fun;
1130  t_m(i) -= b;
1131  ++t_m;
1132  ++t_col_base_fun;
1133  }
1134 
1135  ++t_row_base_fun;
1136  }
1137 
1138  for (; rr != row_nb_base_functions; ++rr)
1139  ++t_row_base_fun;
1140 
1141  ++t_w;
1142  }
1143 
1145 }
1146 
1147 MoFEMErrorCode OpSpatialPhysical_du_dP::integrate(EntData &row_data,
1148  EntData &col_data) {
1150  int nb_integration_pts = row_data.getN().size1();
1151  int row_nb_dofs = row_data.getIndices().size();
1152  int col_nb_dofs = col_data.getIndices().size();
1153  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1155 
1156  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1157 
1158  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1159 
1160  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1161 
1162  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1163 
1164  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1165 
1166  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
1167  };
1168 
1174 
1175  auto v = getVolume();
1176  auto t_w = getFTensor0IntegrationWeight();
1177  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1178  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
1179  &dataAtPts->diffStreachTensorAtPts(0, 0), nb_integration_pts);
1180 
1181  int row_nb_base_functions = row_data.getN().size2();
1182  auto t_row_base_fun = row_data.getFTensor0N();
1183  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1184  double a = v * t_w;
1185 
1186  int rr = 0;
1187  for (; rr != row_nb_dofs / 6; ++rr) {
1188 
1189  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1190  auto t_m = get_ftensor3(K, 6 * rr, 0);
1191 
1192  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1193 
1195  t_RTP_dP(i, j, k) = t_R(k, i) * t_col_base_fun(j) * t_row_base_fun;
1196 
1197  for (int kk = 0; kk != 3; ++kk)
1198  for (int ii = 0; ii != 3; ++ii)
1199  for (int jj = 0; jj <= ii; ++jj) {
1200  for (int mm = 0; mm != 3; ++mm)
1201  for (int nn = 0; nn != 3; ++nn)
1202  t_m(ii, jj, kk) +=
1203  a * t_diff_u(mm, nn, ii, jj) * t_RTP_dP(mm, nn, kk);
1204  }
1205 
1206  ++t_col_base_fun;
1207  ++t_m;
1208  }
1209 
1210  ++t_row_base_fun;
1211  }
1212  for (; rr != row_nb_base_functions; ++rr)
1213  ++t_row_base_fun;
1214  ++t_w;
1215  ++t_R;
1216  ++t_diff_u;
1217  }
1219 }
1220 
1221 MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
1222  EntData &col_data) {
1224  int nb_integration_pts = row_data.getN().size1();
1225  int row_nb_dofs = row_data.getIndices().size();
1226  int col_nb_dofs = col_data.getIndices().size();
1227  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1229  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
1230  &m(r + 5, c));
1231  };
1232 
1236  auto v = getVolume();
1237  auto t_w = getFTensor0IntegrationWeight();
1238  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1239  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
1240  &dataAtPts->diffStreachTensorAtPts(0, 0), nb_integration_pts);
1241 
1242  int row_nb_base_functions = row_data.getN().size2();
1243  auto t_row_base_fun = row_data.getFTensor0N();
1244  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1245  double a = v * t_w;
1246 
1247  int rr = 0;
1248  for (; rr != row_nb_dofs / 6; ++rr) {
1249 
1250  auto t_m = get_ftensor2(K, 6 * rr, 0);
1251 
1252  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1253  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1254 
1255  FTensor::Tensor2<double, 3, 3> t_RTP_dBubble;
1256  t_RTP_dBubble(i, j) =
1257  a * t_row_base_fun * t_R(k, i) * t_col_base_fun(k, j);
1258 
1259  for (int ii = 0; ii != 3; ++ii)
1260  for (int jj = 0; jj <= ii; ++jj)
1261  for (int mm = 0; mm != 3; ++mm)
1262  for (int nn = 0; nn != 3; ++nn)
1263  t_m(ii, jj) += t_diff_u(mm, nn, ii, jj) * t_RTP_dBubble(mm, nn);
1264 
1265  ++t_m;
1266  ++t_col_base_fun;
1267  }
1268 
1269  ++t_row_base_fun;
1270  }
1271  for (; rr != row_nb_base_functions; ++rr)
1272  ++t_row_base_fun;
1273  ++t_w;
1274  ++t_R;
1275  ++t_diff_u;
1276  }
1278 }
1279 
1280 MoFEMErrorCode OpSpatialPhysical_du_du::integrate(EntData &row_data,
1281  EntData &col_data) {
1283 
1284  int nb_integration_pts = row_data.getN().size1();
1285  int row_nb_dofs = row_data.getIndices().size();
1286  int col_nb_dofs = col_data.getIndices().size();
1287  auto get_ftensor4 = [](MatrixDouble &m, const int r, const int c) {
1289 
1290  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1291  &m(r + 0, c + 4), &m(r + 0, c + 5),
1292 
1293  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1294  &m(r + 1, c + 4), &m(r + 1, c + 5),
1295 
1296  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1297  &m(r + 2, c + 4), &m(r + 2, c + 5),
1298 
1299  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1300  &m(r + 3, c + 4), &m(r + 3, c + 5),
1301 
1302  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1303  &m(r + 4, c + 4), &m(r + 4, c + 5),
1304 
1305  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1306  &m(r + 5, c + 4), &m(r + 5, c + 5)
1307 
1308  );
1309  };
1316 
1320 
1321  auto v = getVolume();
1322  auto t_w = getFTensor0IntegrationWeight();
1323  auto t_P_dh0 = getFTensor3FromMat(dataAtPts->P_dh0);
1324  auto t_P_dh1 = getFTensor3FromMat(dataAtPts->P_dh1);
1325  auto t_P_dh2 = getFTensor3FromMat(dataAtPts->P_dh2);
1326  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1327  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1328  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
1329  auto t_dot_log_u =
1330  getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachDotTensorAtPts);
1331  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1332  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
1333  &dataAtPts->diffStreachTensorAtPts(0, 0),
1334  dataAtPts->diffStreachTensorAtPts.size2());
1335  std::array<FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>, 6>
1336  t_diff2_u = {
1337 
1339  &dataAtPts->expLogUHessian[0](0, 0),
1340  dataAtPts->expLogUHessian[0].size2()),
1342  &dataAtPts->expLogUHessian[1](0, 0),
1343  dataAtPts->expLogUHessian[1].size2()),
1345  &dataAtPts->expLogUHessian[2](0, 0),
1346  dataAtPts->expLogUHessian[2].size2()),
1348  &dataAtPts->expLogUHessian[3](0, 0),
1349  dataAtPts->expLogUHessian[3].size2()),
1351  &dataAtPts->expLogUHessian[4](0, 0),
1352  dataAtPts->expLogUHessian[4].size2()),
1354  &dataAtPts->expLogUHessian[5](0, 0),
1355  dataAtPts->expLogUHessian[5].size2())
1356 
1357  };
1358 
1359  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
1361  t_one4(i, j, k, l) = t_kd(j, l) * t_kd(i, k);
1362 
1363  FTensor::Tensor4<double, 3, 3, 3, 3> t_one4_symmetric;
1364  t_one4_symmetric(i, j, k, l) = 0;
1365  for (auto ii : {0, 1, 2})
1366  for (auto jj : {0, 1, 2})
1367  for (auto kk : {0, 1, 2})
1368  for (auto ll : {0, 1, 2}) {
1369 
1370  if (ll != kk)
1371  t_one4_symmetric(ii, jj, kk, ll) =
1372  t_one4(ii, jj, kk, ll) + t_one4(ii, jj, ll, kk);
1373  else
1374  t_one4_symmetric(ii, jj, kk, ll) = t_one4(ii, jj, kk, ll);
1375  }
1376 
1377  int row_nb_base_functions = row_data.getN().size2();
1378  auto t_row_base_fun = row_data.getFTensor0N();
1379 
1381  t_one(i, j) = 0;
1382  for (auto ii : {0, 1, 2})
1383  t_one(ii, ii) = 1;
1384 
1385  const double ts_a = getTSa();
1386  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1387  double a = v * t_w;
1388 
1390  t_P_dh(i, j, N0, k) = t_P_dh0(i, j, k);
1391  t_P_dh(i, j, N1, k) = t_P_dh1(i, j, k);
1392  t_P_dh(i, j, N2, k) = t_P_dh2(i, j, k);
1393 
1395  t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
1397  t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
1398 
1400  t_dot_u(i, j) = t_u(i, k) * t_dot_log_u(k, j);
1401 
1403  t_stress_diff(i, j, k, l) = 0;
1404  for (auto ii : {0, 1, 2})
1405  for (auto jj : {0, 1, 2})
1406  for (auto kk : {0, 1, 2})
1407  for (auto ll : {0, 1, 2})
1408  for (auto mm : {0, 1, 2})
1409  for (auto nn : {0, 1, 2})
1410  t_stress_diff(ii, jj, mm, nn) -=
1411  t_P_dh(ii, jj, kk, ll) * t_diff_u(kk, ll, mm, nn);
1412 
1413  // FTensor::Tensor4<double, 3, 3, 3, 3> t_dot_u_diff;
1414  // t_dot_u_diff(i, j, k, l) = 0;
1415  // for (auto ii : {0, 1, 2})
1416  // for (auto jj : {0, 1, 2})
1417  // for (auto kk : {0, 1, 2})
1418  // for (auto ll : {0, 1, 2})
1419  // for (auto mm : {0, 1, 2})
1420  // t_dot_u_diff(ii, jj, kk, ll) +=
1421  // t_diff_u(ii, mm, kk, ll) * t_dot_log_u(mm, jj) +
1422  // ts_a * t_u(ii, mm) * t_one4_symmetric(mm, jj, kk, ll);
1423 
1424  // auto t_viscous_stress = calculate_viscous_stress(t_dot_u);
1425  // t_viscous_stress(i, j) = t_dot_u(i, j);
1426  // t_viscous_stress(i, j) *= alpha_u * t_dot_log_u(i, j);
1427 
1428  // auto t_viscous_stress_diff =
1429  // calculate_diff_viscous_stress(t_dot_u, t_dot_u_diff);
1430  // t_viscous_stress_diff(i, j, k, l) = t_dot_u_diff(i, j, k, l);
1431  // t_viscous_stress_diff(i, j, k, l) *= alpha_u * t_one4(i, j, k, l);
1432 
1433  FTensor::Tensor2_symmetric<double, 3> t_viscous_stress;
1434  t_viscous_stress(i, j) = alpha_u * t_dot_log_u(i, j);
1435  FTensor::Tensor4<double, 3, 3, 3, 3> t_viscous_stress_diff;
1436  t_viscous_stress_diff(i, j, k, l) = t_one4_symmetric(i, j, k, l);
1437  t_viscous_stress_diff(i, j, k, l) *= (ts_a * alpha_u);
1438 
1439  int rr = 0;
1440  for (; rr != row_nb_dofs / 6; ++rr) {
1441 
1442  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1443  auto t_m = get_ftensor4(K, 6 * rr, 0);
1444 
1445  for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1446 
1447  const double b = a * t_row_base_fun * t_col_base_fun;
1448 
1449  for (int ii : {0, 1, 2})
1450  for (int jj = 0; jj <= ii; ++jj)
1451  for (int kk : {0, 1, 2})
1452  for (int ll = 0; ll <= kk; ++ll)
1453  for (int mm : {0, 1, 2})
1454  for (int nn : {0, 1, 2})
1455  t_m(ii, jj, kk, ll) +=
1456  b * t_diff_u(mm, nn, ii, jj) *
1457  (t_stress_diff(mm, nn, kk, ll) -
1458  t_viscous_stress_diff(mm, nn, kk, ll));
1459 
1460  for (int ii : {0, 1, 2})
1461  for (int jj = 0; jj <= ii; ++jj) {
1462  int ss = 0;
1463  for (int kk : {0, 1, 2})
1464  for (int ll = 0; ll <= kk; ++ll, ++ss)
1465  for (int mm : {0, 1, 2})
1466  for (int nn : {0, 1, 2})
1467  t_m(ii, jj, kk, ll) +=
1468  b * (t_diff2_u[ss])(mm, nn, ii, jj) *
1469  (t_deltaP(mm, nn) - t_viscous_stress(mm, nn));
1470  }
1471 
1472  ++t_m;
1473  ++t_col_base_fun;
1474  }
1475 
1476  ++t_row_base_fun;
1477  }
1478  for (; rr != row_nb_base_functions; ++rr)
1479  ++t_row_base_fun;
1480 
1481  ++t_w;
1482  ++t_P_dh0;
1483  ++t_P_dh1;
1484  ++t_P_dh2;
1485  ++t_R;
1486  ++t_approx_P;
1487  ++t_P;
1488  ++t_dot_log_u;
1489  ++t_u;
1490  ++t_diff_u;
1491  for (auto ll : {0, 1, 2, 3, 4, 5})
1492  ++(t_diff2_u[ll]);
1493  }
1495 }
1496 
1497 MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
1498  EntData &col_data) {
1500  int row_nb_dofs = row_data.getIndices().size();
1501  int col_nb_dofs = col_data.getIndices().size();
1502  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1504 
1505  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1506 
1507  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1508 
1509  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1510 
1511  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1512 
1513  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1514 
1515  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1516 
1517  );
1518  };
1523 
1524  auto v = getVolume();
1525  auto t_w = getFTensor0IntegrationWeight();
1526  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1527  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1528  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
1529  &dataAtPts->diffStreachTensorAtPts(0, 0),
1530  dataAtPts->diffStreachTensorAtPts.size2());
1531 
1532  int row_nb_base_functions = row_data.getN().size2();
1533  auto t_row_base_fun = row_data.getFTensor0N();
1534 
1535  int nb_integration_pts = row_data.getN().size1();
1536 
1537  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1538  double a = v * t_w;
1539 
1540  FTensor::Tensor3<double, 3, 3, 3> t_RT_domega_P;
1541  t_RT_domega_P(i, j, m) = t_diff_R(k, i, m) * t_P(k, j);
1543  t(i, j, k) = 0;
1544  for (int ii = 0; ii != 3; ++ii)
1545  for (int jj = 0; jj <= ii; ++jj)
1546  for (int mm = 0; mm != 3; ++mm)
1547  for (int nn = 0; nn != 3; ++nn)
1548  for (int kk = 0; kk != 3; ++kk)
1549  t(ii, jj, kk) +=
1550  t_diff_u(mm, nn, ii, jj) * t_RT_domega_P(mm, nn, kk);
1551 
1552  int rr = 0;
1553  for (; rr != row_nb_dofs / 6; ++rr) {
1554  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1555  auto t_m = get_ftensor3(K, 6 * rr, 0);
1556  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1557 
1558  double v = a * t_row_base_fun * t_col_base_fun;
1559  t_m(i, j, k) += v * t(i, j, k);
1560 
1561  ++t_m;
1562  ++t_col_base_fun;
1563  }
1564 
1565  ++t_row_base_fun;
1566  }
1567 
1568  for (; rr != row_nb_base_functions; ++rr)
1569  ++t_row_base_fun;
1570 
1571  ++t_w;
1572  ++t_P;
1573  ++t_diff_R;
1574  ++t_diff_u;
1575  }
1576 
1578 }
1579 
1580 MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
1581  EntData &col_data) {
1583  int nb_integration_pts = row_data.getN().size1();
1584  int row_nb_dofs = row_data.getIndices().size();
1585  int col_nb_dofs = col_data.getIndices().size();
1586  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1588 
1589  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1590 
1591  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1592 
1593  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1594 
1595  );
1596  };
1603 
1604  auto v = getVolume();
1605  auto t_w = getFTensor0IntegrationWeight();
1606  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1607 
1608  int row_nb_base_functions = row_data.getN().size2();
1609  auto t_row_base_fun = row_data.getFTensor0N();
1610 
1611  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1612  double a = v * t_w;
1613 
1614  int rr = 0;
1615  for (; rr != row_nb_dofs / 3; ++rr) {
1616  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1617  auto t_m = get_ftensor2(K, 3 * rr, 0);
1619  t_levi_R(i, j, k) =
1620  (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1621  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1622  FTensor::Tensor1<double, 3> t_pushed_base;
1623  t_m(k, i) += t_levi_R(i, j, k) * t_col_base_fun(j);
1624  ++t_m;
1625  ++t_col_base_fun;
1626  }
1627 
1628  ++t_row_base_fun;
1629  }
1630 
1631  for (; rr != row_nb_base_functions; ++rr)
1632  ++t_row_base_fun;
1633  ++t_w;
1634  ++t_R;
1635  }
1637 }
1638 
1639 MoFEMErrorCode OpSpatialRotation_domega_dBubble::integrate(EntData &row_data,
1640  EntData &col_data) {
1642  int nb_integration_pts = row_data.getN().size1();
1643  int row_nb_dofs = row_data.getIndices().size();
1644  int col_nb_dofs = col_data.getIndices().size();
1645  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1647  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
1648  };
1654 
1655  auto v = getVolume();
1656  auto t_w = getFTensor0IntegrationWeight();
1657  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1658 
1659  int row_nb_base_functions = row_data.getN().size2();
1660  auto t_row_base_fun = row_data.getFTensor0N();
1661 
1662  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1663  double a = v * t_w;
1664 
1665  int rr = 0;
1666  for (; rr != row_nb_dofs / 3; ++rr) {
1667 
1668  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1669  auto t_m = get_ftensor1(K, 3 * rr, 0);
1670 
1672  t_levi_r(i, j, k) =
1673  (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1674 
1675  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1676  FTensor::Tensor2<double, 3, 3> t_pushed_base;
1677  t_m(k) += t_levi_r(i, j, k) * t_col_base_fun(i, j);
1678  ++t_m;
1679  ++t_col_base_fun;
1680  }
1681 
1682  ++t_row_base_fun;
1683  }
1684 
1685  for (; rr != row_nb_base_functions; ++rr)
1686  ++t_row_base_fun;
1687  ++t_w;
1688  ++t_R;
1689  }
1691 }
1692 
1693 MoFEMErrorCode OpSpatialRotation_domega_domega::integrate(EntData &row_data,
1694  EntData &col_data) {
1696  int nb_integration_pts = row_data.getN().size1();
1697  int row_nb_dofs = row_data.getIndices().size();
1698  int col_nb_dofs = col_data.getIndices().size();
1699  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1701 
1702  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1703 
1704  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1705 
1706  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1707 
1708  );
1709  };
1716 
1717  auto v = getVolume();
1718  auto t_w = getFTensor0IntegrationWeight();
1719  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1720  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1721  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1722 
1723  int row_nb_base_functions = row_data.getN().size2();
1724  auto t_row_base_fun = row_data.getFTensor0N();
1725 
1726  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1727  double a = v * t_w;
1728 
1730  t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
1731 
1732  FTensor::Tensor2<double, 3, 3> t_leviPRT_domega;
1733  t_leviPRT_domega(n, l) =
1734  levi_civita(i, j, n) * (t_approx_P(i, k) * t_diff_R(j, k, l));
1735 
1736  int rr = 0;
1737  for (; rr != row_nb_dofs / 3; ++rr) {
1738 
1739  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1740  auto t_m = get_ftensor2(K, 3 * rr, 0);
1741  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1742 
1743  const double b = a * t_row_base_fun * t_col_base_fun;
1744 
1745  t_m(k, m) += b * t_leviPRT_domega(k, m);
1746 
1747  ++t_m;
1748  ++t_col_base_fun;
1749  }
1750 
1751  ++t_row_base_fun;
1752  }
1753 
1754  for (; rr != row_nb_base_functions; ++rr)
1755  ++t_row_base_fun;
1756  ++t_w;
1757  ++t_approx_P;
1758  ++t_R;
1759  ++t_diff_R;
1760  }
1762 }
1763 
1764 MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
1765  EntData &col_data) {
1767 
1768  int nb_integration_pts = row_data.getN().size1();
1769  int row_nb_dofs = row_data.getIndices().size();
1770  int col_nb_dofs = col_data.getIndices().size();
1771 
1772  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1774 
1775  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1776 
1777  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1778 
1779  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1780 
1781  );
1782  };
1783 
1789 
1790  auto v = getVolume();
1791  auto t_w = getFTensor0IntegrationWeight();
1792  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1793  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1794  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1795  int row_nb_base_functions = row_data.getN().size2() / 3;
1796  auto t_row_base_fun = row_data.getFTensor1N<3>();
1797 
1798  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1799  double a = v * t_w;
1800 
1801  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
1802 
1803  int rr = 0;
1804  for (; rr != row_nb_dofs / 3; ++rr) {
1805 
1807  t_PRT(i, k) = t_row_base_fun(j) * (t_diff_R(i, l, k) * t_u(l, j));
1808 
1809  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1810  auto t_m = get_ftensor2(K, 3 * rr, 0);
1811  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1812  t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
1813  ++t_m;
1814  ++t_col_base_fun;
1815  }
1816 
1817  ++t_row_base_fun;
1818  }
1819 
1820  for (; rr != row_nb_base_functions; ++rr)
1821  ++t_row_base_fun;
1822  ++t_w;
1823  ++t_R;
1824  ++t_diff_R;
1825  ++t_u;
1826  }
1828 }
1829 
1831 OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
1832  EntData &col_data) {
1834 
1835  int nb_integration_pts = row_data.getN().size1();
1836  int row_nb_dofs = row_data.getIndices().size();
1837  int col_nb_dofs = col_data.getIndices().size();
1838 
1839  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1841  &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
1842  };
1843 
1849 
1850  auto v = getVolume();
1851  auto t_w = getFTensor0IntegrationWeight();
1852  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1853  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1854  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1855  int row_nb_base_functions = row_data.getN().size2() / 9;
1856  auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
1857 
1858  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1859  double a = v * t_w;
1860 
1861  int rr = 0;
1862  for (; rr != row_nb_dofs; ++rr) {
1863 
1865  t_PRT(k) = t_row_base_fun(i, j) * (t_diff_R(i, l, k) * t_u(l, j));
1866 
1867  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1868  auto t_m = get_ftensor2(K, rr, 0);
1869  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1870  t_m(j) += (a * t_col_base_fun) * t_PRT(j);
1871  ++t_m;
1872  ++t_col_base_fun;
1873  }
1874 
1875  ++t_row_base_fun;
1876  }
1877 
1878  for (; rr != row_nb_base_functions; ++rr)
1879  ++t_row_base_fun;
1880 
1881  ++t_w;
1882  ++t_R;
1883  ++t_diff_R;
1884  ++t_u;
1885  }
1887 }
1888 
1889 MoFEMErrorCode OpSpatialSchurBegin::assemble(int row_side, EntityType row_type,
1890  EntData &data) {
1892  if (row_type != MBTET)
1894  auto &bmc = dataAtPts->blockMatContainor;
1895  for (auto &bit : bmc)
1896  bit.unSetAtElement();
1897  // bmc.clear();
1899 }
1900 
1901 MoFEMErrorCode OpSpatialPreconditionMass::integrate(EntData &row_data) {
1903  const int nb_integration_pts = row_data.getN().size1();
1904  const int row_nb_dofs = row_data.getIndices().size();
1905  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1907  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
1908 
1909  );
1910  };
1912 
1913  auto v = getVolume();
1914  auto t_w = getFTensor0IntegrationWeight();
1915  auto &m = *mPtr;
1916 
1917  m.resize(row_nb_dofs, row_nb_dofs, false);
1918  m.clear();
1919 
1920  int row_nb_base_functions = row_data.getN().size2();
1921  auto t_row_base_fun = row_data.getFTensor0N();
1922 
1923  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1924  double a = v * t_w;
1925 
1926  int rr = 0;
1927  for (; rr != row_nb_dofs / 3; ++rr) {
1928 
1929  auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
1930  auto t_m = get_ftensor1(m, 3 * rr, 0);
1931  for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1932  const double b = a * t_row_base_fun * t_col_base_fun;
1933  t_m(i) += b;
1934  ++t_m;
1935  ++t_col_base_fun;
1936  }
1937 
1938  ++t_row_base_fun;
1939  }
1940 
1941  for (; rr != row_nb_base_functions; ++rr)
1942  ++t_row_base_fun;
1943 
1944  ++t_w;
1945  }
1946 
1948 }
1949 
1950 MoFEMErrorCode OpSpatialSchurEnd::assemble(int row_side, EntityType row_type,
1951  EntData &data) {
1953  if (row_type != MBTET)
1955 
1956  if (auto ep_fe_vol_ptr =
1957  dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
1958  getFEMethod())) {
1959 
1960  AO aoSuu = PETSC_NULL;
1961  AO aoSBB = PETSC_NULL;
1962  AO aoSOO = PETSC_NULL;
1963  AO aoSww = PETSC_NULL;
1964 
1965  Mat Suu = PETSC_NULL;
1966  Mat SBB = PETSC_NULL;
1967  Mat SOO = PETSC_NULL;
1968  Mat Sww = PETSC_NULL;
1969 
1970  auto set_data = [&](auto fe) {
1971  aoSuu = fe->aoSuu;
1972  aoSBB = fe->aoSBubble;
1973  aoSOO = fe->aoSOmega;
1974  aoSww = fe->aoSw;
1975 
1976  Suu = fe->Suu;
1977  SBB = fe->SBubble;
1978  SOO = fe->SOmega;
1979  Sww = fe->Sw;
1980  };
1981  set_data(ep_fe_vol_ptr);
1982 
1983  if (Suu) {
1984 
1985  auto find_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1986  auto &row_name, auto &col_name, auto row_side,
1987  auto col_side, auto row_type, auto col_type) {
1988  return add_bmc.get<0>().find(boost::make_tuple(
1989  row_name, col_name, row_type, col_type, row_side, col_side));
1990  };
1991 
1992  auto set_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1993  auto &row_name, auto &col_name, auto row_side,
1994  auto col_side, auto row_type, auto col_type,
1995  const auto &m, const auto &row_ind,
1996  const auto &col_ind) {
1997  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
1998  row_type, col_type);
1999  if (it != add_bmc.get<0>().end()) {
2000  it->setMat(m);
2001  it->setInd(row_ind, col_ind);
2002  it->setSetAtElement();
2003  return it;
2004  } else {
2005  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2006  row_name, col_name, row_type, col_type, row_side, col_side, m,
2007  row_ind, col_ind));
2008  return p_it.first;
2009  }
2010  };
2011 
2012  auto add_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2013  auto &row_name, auto &col_name, auto row_side,
2014  auto col_side, auto row_type, auto col_type,
2015  const auto &m, const auto &row_ind,
2016  const auto &col_ind) {
2017  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2018  row_type, col_type);
2019  if (it != add_bmc.get<0>().end()) {
2020  it->addMat(m);
2021  it->setInd(row_ind, col_ind);
2022  it->setSetAtElement();
2023  return it;
2024  } else {
2025  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2026  row_name, col_name, row_type, col_type, row_side, col_side, m,
2027  row_ind, col_ind));
2028  return p_it.first;
2029  }
2030  };
2031 
2032  auto assemble_block = [&](auto &bit, Mat S) {
2034  const VectorInt &rind = bit.rowInd;
2035  const VectorInt &cind = bit.colInd;
2036  const MatrixDouble &m = bit.M;
2037  CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
2038  &*cind.begin(), &*m.data().begin(), ADD_VALUES);
2039 
2041  };
2042 
2043  auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
2045  const int nb = m.size1();
2046 
2047  inv.resize(nb, nb, false);
2048  inv.clear();
2049  for (int cc = 0; cc != nb; ++cc)
2050  inv(cc, cc) = -1;
2051 
2052  iPIV.resize(nb, false);
2053  lapackWork.resize(nb * nb, false);
2054  const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
2055  &*iPIV.begin(), &*inv.data().begin(), nb,
2056  &*lapackWork.begin(), nb * nb);
2057  // if (info != 0)
2058  // SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2059  // "Can not invert matrix info = %d", info);
2060 
2062  };
2063 
2064  auto invert_symm_schur = [&](DataAtIntegrationPts::BlockMatContainor &bmc,
2065  std::string field, auto &inv) {
2067 
2068  auto bit =
2069  bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
2070  if (bit != bmc.get<1>().end()) {
2071 
2072  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2073  CHKERR invert_symm_mat(m, inv);
2074 
2075  } else
2076  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2077  "%s matrix not found", field.c_str());
2078 
2080  };
2081 
2082  auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
2084 
2085  const int nb = m.size1();
2086 
2087  MatrixDouble trans_m = trans(m);
2088  MatrixDouble trans_inv;
2089  trans_inv.resize(nb, nb, false);
2090  trans_inv.clear();
2091  for (int c = 0; c != nb; ++c)
2092  trans_inv(c, c) = -1;
2093 
2094  iPIV.resize(nb, false);
2095  const auto info =
2096  lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb, &*iPIV.begin(),
2097  &*trans_inv.data().begin(), nb);
2098  if (info != 0)
2099  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2100  "Can not invert matrix info = %d", info);
2101 
2102  inv.resize(nb, nb, false);
2103  noalias(inv) = trans(trans_inv);
2104 
2106  };
2107 
2108  auto invert_nonsymm_schur =
2109  [&](DataAtIntegrationPts::BlockMatContainor &bmc, std::string field,
2110  auto &inv, const bool debug = false) {
2112 
2113  auto bit = bmc.get<1>().find(
2114  boost::make_tuple(field, field, MBTET, MBTET));
2115  if (bit != bmc.get<1>().end()) {
2116 
2117  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2118  CHKERR invert_nonsymm_mat(m, inv);
2119 
2120  if (debug) {
2121  std::cerr << prod(m, inv) << endl;
2122  std::cerr << endl;
2123  }
2124 
2125  } else
2126  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2127  "%s matrix not found", field.c_str());
2128 
2130  };
2131 
2132  auto create_block_schur =
2135  std::string field, AO ao, MatrixDouble &inv) {
2137 
2138  for (auto &bit : add_bmc) {
2139  bit.unSetAtElement();
2140  bit.clearMat();
2141  }
2142 
2143  for (auto &bit : bmc) {
2144  if (bit.setAtElement && bit.rowField != field &&
2145  bit.colField != field) {
2146  VectorInt rind = bit.rowInd;
2147  VectorInt cind = bit.colInd;
2148  const MatrixDouble &m = bit.M;
2149  if (ao) {
2150  CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
2151  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2152  }
2153  auto it = set_block(add_bmc, bit.rowField, bit.colField,
2154  bit.rowSide, bit.colSide, bit.rowType,
2155  bit.colType, m, rind, cind);
2156  }
2157  }
2158 
2159  for (auto &bit_col : bmc) {
2160  if (bit_col.setAtElement && bit_col.rowField == field &&
2161  bit_col.colField != field) {
2162  const MatrixDouble &cm = bit_col.M;
2163  VectorInt cind = bit_col.colInd;
2164  invMat.resize(inv.size1(), cm.size2(), false);
2165  noalias(invMat) = prod(inv, cm);
2166  if (ao)
2167  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2168  for (auto &bit_row : bmc) {
2169  if (bit_row.setAtElement && bit_row.rowField != field &&
2170  bit_row.colField == field) {
2171  const MatrixDouble &rm = bit_row.M;
2172  VectorInt rind = bit_row.rowInd;
2173  K.resize(rind.size(), cind.size(), false);
2174  noalias(K) = prod(rm, invMat);
2175  if (ao)
2176  CHKERR AOApplicationToPetsc(ao, rind.size(),
2177  &*rind.begin());
2178  auto it = add_block(add_bmc, bit_row.rowField,
2179  bit_col.colField, bit_row.rowSide,
2180  bit_col.colSide, bit_row.rowType,
2181  bit_col.colType, K, rind, cind);
2182  }
2183  }
2184  }
2185  }
2186 
2188  };
2189 
2190  auto assemble_schur =
2191  [&](DataAtIntegrationPts::BlockMatContainor &add_bmc, Mat S,
2192  bool debug = false) {
2194  for (auto &bit : add_bmc) {
2195  if (bit.setAtElement)
2196  CHKERR assemble_block(bit, S);
2197  }
2198  if (debug) {
2199  for (auto &bit : add_bmc) {
2200  if (bit.setAtElement) {
2201  std::cerr << "assemble: " << bit.rowField << " "
2202  << bit.colField << endl;
2203  std::cerr << bit.M << endl;
2204  }
2205  }
2206  std::cerr << std::endl;
2207  }
2209  };
2210 
2211  auto precondition_schur =
2214  const std::string field, const MatrixDouble &diag_mat,
2215  const double eps) {
2217 
2218  for (auto &bit : add_bmc) {
2219  bit.unSetAtElement();
2220  bit.clearMat();
2221  }
2222 
2223  for (auto &bit : bmc) {
2224  if (bit.setAtElement) {
2225  if (bit.rowField != field || bit.colField != field)
2226  auto it =
2227  set_block(add_bmc, bit.rowField, bit.colField,
2228  bit.rowSide, bit.colSide, bit.rowType,
2229  bit.colType, bit.M, bit.rowInd, bit.colInd);
2230  }
2231  }
2232 
2233  auto bit = bmc.get<1>().find(
2234  boost::make_tuple(field, field, MBTET, MBTET));
2235  if (bit->setAtElement && bit != bmc.get<1>().end()) {
2236  auto it =
2237  set_block(add_bmc, bit->rowField, bit->colField, bit->rowSide,
2238  bit->colSide, bit->rowType, bit->colType, bit->M,
2239  bit->rowInd, bit->colInd);
2240  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2241  m += eps * diag_mat;
2242  } else {
2243  auto row_it = bmc.get<3>().lower_bound(field);
2244  for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
2245  if (row_it->setAtElement) {
2246  auto it = set_block(add_bmc, field, field, 0, 0, MBTET, MBTET,
2247  diag_mat, row_it->rowInd, row_it->rowInd);
2248  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2249  m *= eps;
2250  break;
2251  }
2252  }
2253  if (row_it == bmc.get<3>().end())
2254  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2255  "row field not found %s", field.c_str());
2256  }
2257 
2259  };
2260 
2261  CHKERR invert_symm_schur(dataAtPts->blockMatContainor, "u",
2262  invBlockMat["uu"]);
2263  CHKERR create_block_schur(dataAtPts->blockMatContainor, blockMat["uu"],
2264  "u", aoSuu, invBlockMat["uu"]);
2265  CHKERR assemble_schur(blockMat["uu"], Suu);
2266 
2267  if (SBB) {
2268  CHKERR invert_symm_schur(blockMat["uu"], "BUBBLE", invBlockMat["BB"]);
2269  CHKERR create_block_schur(blockMat["uu"], blockMat["BB"], "BUBBLE",
2270  aoSBB, invBlockMat["BB"]);
2271  CHKERR precondition_schur(blockMat["BB"], blockMat["precBBOO"], "omega",
2272  *dataAtPts->ooMatPtr, eps);
2273  CHKERR assemble_schur(blockMat["precBBOO"], SBB);
2274 
2275  if (SOO) {
2276  CHKERR invert_nonsymm_schur(blockMat["precBBOO"], "omega",
2277  invBlockMat["OO"]);
2278  CHKERR create_block_schur(blockMat["precBBOO"], blockMat["OO"],
2279  "omega", aoSOO, invBlockMat["OO"]);
2280  if (dataAtPts->wwMatPtr) {
2281  CHKERR precondition_schur(blockMat["OO"], blockMat["precOOww"], "w",
2282  *dataAtPts->wwMatPtr, -eps);
2283  } else {
2284  blockMat["precOOww"] = blockMat["OO"];
2285  }
2286  CHKERR assemble_schur(blockMat["precOOww"], SOO);
2287 
2288  if (Sww) {
2289  CHKERR invert_symm_schur(blockMat["precOOww"], "w",
2290  invBlockMat["ww"]);
2291  CHKERR create_block_schur(blockMat["precOOww"], blockMat["ww"], "w",
2292  aoSww, invBlockMat["ww"]);
2293  CHKERR assemble_schur(blockMat["ww"], Sww);
2294  }
2295  }
2296  }
2297  }
2298  }
2299 
2301 }
2302 
2303 MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
2304  EntData &data) {
2306  if (type != MBVERTEX)
2308 
2309  auto create_tag = [this](const std::string tag_name, const int size) {
2310  double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2311  Tag th;
2312  CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
2313  th, MB_TAG_CREAT | MB_TAG_SPARSE,
2314  def_VAL);
2315  return th;
2316  };
2317 
2318  Tag th_w = create_tag("SpatialDisplacement", 3);
2319  Tag th_omega = create_tag("Omega", 3);
2320  Tag th_approxP = create_tag("Piola1Stress", 9);
2321  Tag th_sigma = create_tag("CauchyStress", 9);
2322  Tag th_log_u = create_tag("LogSpatialStreach", 9);
2323  Tag th_u = create_tag("SpatialStreach", 9);
2324  Tag th_h = create_tag("h", 9);
2325  Tag th_X = create_tag("X", 3);
2326  Tag th_detF = create_tag("detF", 1);
2327  Tag th_angular_momentum = create_tag("AngularMomentum", 3);
2328 
2329  Tag th_u_eig_vec = create_tag("SpatialStreachEigenVec", 9);
2330  Tag th_u_eig_vals = create_tag("SpatialStreachEigenVals", 3);
2331 
2332  int nb_gauss_pts = data.getN().size1();
2333  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
2334  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2335  "Nb. of integration points is not equal to number points on "
2336  "post-processing mesh");
2337  }
2338 
2339  auto t_w = getFTensor1FromMat<3>(dataAtPts->wAtPts);
2340  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
2341  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2342  auto t_log_u =
2343  getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachTensorAtPts);
2344  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
2345  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
2346  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2347  auto t_coords = getFTensor1CoordsAtGaussPts();
2348 
2353 
2354  // vectors
2355  VectorDouble3 v(3);
2356  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
2357  auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
2359  t_v(i) = t_d(i);
2360  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2361  &*v.data().begin());
2363  };
2364 
2365  MatrixDouble3by3 m(3, 3);
2367  &m(0, 0), &m(0, 1), &m(0, 2),
2368 
2369  &m(1, 0), &m(1, 1), &m(1, 2),
2370 
2371  &m(2, 0), &m(2, 1), &m(2, 2));
2372 
2373  auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
2375  t_m(i, j) = t_d(i, j);
2376  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2377  &*m.data().begin());
2379  };
2380 
2381  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2382 
2383  // vetors
2384  CHKERR save_vec_tag(th_w, t_w, gg);
2385  CHKERR save_vec_tag(th_X, t_coords, gg);
2386  CHKERR save_vec_tag(th_omega, t_omega, gg);
2387 
2388  // tensors
2389  CHKERR save_mat_tag(th_h, t_h, gg);
2390 
2391  FTensor::Tensor2<double, 3, 3> t_log_u_tmp;
2392  for (int ii = 0; ii != 3; ++ii)
2393  for (int jj = 0; jj != 3; ++jj)
2394  t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
2395 
2396  CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
2397 
2399  for (int ii = 0; ii != 3; ++ii)
2400  for (int jj = 0; jj != 3; ++jj)
2401  t_u_tmp(ii, jj) = t_u(ii, jj);
2402 
2403  CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2404  CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2405 
2406  const double jac = dEterminant(t_h);
2408  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
2409  CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2410  CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
2411 
2413  t_PhT(i, k) = t_approx_P(i, j) * t_R(k, j);
2414  FTensor::Tensor1<double, 3> t_leviPRT;
2415  t_leviPRT(k) = levi_civita(i, l, k) * t_PhT(i, l);
2416 
2417  CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
2418  &t_leviPRT(0));
2419 
2420  auto get_eiegn_vector_symmetric = [&](auto &t_u) {
2422 
2423  for (int ii = 0; ii != 3; ++ii)
2424  for (int jj = 0; jj != 3; ++jj)
2425  t_m(ii, jj) = t_u(ii, jj);
2426 
2427  VectorDouble3 eigen_values(3);
2428 
2429  // LAPACK - eigenvalues and vectors. Applied twice for initial creates
2430  // memory space
2431  int n = 3, lda = 3, info, lwork = -1;
2432  double wkopt;
2433  info = lapack_dsyev('V', 'U', n, &(m.data()[0]), lda,
2434  &(eigen_values.data()[0]), &wkopt, lwork);
2435  if (info != 0)
2436  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2437  "is something wrong with lapack_dsyev info = %d", info);
2438  lwork = (int)wkopt;
2439  double work[lwork];
2440  info = lapack_dsyev('V', 'U', n, &(m.data()[0]), lda,
2441  &(eigen_values.data()[0]), work, lwork);
2442  if (info != 0)
2443  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2444  "is something wrong with lapack_dsyev info = %d", info);
2445 
2446  CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
2447  &*m.data().begin());
2448  CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
2449  &*eigen_values.data().begin());
2450 
2452  };
2453 
2454  CHKERR get_eiegn_vector_symmetric(t_u);
2455 
2456  ++t_w;
2457  ++t_h;
2458  ++t_log_u;
2459  ++t_u;
2460  ++t_omega;
2461  ++t_R;
2462  ++t_approx_P;
2463  ++t_coords;
2464  }
2465 
2467 }
2468 
2470 OpCalculateStrainEnergy::doWork(int side, EntityType type,
2473  int nb_integration_pts = data.getN().size1();
2474  auto v = getVolume();
2475  auto t_w = getFTensor0IntegrationWeight();
2476  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2477  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2478 
2481 
2482  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2483  const double a = t_w * v;
2484  (*energy) += a * t_P(i, J) * t_h(i, J);
2485  ++t_w;
2486  ++t_P;
2487  ++t_h;
2488  }
2489 
2491 }
2492 
2493 MoFEMErrorCode OpL2Transform::doWork(int side, EntityType type,
2496  if (type == MBTET) {
2497  auto v = getVolume();
2498  for (auto &base : data.getN(AINSWORTH_LEGENDRE_BASE).data())
2499  base /= v;
2500  }
2502 }
2503 
2504 
2505 } // namespace EshelbianPlasticity
static MoFEMErrorCode matrix_exponent_derivative(MatrixDouble &log_u, MatrixDouble &diff_u, std::array< MatrixDouble, 6 > &diff2_u)
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
FTensor::Index< 'l', 3 > l
static MoFEMErrorCode matrix_rotation_derivative(MatrixDouble &omega, std::array< MatrixDouble, 3 > &diff2_omega)
static FTensor::Tensor3< typename TensorTypeExtractor< T >::Type, 3, 3, 3 > get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
multi_index_container< BlockMatData, indexed_by< ordered_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType >, member< BlockMatData, int, &BlockMatData::rowSide >, member< BlockMatData, int, &BlockMatData::colSide > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField > > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::rowField > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::colField > > > > BlockMatContainor
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
Kronecker Delta class.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
JSON compatible output.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
static MoFEMErrorCode tensor_exponent(FTensor::Tensor2_symmetric< T, 3 > &t_log_s_u, FTensor::Tensor2_symmetric< T, 3 > &t_s_u, FTensor::Tensor4< T, 3, 3, 3, 3 > &t_diff_s_u)
static Index< 'n', 3 > n
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:273
const VectorInt & getIndices() const
Get global indices of dofs on entity.
static Index< 'm', 3 > m
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
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:223
double tol
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
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
static Number< 2 > N2
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:684
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
static const bool debug
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
FTensor::Index< 'j', 3 > j
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:102
#define CHKERR
Inline error check.
Definition: definitions.h:604
const double r
rate factor
static Number< 1 > N1
constexpr double omega
static FTensor::Tensor4< typename TensorTypeExtractor< T >::Type, 3, 3, 3, 3 > calculate_diff_viscous_stress(FTensor::Tensor2< T, 3, 3 > &t, FTensor::Tensor4< T, 3, 3, 3, 3 > &t_diff)
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Eshelbian plasticity interface.
Data on single entity (This is passed as argument to DataOperator::doWork)
static Index< 'K', 3 > K
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
static const double eps
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
static Number< 0 > N0
const int N
Definition: speed_test.cpp:3
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > calculate_viscous_stress(FTensor::Tensor2< T, 3, 3 > &t)
static __CLPK_integer lapack_dsysv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:232
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:188
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
static Index< 'J', 3 > J
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