v0.9.1
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  t_R(i, j) = 0.;
36  t_R(0, 0) = 1.;
37  t_R(1, 1) = 1.;
38  t_R(2, 2) = 1.;
39 
40  const typename TensorTypeExtractor<T>::Type angle =
41  sqrt(t_omega(i) * t_omega(i));
42 
43  constexpr double tol = 1e-18;
44  if (std::abs(angle) < tol) {
45  return t_R;
46  }
47 
49  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
50  const typename TensorTypeExtractor<T>::Type a = sin(angle) / angle;
51  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
52  const typename TensorTypeExtractor<T>::Type b =
53  2. * ss_2 * ss_2 / (angle * angle);
54  t_R(i, j) += a * t_Omega(i, j);
55  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
56 
57  return t_R;
58 }
59 
60 template <typename T>
63  &t_omega) {
64 
69 
70  const typename TensorTypeExtractor<T>::Type angle =
71  sqrt(t_omega(i) * t_omega(i));
72 
73  constexpr double tol = 1e-18;
74  if (std::abs(angle) < tol) {
76  auto t_R = get_rotation_form_vector(t_omega);
77  t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, l, k) * t_R(l, j);
78  return t_diff_R;
79  }
80 
83  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
84 
85  const typename TensorTypeExtractor<T>::Type angle2 = angle * angle;
86  const typename TensorTypeExtractor<T>::Type ss = sin(angle);
87  const typename TensorTypeExtractor<T>::Type a = ss / angle;
88  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
89  const typename TensorTypeExtractor<T>::Type b = 2. * ss_2 * ss_2 / angle2;
90  t_diff_R(i, j, k) = 0;
91  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
92  t_diff_R(i, j, k) +=
93  b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
94  FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
95 
96  const typename TensorTypeExtractor<T>::Type cc = cos(angle);
97  const typename TensorTypeExtractor<T>::Type cc_2 = cos(angle / 2.);
98  const typename TensorTypeExtractor<T>::Type diff_a =
99  (angle * cc - ss) / angle2;
100 
101  const typename TensorTypeExtractor<T>::Type diff_b =
102  (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
103  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
104  t_diff_R(i, j, k) +=
105  diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
106 
107  return t_diff_R;
108 }
109 
110 template <typename T>
111 static MoFEMErrorCode
114  FTensor::Tensor4<T, 3, 3, 3, 3> &t_diff_s_u) {
116 
117  typedef typename TensorTypeExtractor<T>::Type A;
118 
119  auto calc_exponent = [](FTensor::Tensor2<A, 3, 3> &t_log_u,
123 
124  constexpr int N = 800;
125  constexpr double eps = 1e-12;
126 
132 
137 
140 
141  t_one2(i, j) = 0;
142  for (auto ii : {0, 1, 2})
143  t_one2(ii, ii) = 1;
144  t_one4(i, j, k, l) = t_one2(j, l) * t_one2(i, k);
145 
146  t_p(i, j) = t_log_u(i, j);
147  t_diff_p(i, j, k, l) = t_one4(i, j, k, l);
148 
149  t_u(i, j) = t_p(i, j);
150  for (auto ii : {0, 1, 2})
151  t_u(ii, ii) += 1;
152  t_diff_u(i, j, k, l) = t_diff_p(i, j, k, l);
153 
154  for (int nn = 2; nn <= N; ++nn) {
155 
156  t_diff_tmp(i, j, k, l) =
157  (t_p(i, m) / static_cast<double>(nn)) * t_one4(m, j, k, l) +
158  t_diff_p(i, m, k, l) * (t_log_u(m, j) / static_cast<double>(nn));
159 
160  t_tmp(i, j) = t_p(i, k) * (t_log_u(k, j) / static_cast<double>(nn));
161 
162  t_p(i, j) = t_tmp(i, j);
163  t_diff_p(i, j, k, l) = t_diff_tmp(i, j, k, l);
164 
165  t_u(i, j) += t_p(i, j);
166  t_diff_u(i, j, k, l) += t_diff_p(i, j, k, l);
167 
168  double e = std::abs(sqrt(t_p(i, j) * t_p(i, j)));
169  if (e < eps)
171  }
173  };
174 
177 
179  for (auto ii : {0, 1, 2})
180  for (auto jj : {0, 1, 2})
181  t_log_u(ii, jj) = t_log_s_u(ii, jj);
182 
183  CHKERR calc_exponent(t_log_u, t_u, t_diff_u);
184 
185  for (auto ii : {0, 1, 2})
186  for (auto jj = 0; jj != 3; ++jj) {
187  t_s_u(ii, jj) = (t_u(ii, jj) + t_u(jj, ii)) / 2.;
188  for (auto kk : {0, 1, 2})
189  for (auto ll = 0; ll <= kk; ++ll) {
190 
191  if (ll != kk)
192  t_diff_s_u(ii, jj, kk, ll) =
193 
194  t_diff_u(ii, jj, kk, ll) + t_diff_u(ii, jj, ll, kk);
195 
196  else
197  t_diff_s_u(ii, jj, kk, ll) = t_diff_u(ii, jj, kk, ll);
198  }
199  }
200 
202 }
203 
204 static MoFEMErrorCode
206  std::array<MatrixDouble, 6> &diff2_u) {
212 
213  for (auto ll : {0, 1, 2, 3, 4, 5})
214  diff2_u[ll].resize(81, log_u.size2(), false);
215 
216  std::array<FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>, 6>
217  t_diff2_u = {
218 
220  &diff2_u[0](0, 0), diff2_u[0].size2()),
222  &diff2_u[1](0, 0), diff2_u[1].size2()),
224  &diff2_u[2](0, 0), diff2_u[2].size2()),
226  &diff2_u[3](0, 0), diff2_u[3].size2()),
228  &diff2_u[4](0, 0), diff2_u[4].size2()),
230  &diff2_u[5](0, 0), diff2_u[5].size2())
231 
232  };
233 
234  constexpr double h = 1e-18;
235  const std::complex<double> ih = h * 1i;
236 
237  int ss = 0;
238  for (auto ii : {0, 1, 2})
239  for (auto jj = 0; jj <= ii; ++jj, ++ss) {
240 
241  auto t_log_u = getFTensor2SymmetricFromMat<3>(log_u);
242  auto t_diff_u =
244  &diff_u(0, 0), diff_u.size2());
245 
247  tc_h(i, j) = 0;
248  tc_h(ii, jj) = ih;
249 
250  for (int gg = 0; gg != log_u.size2(); ++gg) {
251 
254  FTensor::Tensor4<std::complex<double>, 3, 3, 3, 3> tc_diff_u;
255 
256  tc_log_u(i, j) = t_log_u(i, j) + tc_h(i, j);
257  CHKERR tensor_exponent(tc_log_u, tc_u, tc_diff_u);
258 
259  (t_diff2_u[ss])(i, j, k, l) = 0;
260  for (auto kk : {0, 1, 2})
261  for (auto ll : {0, 1, 2})
262  for (auto mm : {0, 1, 2})
263  for (auto nn : {0, 1, 2}) {
264  (t_diff2_u[ss])(kk, ll, mm, nn) =
265  std::imag(tc_diff_u(kk, ll, mm, nn)) / h;
266  }
267 
268  // for (auto kk : {0, 1, 2})
269  // for (auto ll = 0; ll <= 2; ++ll) {
270 
271  // constexpr double eps = 1e-6;
272  // if (std::abs(t_diff_u(kk, ll, ii, jj) -
273  // std::imag(tc_u(kk, ll)) / h) > eps) {
274  // std::ostringstream ss;
275  // ss << "Error " << kk << " " << ll << " : "
276  // << t_diff_u(kk, ll, ii, jj) << " "
277  // << std::imag(tc_u(kk, ll)) / h << " error "
278  // << t_diff_u(kk, ll, ii, jj) - std::imag(tc_u(kk, ll)) / h
279  // << " "
280  // << t_diff_u(kk, ll, ii, jj) / (std::imag(tc_u(kk, ll)) / h)
281  // << endl;
282  // SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "%s",
283  // ss.str().c_str());
284  // }
285  // }
286 
287  ++t_log_u;
288  ++t_diff_u;
289  ++(t_diff2_u[ss]);
290  }
291  }
292 
294 }
295 
296 static MoFEMErrorCode
298  std::array<MatrixDouble, 3> &diff2_omega) {
303 
304  for (auto ll : {0, 1, 2})
305  diff2_omega[ll].resize(27, omega.size2(), false);
306 
307  std::array<FTensor::Tensor3<FTensor::PackPtr<double *, 1>, 3, 3, 3>, 3>
308  t_diff2_omega = {
309 
310  getFTensor3FromMat(diff2_omega[0]),
311 
312  getFTensor3FromMat(diff2_omega[1]),
313 
314  getFTensor3FromMat(diff2_omega[2])
315 
316  };
317 
318  constexpr double h = 1e-18;
319  const std::complex<double> ih = h * 1i;
320 
321  for (auto ii : {0, 1, 2}) {
322 
323  auto t_omega = getFTensor1FromMat<3>(omega);
324 
326  tc_h(i) = 0;
327  tc_h(ii) = ih;
328 
329  for (int gg = 0; gg != omega.size2(); ++gg) {
330 
332  tc_omega(i) = t_omega(i) + tc_h(i);
333 
334  auto tc_diff_R = get_diff_rotation_form_vector(tc_omega);
335 
336  (t_diff2_omega[ii])(i, j, k) = 0;
337  for (auto kk : {0, 1, 2})
338  for (auto ll : {0, 1, 2})
339  for (auto mm : {0, 1, 2})
340  (t_diff2_omega[ii])(kk, ll, mm) =
341  std::imag(tc_diff_R(kk, ll, mm)) / h;
342 
343  ++t_omega;
344  ++(t_diff2_omega[ii]);
345  }
346  }
348 }
349 
350 template <typename T>
353 
357 
358  auto get_dot_fun = [i, j, k](auto &t) {
360  t_c(i, j) = t(i, j);
361  return t_c;
362  };
363 
364  auto get_trace = [](auto &t) {
365  typename TensorTypeExtractor<T>::Type tr = 0;
366  for (auto ii = 0; ii != 3; ++ii)
367  tr += t(ii, ii) / 3.;
368  return tr;
369  };
370 
371  auto get_deviator = [i, j, get_trace](auto &&t) {
372  const auto tr = get_trace(t);
374  t_dev(i, j) = t(i, j);
375  for (auto ii = 0; ii != 3; ++ii)
376  t_dev(ii, ii) -= tr;
377  return t_dev;
378  };
379 
380  return get_deviator(get_dot_fun(t));
381 }
382 
383 template <typename T>
387 
393 
394  auto get_dot_fun = [i, j, k, m, n](auto &t, auto &t_diff) {
396  t_diff_c;
397  t_diff_c(i, j, m, n) = t_diff(i, j, m, n);
398  return t_diff_c;
399  };
400 
401  auto get_trace = [i, j](auto &t_diff) {
403  t_diff_tr(i, j) = 0;
404  for (auto ii = 0; ii != 3; ++ii)
405  for (auto mm = 0; mm != 3; ++mm)
406  for (auto nn = 0; nn != 3; ++nn)
407  t_diff_tr(mm, nn) += t_diff(ii, ii, mm, nn) / 3.;
408  return t_diff_tr;
409  };
410 
411  auto get_deviator = [i, j, m, n, get_trace](auto &&t_diff) {
412  auto t_diff_tr = get_trace(t_diff);
414  t_diff_dev;
415  t_diff_dev(i, j, m, n) = t_diff(i, j, m, n);
416  for (auto mm = 0; mm != 3; ++mm)
417  for (auto nn = 0; nn != 3; ++nn)
418  for (auto ii = 0; ii != 3; ++ii)
419  t_diff_dev(ii, ii, mm, nn) -= t_diff_tr(mm, nn);
420  return t_diff_dev;
421  };
422 
423  return get_deviator(get_dot_fun(t, t_diff));
424 }
425 
426 MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
427  EntityType type,
428  EntData &data) {
430  if (type != MBTET)
432  int nb_integration_pts = data.getN().size1();
437 
438  dataAtPts->hAtPts->resize(9, nb_integration_pts, false);
439  dataAtPts->rotMatAtPts->resize(9, nb_integration_pts, false);
440  dataAtPts->diffRotMatAtPts->resize(27, nb_integration_pts, false);
441  dataAtPts->streachTensorAtPts->resize(6, nb_integration_pts, false);
442  dataAtPts->diffStreachTensorAtPts->resize(81, nb_integration_pts, false);
443 
444  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
445  auto t_omega = getFTensor1FromMat<3>(*(dataAtPts->rotAxisAtPts));
446  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
447  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
448  auto t_log_u =
449  getFTensor2SymmetricFromMat<3>(*(dataAtPts->logStreachTensorAtPts));
450  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
451  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
452  &(*dataAtPts->diffStreachTensorAtPts)(0, 0), nb_integration_pts);
453 
454  for (int gg = 0; gg != nb_integration_pts; ++gg) {
455 
456  auto t0_diff = get_diff_rotation_form_vector(t_omega);
457  auto t0 = get_rotation_form_vector(t_omega);
458 
459  t_diff_R(i, j, k) = t0_diff(i, j, k);
460  t_R(i, j) = t0(i, j);
461 
462  CHKERR tensor_exponent(t_log_u, t_u, t_diff_u);
463 
464  t_h(i, j) = t_R(i, k) * t_u(k, j);
465 
466  ++t_h;
467  ++t_R;
468  ++t_diff_R;
469  ++t_log_u;
470  ++t_u;
471  ++t_diff_u;
472  ++t_omega;
473  }
474 
475  CHKERR matrix_exponent_derivative(*(dataAtPts->logStreachTensorAtPts),
476  *(dataAtPts->diffStreachTensorAtPts),
477  dataAtPts->expLogUHessian);
478 
479  // CHKERR matrix_rotation_derivative(*(dataAtPts->rotAxisAtPts),
480  // dataAtPts->rotationHessian);
481 
483 }
484 
485 MoFEMErrorCode OpSpatialEquilibrium::integrate(EntData &data) {
487  int nb_dofs = data.getIndices().size();
488  int nb_integration_pts = data.getN().size1();
489  auto v = getVolume();
490  auto t_w = getFTensor0IntegrationWeight();
491  auto t_div_P = getFTensor1FromMat<3>(*(dataAtPts->divPAtPts));
492  auto t_s_dot_w = getFTensor1FromMat<3>(*(dataAtPts->wDotAtPts));
493  int nb_base_functions = data.getN().size2();
494  auto t_row_base_fun = data.getFTensor0N();
495 
497  auto get_ftensor1 = [](auto &v) {
498  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
499  &v[2]);
500  };
501 
502  for (int gg = 0; gg != nb_integration_pts; ++gg) {
503  double a = v * t_w;
504  auto t_nf = get_ftensor1(nF);
505  int bb = 0;
506  for (; bb != nb_dofs / 3; ++bb) {
507  t_nf(i) += a * t_row_base_fun * t_div_P(i);
508  t_nf(i) -= a * t_row_base_fun * alpha_w * t_s_dot_w(i);
509  ++t_nf;
510  ++t_row_base_fun;
511  }
512  for (; bb != nb_base_functions; ++bb)
513  ++t_row_base_fun;
514  ++t_w;
515  ++t_div_P;
516  ++t_s_dot_w;
517  }
518 
520 }
521 
522 MoFEMErrorCode OpSpatialRotation::integrate(EntData &data) {
524  int nb_dofs = data.getIndices().size();
525  int nb_integration_pts = data.getN().size1();
526  auto v = getVolume();
527  auto t_w = getFTensor0IntegrationWeight();
528  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
529  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
530 
531  int nb_base_functions = data.getN().size2();
532  auto t_row_base_fun = data.getFTensor0N();
538  auto get_ftensor1 = [](auto &v) {
539  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
540  &v[2]);
541  };
542 
543  for (int gg = 0; gg != nb_integration_pts; ++gg) {
544  double a = v * t_w;
545  auto t_nf = get_ftensor1(nF);
547  t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
548  FTensor::Tensor1<double, 3> t_leviPRT;
549  t_leviPRT(k) = levi_civita(i, m, k) * t_PRT(i, m);
550  int bb = 0;
551  for (; bb != nb_dofs / 3; ++bb) {
552  t_nf(k) += a * t_row_base_fun * t_leviPRT(k);
553  ++t_nf;
554  ++t_row_base_fun;
555  }
556  for (; bb != nb_base_functions; ++bb)
557  ++t_row_base_fun;
558  ++t_w;
559  ++t_approx_P;
560  ++t_R;
561  }
563 }
564 
565 MoFEMErrorCode OpSpatialPhysical::integrate(EntData &data) {
567  int nb_dofs = data.getIndices().size();
568  int nb_integration_pts = data.getN().size1();
569  auto v = getVolume();
570  auto t_w = getFTensor0IntegrationWeight();
571  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
572  auto t_P = getFTensor2FromMat<3, 3>(*(dataAtPts->PAtPts));
573  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
574  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
575  auto t_dot_log_u =
576  getFTensor2SymmetricFromMat<3>(*(dataAtPts->logStreachDotTensorAtPts));
577  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
578  &(*dataAtPts->diffStreachTensorAtPts)(0, 0),
579  (*dataAtPts->diffStreachTensorAtPts).size2());
580 
585  auto get_ftensor2 = [](auto &v) {
587  &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
588  };
589 
590  int nb_base_functions = data.getN().size2();
591  auto t_row_base_fun = data.getFTensor0N();
592  for (int gg = 0; gg != nb_integration_pts; ++gg) {
593  double a = v * t_w;
594  auto t_nf = get_ftensor2(nF);
595 
597  t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
599  t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
600 
602  t_dot_u(i, j) = t_u(i, k) * t_dot_log_u(k, j);
603 
604  auto t_viscous_stress = calculate_dot_deviator(t_dot_u);
605  // t_viscous_stress(i, j) = t_dot_u(i, j);
606  t_viscous_stress(i, j) *= alpha_u;
607 
609  t1(i, j) = 0;
610  for (int ii : {0, 1, 2})
611  for (int jj = 0; jj <= ii; ++jj)
612  for (int mm : {0, 1, 2})
613  for (int nn : {0, 1, 2})
614  t1(ii, jj) += a * t_diff_u(mm, nn, ii, jj) *
615  (t_deltaP(mm, nn) - t_viscous_stress(mm, nn));
616 
617  int bb = 0;
618  for (; bb != nb_dofs / 6; ++bb) {
619 
620  for (int ii : {0, 1, 2})
621  for (int jj = 0; jj <= ii; ++jj) {
622  t_nf(ii, jj) += t_row_base_fun * t1(ii, jj);
623  }
624 
625  ++t_nf;
626  ++t_row_base_fun;
627  }
628  for (; bb != nb_base_functions; ++bb)
629  ++t_row_base_fun;
630 
631  ++t_w;
632  ++t_approx_P;
633  ++t_P;
634  ++t_R;
635  ++t_dot_log_u;
636  ++t_u;
637  ++t_diff_u;
638  }
640 }
641 
642 MoFEMErrorCode OpSpatialConsistencyP::integrate(EntData &data) {
644  int nb_dofs = data.getIndices().size();
645  int nb_integration_pts = data.getN().size1();
646  auto v = getVolume();
647  auto t_w = getFTensor0IntegrationWeight();
648  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
649 
650  int nb_base_functions = data.getN().size2() / 3;
651  auto t_row_base_fun = data.getFTensor1N<3>();
655  auto get_ftensor1 = [](auto &v) {
656  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
657  &v[2]);
658  };
659 
660  for (int gg = 0; gg != nb_integration_pts; ++gg) {
661  double a = v * t_w;
662  auto t_nf = get_ftensor1(nF);
663  int bb = 0;
664  for (; bb != nb_dofs / 3; ++bb) {
665  t_nf(i) += a * t_row_base_fun(j) * t_h(i, j);
666  t_nf(i) -= a * t_row_base_fun(i);
667  ++t_nf;
668  ++t_row_base_fun;
669  }
670  for (; bb != nb_base_functions; ++bb) {
671  ++t_row_base_fun;
672  }
673  ++t_w;
674  ++t_h;
675  }
676 
678 }
679 
680 MoFEMErrorCode OpSpatialConsistencyBubble::integrate(EntData &data) {
682  int nb_dofs = data.getIndices().size();
683  int nb_integration_pts = data.getN().size1();
684  auto v = getVolume();
685  auto t_w = getFTensor0IntegrationWeight();
686  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
687 
688  int nb_base_functions = data.getN().size2() / 9;
689  auto t_row_base_fun = data.getFTensor2N<3, 3>();
693  auto get_ftensor0 = [](auto &v) {
695  };
696 
697  for (int gg = 0; gg != nb_integration_pts; ++gg) {
698  double a = v * t_w;
699  auto t_nf = get_ftensor0(nF);
700  int bb = 0;
701  for (; bb != nb_dofs; ++bb) {
702  t_nf += a * t_row_base_fun(i, j) * t_h(i, j);
703  t_nf -= a * t_row_base_fun(i, i);
704  ++t_nf;
705  ++t_row_base_fun;
706  }
707  for (; bb != nb_base_functions; ++bb) {
708  ++t_row_base_fun;
709  }
710  ++t_w;
711  ++t_h;
712  }
713 
715 }
716 
717 MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
719  int nb_dofs = data.getIndices().size();
720  int nb_integration_pts = data.getN().size1();
721  auto v = getVolume();
722  auto t_w = getFTensor0IntegrationWeight();
723  auto t_s_w = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
724  int nb_base_functions = data.getN().size2() / 3;
725  auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
727  auto get_ftensor1 = [](auto &v) {
728  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
729  &v[2]);
730  };
731 
732  for (int gg = 0; gg != nb_integration_pts; ++gg) {
733  double a = v * t_w;
734  auto t_nf = get_ftensor1(nF);
735  int bb = 0;
736  for (; bb != nb_dofs / 3; ++bb) {
737  double div_row_base = t_row_diff_base_fun(i, i);
738  t_nf(i) += a * div_row_base * t_s_w(i);
739  ++t_nf;
740  ++t_row_diff_base_fun;
741  }
742  for (; bb != nb_base_functions; ++bb) {
743  ++t_row_diff_base_fun;
744  }
745  ++t_w;
746  ++t_s_w;
747  }
748 
750 }
751 
752 MoFEMErrorCode OpDispBc::integrate(EntData &data) {
754  // get entity of face
755  EntityHandle fe_ent = getFEEntityHandle();
756  // interate over all boundary data
757  for (auto &bc : (*bcDispPtr)) {
758  // check if finite element entity is part of boundary condition
759  if (bc.faces.find(fe_ent) != bc.faces.end()) {
760  int nb_dofs = data.getIndices().size();
761  int nb_integration_pts = data.getN().size1();
762  auto t_normal = getFTensor1Normal();
763  auto t_w = getFTensor0IntegrationWeight();
764  int nb_base_functions = data.getN().size2() / 3;
765  auto t_row_base_fun = data.getFTensor1N<3>();
766 
769  auto get_ftensor1 = [](auto &v) {
770  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
771  &v[2]);
772  };
773 
774  // get bc data
775  FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
776  t_bc_disp(i) *= getFEMethod()->ts_t;
777 
778  for (int gg = 0; gg != nb_integration_pts; ++gg) {
780  t_bc_res(i) = t_bc_disp(i);
781  auto t_nf = get_ftensor1(nF);
782  int bb = 0;
783  for (; bb != nb_dofs / 3; ++bb) {
784  t_nf(i) -=
785  t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
786  ++t_nf;
787  ++t_row_base_fun;
788  }
789  for (; bb != nb_base_functions; ++bb)
790  ++t_row_base_fun;
791 
792  ++t_w;
793  }
794  }
795  }
797 }
798 
799 MoFEMErrorCode OpRotationBc::integrate(EntData &data) {
801  // get entity of face
802  EntityHandle fe_ent = getFEEntityHandle();
803  // interate over all boundary data
804  for (auto &bc : (*bcRotPtr)) {
805  // check if finite element entity is part of boundary condition
806  if (bc.faces.find(fe_ent) != bc.faces.end()) {
807  int nb_dofs = data.getIndices().size();
808  int nb_integration_pts = data.getN().size1();
809  auto t_normal = getFTensor1Normal();
810  auto t_w = getFTensor0IntegrationWeight();
811 
812  int nb_base_functions = data.getN().size2() / 3;
813  auto t_row_base_fun = data.getFTensor1N<3>();
817  auto get_ftensor1 = [](auto &v) {
818  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
819  &v[2]);
820  };
821 
822  // get bc data
823  FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
824 
825  double theta = bc.theta;
826  theta *= getFEMethod()->ts_t;
827 
829  const double a = sqrt(t_normal(i) * t_normal(i));
830  t_omega(i) = t_normal(i) * (theta / a);
831  auto t_R = get_rotation_form_vector(t_omega);
832  auto t_coords = getFTensor1CoordsAtGaussPts();
833 
834  for (int gg = 0; gg != nb_integration_pts; ++gg) {
836  t_delta(i) = t_center(i) - t_coords(i);
838  t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
839 
840  auto t_nf = get_ftensor1(nF);
841  int bb = 0;
842  for (; bb != nb_dofs / 3; ++bb) {
843  t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
844  ++t_nf;
845  ++t_row_base_fun;
846  }
847  for (; bb != nb_base_functions; ++bb)
848  ++t_row_base_fun;
849 
850  ++t_w;
851  ++t_coords;
852  }
853  }
854  }
856 }
857 
858 MoFEMErrorCode FeTractionBc::preProcess() {
860 
861  struct FaceRule {
862  int operator()(int p_row, int p_col, int p_data) const {
863  return 2 * (p_data + 1);
864  }
865  };
866 
867  if (ts_ctx == CTX_TSSETIFUNCTION) {
868 
869  // Loop boundary elements with traction boundary conditions
871  fe.getOpPtrVector().push_back(
872  new OpTractionBc(std::string("P") /* + "_RT"*/, *this));
873  fe.getRuleHook = FaceRule();
874  fe.ts_t = ts_t;
875  CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
876  fe);
877  }
878 
880 }
881 
882 MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type, EntData &data) {
884 
887 
888  auto t_normal = getFTensor1Normal();
889  const double nrm2 = sqrt(t_normal(i) * t_normal(i));
890  FTensor::Tensor1<double, 3> t_unit_normal;
891  t_unit_normal(i) = t_normal(i) / nrm2;
892  int nb_dofs = data.getFieldData().size();
893  int nb_integration_pts = data.getN().size1();
894  int nb_base_functions = data.getN().size2() / 3;
895  double ts_t = getFEMethod()->ts_t;
896 
897  auto integrate_matrix = [&]() {
899 
900  auto t_w = getFTensor0IntegrationWeight();
901  matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
902  matPP.clear();
903 
904  auto t_row_base_fun = data.getFTensor1N<3>();
905  for (int gg = 0; gg != nb_integration_pts; ++gg) {
906 
907  int rr = 0;
908  for (; rr != nb_dofs / 3; ++rr) {
909  const double a = t_row_base_fun(i) * t_unit_normal(i);
910  auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
911  for (int cc = 0; cc != nb_dofs / 3; ++cc) {
912  const double b = t_col_base_fun(i) * t_unit_normal(i);
913  matPP(rr, cc) += t_w * a * b;
914  ++t_col_base_fun;
915  }
916  ++t_row_base_fun;
917  }
918 
919  for (; rr != nb_base_functions; ++rr)
920  ++t_row_base_fun;
921 
922  ++t_w;
923  }
924 
926  };
927 
928  auto integrate_rhs = [&](auto &bc) {
930 
931  auto t_w = getFTensor0IntegrationWeight();
932  vecPv.resize(3, nb_dofs / 3, false);
933  vecPv.clear();
934 
935  auto t_row_base_fun = data.getFTensor1N<3>();
936  double ts_t = getFEMethod()->ts_t;
937 
938  for (int gg = 0; gg != nb_integration_pts; ++gg) {
939  int rr = 0;
940  for (; rr != nb_dofs / 3; ++rr) {
941  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
942  for (int dd = 0; dd != 3; ++dd)
943  if (bc.flags[dd])
944  vecPv(dd, rr) += t * bc.vals[dd];
945  ++t_row_base_fun;
946  }
947  for (; rr != nb_base_functions; ++rr)
948  ++t_row_base_fun;
949  ++t_w;
950  }
952  };
953 
954  auto integrate_rhs_cook = [&](auto &bc) {
956 
957  vecPv.resize(3, nb_dofs / 3, false);
958  vecPv.clear();
959 
960  auto t_w = getFTensor0IntegrationWeight();
961  auto t_row_base_fun = data.getFTensor1N<3>();
962  auto t_coords = getFTensor1CoordsAtGaussPts();
963 
964  for (int gg = 0; gg != nb_integration_pts; ++gg) {
965 
966  auto calc_tau = [](double y) {
967  y -= 44;
968  y /= (60 - 44);
969  return -y * (y - 1) / 0.25;
970  };
971 
972  const double tau = calc_tau(t_coords(1));
973 
974  int rr = 0;
975  for (; rr != nb_dofs / 3; ++rr) {
976  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
977  for (int dd = 0; dd != 3; ++dd)
978  if (bc.flags[dd])
979  vecPv(dd, rr) += t * tau * bc.vals[dd];
980  ++t_row_base_fun;
981  }
982 
983  for (; rr != nb_base_functions; ++rr)
984  ++t_row_base_fun;
985  ++t_w;
986  ++t_coords;
987  }
989  };
990 
991  // get entity of face
992  EntityHandle fe_ent = getFEEntityHandle();
993  for (auto &bc : *(bcFe.bcData)) {
994  if (bc.faces.find(fe_ent) != bc.faces.end()) {
995 
996  int nb_dofs = data.getFieldData().size();
997  if (nb_dofs) {
998 
999  CHKERR integrate_matrix();
1000  if (bc.blockName.compare(20, 4, "COOK") == 0)
1001  CHKERR integrate_rhs_cook(bc);
1002  else
1003  CHKERR integrate_rhs(bc);
1004 
1005  const auto info =
1006  lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
1007  nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
1008  if (info > 0)
1009  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1010  "The leading minor of order %i is "
1011  "not positive; definite;\nthe "
1012  "solution could not be computed",
1013  info);
1014 
1015  for (int dd = 0; dd != 3; ++dd)
1016  if (bc.flags[dd])
1017  for (int rr = 0; rr != nb_dofs / 3; ++rr)
1018  data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
1019 
1020  // int nb_integration_pts = data.getN().size1();
1021  // auto t_row_base_fun = data.getFTensor1N<3>();
1022  // auto t_coords = getFTensor1CoordsAtGaussPts();
1023  // for (int gg = 0; gg != nb_integration_pts; ++gg) {
1024 
1025  // FTensor::Tensor1<double, 3> t_v;
1026  // t_v(i) = 0;
1027 
1028  // int rr = 0;
1029  // for (; rr != nb_dofs / 3; ++rr) {
1030  // for (int dd = 0; dd != 3; ++dd) {
1031  // const double v = data.getFieldDofs()[3 * rr +
1032  // dd]->getFieldData(); t_v(dd) += t_row_base_fun(i) *
1033  // t_unit_normal(i) * v;
1034  // }
1035  // ++t_row_base_fun;
1036  // }
1037  // for (; rr != nb_base_functions; ++rr)
1038  // ++t_row_base_fun;
1039  // ++t_coords;
1040 
1041  // cerr << t_v << endl;
1042  // }
1043  }
1044  }
1045  }
1046 
1048 }
1049 
1050 MoFEMErrorCode OpSpatialEquilibrium_dw_dP::integrate(EntData &row_data,
1051  EntData &col_data) {
1053  int nb_integration_pts = row_data.getN().size1();
1054  int row_nb_dofs = row_data.getIndices().size();
1055  int col_nb_dofs = col_data.getIndices().size();
1056  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1058  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
1059  };
1061  auto v = getVolume();
1062  auto t_w = getFTensor0IntegrationWeight();
1063  int row_nb_base_functions = row_data.getN().size2();
1064  auto t_row_base_fun = row_data.getFTensor0N();
1065  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1066  double a = v * t_w;
1067  int rr = 0;
1068  for (; rr != row_nb_dofs / 3; ++rr) {
1069  auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
1070  auto t_m = get_ftensor1(K, 3 * rr, 0);
1071  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1072  double div_col_base = t_col_diff_base_fun(i, i);
1073  t_m(i) += a * t_row_base_fun * div_col_base;
1074  ++t_m;
1075  ++t_col_diff_base_fun;
1076  }
1077  ++t_row_base_fun;
1078  }
1079  for (; rr != row_nb_base_functions; ++rr)
1080  ++t_row_base_fun;
1081  ++t_w;
1082  }
1084 }
1085 
1086 MoFEMErrorCode OpSpatialEquilibrium_dw_dw::integrate(EntData &row_data,
1087  EntData &col_data) {
1089 
1090  if (alpha_w < std::numeric_limits<double>::epsilon())
1092 
1093  const int nb_integration_pts = row_data.getN().size1();
1094  const int row_nb_dofs = row_data.getIndices().size();
1095  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1097  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
1098 
1099  );
1100  };
1102 
1103  auto v = getVolume();
1104  auto t_w = getFTensor0IntegrationWeight();
1105 
1106  int row_nb_base_functions = row_data.getN().size2();
1107  auto t_row_base_fun = row_data.getFTensor0N();
1108 
1109  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1110  double a = v * t_w * alpha_w * getTSa();
1111 
1112  int rr = 0;
1113  for (; rr != row_nb_dofs / 3; ++rr) {
1114 
1115  auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
1116  auto t_m = get_ftensor1(K, 3 * rr, 0);
1117  for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1118  const double b = a * t_row_base_fun * t_col_base_fun;
1119  t_m(i) -= b;
1120  ++t_m;
1121  ++t_col_base_fun;
1122  }
1123 
1124  ++t_row_base_fun;
1125  }
1126 
1127  for (; rr != row_nb_base_functions; ++rr)
1128  ++t_row_base_fun;
1129 
1130  ++t_w;
1131  }
1132 
1134 }
1135 
1136 MoFEMErrorCode OpSpatialPhysical_du_dP::integrate(EntData &row_data,
1137  EntData &col_data) {
1139  int nb_integration_pts = row_data.getN().size1();
1140  int row_nb_dofs = row_data.getIndices().size();
1141  int col_nb_dofs = col_data.getIndices().size();
1142  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1144 
1145  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1146 
1147  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1148 
1149  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1150 
1151  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1152 
1153  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1154 
1155  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
1156  };
1157 
1163 
1164  auto v = getVolume();
1165  auto t_w = getFTensor0IntegrationWeight();
1166  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1167  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
1168  &(*dataAtPts->diffStreachTensorAtPts)(0, 0), nb_integration_pts);
1169 
1170  int row_nb_base_functions = row_data.getN().size2();
1171  auto t_row_base_fun = row_data.getFTensor0N();
1172  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1173  double a = v * t_w;
1174 
1175  int rr = 0;
1176  for (; rr != row_nb_dofs / 6; ++rr) {
1177 
1178  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1179  auto t_m = get_ftensor3(K, 6 * rr, 0);
1180 
1181  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1182 
1184  t_RTP_dP(i, j, k) = t_R(k, i) * t_col_base_fun(j) * t_row_base_fun;
1185 
1186  for (int kk = 0; kk != 3; ++kk)
1187  for (int ii = 0; ii != 3; ++ii)
1188  for (int jj = 0; jj <= ii; ++jj) {
1189  for (int mm = 0; mm != 3; ++mm)
1190  for (int nn = 0; nn != 3; ++nn)
1191  t_m(ii, jj, kk) +=
1192  a * t_diff_u(mm, nn, ii, jj) * t_RTP_dP(mm, nn, kk);
1193  }
1194 
1195  ++t_col_base_fun;
1196  ++t_m;
1197  }
1198 
1199  ++t_row_base_fun;
1200  }
1201  for (; rr != row_nb_base_functions; ++rr)
1202  ++t_row_base_fun;
1203  ++t_w;
1204  ++t_R;
1205  ++t_diff_u;
1206  }
1208 }
1209 
1210 MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
1211  EntData &col_data) {
1213  int nb_integration_pts = row_data.getN().size1();
1214  int row_nb_dofs = row_data.getIndices().size();
1215  int col_nb_dofs = col_data.getIndices().size();
1216  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1218  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
1219  &m(r + 5, c));
1220  };
1221 
1225  auto v = getVolume();
1226  auto t_w = getFTensor0IntegrationWeight();
1227  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1228  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
1229  &(*dataAtPts->diffStreachTensorAtPts)(0, 0), nb_integration_pts);
1230 
1231  int row_nb_base_functions = row_data.getN().size2();
1232  auto t_row_base_fun = row_data.getFTensor0N();
1233  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1234  double a = v * t_w;
1235 
1236  int rr = 0;
1237  for (; rr != row_nb_dofs / 6; ++rr) {
1238 
1239  auto t_m = get_ftensor2(K, 6 * rr, 0);
1240 
1241  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1242  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1243 
1244  FTensor::Tensor2<double, 3, 3> t_RTP_dBubble;
1245  t_RTP_dBubble(i, j) =
1246  a * t_row_base_fun * t_R(k, i) * t_col_base_fun(k, j);
1247 
1248  for (int ii = 0; ii != 3; ++ii)
1249  for (int jj = 0; jj <= ii; ++jj)
1250  for (int mm = 0; mm != 3; ++mm)
1251  for (int nn = 0; nn != 3; ++nn)
1252  t_m(ii, jj) += t_diff_u(mm, nn, ii, jj) * t_RTP_dBubble(mm, nn);
1253 
1254  ++t_m;
1255  ++t_col_base_fun;
1256  }
1257 
1258  ++t_row_base_fun;
1259  }
1260  for (; rr != row_nb_base_functions; ++rr)
1261  ++t_row_base_fun;
1262  ++t_w;
1263  ++t_R;
1264  ++t_diff_u;
1265  }
1267 }
1268 
1269 MoFEMErrorCode OpSpatialPhysical_du_du::integrate(EntData &row_data,
1270  EntData &col_data) {
1272 
1273  int nb_integration_pts = row_data.getN().size1();
1274  int row_nb_dofs = row_data.getIndices().size();
1275  int col_nb_dofs = col_data.getIndices().size();
1276  auto get_ftensor4 = [](MatrixDouble &m, const int r, const int c) {
1278 
1279  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1280  &m(r + 0, c + 4), &m(r + 0, c + 5),
1281 
1282  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1283  &m(r + 1, c + 4), &m(r + 1, c + 5),
1284 
1285  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1286  &m(r + 2, c + 4), &m(r + 2, c + 5),
1287 
1288  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1289  &m(r + 3, c + 4), &m(r + 3, c + 5),
1290 
1291  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1292  &m(r + 4, c + 4), &m(r + 4, c + 5),
1293 
1294  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1295  &m(r + 5, c + 4), &m(r + 5, c + 5)
1296 
1297  );
1298  };
1305 
1306  FTensor::Number<0> N0;
1307  FTensor::Number<1> N1;
1308  FTensor::Number<2> N2;
1309 
1310  auto v = getVolume();
1311  auto t_w = getFTensor0IntegrationWeight();
1312  auto t_P_dh0 = getFTensor3FromMat(*(dataAtPts->P_dh0));
1313  auto t_P_dh1 = getFTensor3FromMat(*(dataAtPts->P_dh1));
1314  auto t_P_dh2 = getFTensor3FromMat(*(dataAtPts->P_dh2));
1315  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1316  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
1317  auto t_P = getFTensor2FromMat<3, 3>(*(dataAtPts->PAtPts));
1318  auto t_dot_log_u =
1319  getFTensor2SymmetricFromMat<3>(*(dataAtPts->logStreachDotTensorAtPts));
1320  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
1321  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
1322  &(*dataAtPts->diffStreachTensorAtPts)(0, 0),
1323  (*dataAtPts->diffStreachTensorAtPts).size2());
1324  std::array<FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>, 6>
1325  t_diff2_u = {
1326 
1328  &dataAtPts->expLogUHessian[0](0, 0),
1329  dataAtPts->expLogUHessian[0].size2()),
1331  &dataAtPts->expLogUHessian[1](0, 0),
1332  dataAtPts->expLogUHessian[1].size2()),
1334  &dataAtPts->expLogUHessian[2](0, 0),
1335  dataAtPts->expLogUHessian[2].size2()),
1337  &dataAtPts->expLogUHessian[3](0, 0),
1338  dataAtPts->expLogUHessian[3].size2()),
1340  &dataAtPts->expLogUHessian[4](0, 0),
1341  dataAtPts->expLogUHessian[4].size2()),
1343  &dataAtPts->expLogUHessian[5](0, 0),
1344  dataAtPts->expLogUHessian[5].size2())
1345 
1346  };
1347 
1350  t_one2(i, j) = 0;
1351  for (auto ii : {0, 1, 2})
1352  t_one2(ii, ii) = 1;
1353  t_one4(i, j, k, l) = t_one2(j, l) * t_one2(i, k);
1354 
1355  FTensor::Tensor4<double, 3, 3, 3, 3> t_one4_symmetric;
1356  t_one4_symmetric(i, j, k, l) = 0;
1357  for (auto ii : {0, 1, 2})
1358  for (auto jj : {0, 1, 2})
1359  for (auto kk : {0, 1, 2})
1360  for (auto ll : {0, 1, 2}) {
1361 
1362  if (ll != kk)
1363  t_one4_symmetric(ii, jj, kk, ll) =
1364  t_one4(ii, jj, kk, ll) + t_one4(ii, jj, ll, kk);
1365  else
1366  t_one4_symmetric(ii, jj, kk, ll) = t_one4(ii, jj, kk, ll);
1367  }
1368 
1369  int row_nb_base_functions = row_data.getN().size2();
1370  auto t_row_base_fun = row_data.getFTensor0N();
1371 
1373  t_one(i, j) = 0;
1374  for (auto ii : {0, 1, 2})
1375  t_one(ii, ii) = 1;
1376 
1377  const double ts_a = getTSa();
1378  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1379  double a = v * t_w;
1380 
1382  t_P_dh(i, j, N0, k) = t_P_dh0(i, j, k);
1383  t_P_dh(i, j, N1, k) = t_P_dh1(i, j, k);
1384  t_P_dh(i, j, N2, k) = t_P_dh2(i, j, k);
1385 
1387  t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
1389  t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
1390 
1392  t_dot_u(i, j) = t_u(i, k) * t_dot_log_u(k, j);
1393 
1395  t_stress_diff(i, j, k, l) = 0;
1396  for (auto ii : {0, 1, 2})
1397  for (auto jj : {0, 1, 2})
1398  for (auto kk : {0, 1, 2})
1399  for (auto ll : {0, 1, 2})
1400  for (auto mm : {0, 1, 2})
1401  for (auto nn : {0, 1, 2})
1402  t_stress_diff(ii, jj, mm, nn) -=
1403  t_P_dh(ii, jj, kk, ll) * t_diff_u(kk, ll, mm, nn);
1404 
1406  t_dot_u_diff(i, j, k, l) = 0;
1407  for (auto ii : {0, 1, 2})
1408  for (auto jj : {0, 1, 2})
1409  for (auto kk : {0, 1, 2})
1410  for (auto ll : {0, 1, 2})
1411  for (auto mm : {0, 1, 2})
1412  t_dot_u_diff(ii, jj, kk, ll) +=
1413  t_diff_u(ii, mm, kk, ll) * t_dot_log_u(mm, jj) +
1414  ts_a * t_u(ii, mm) * t_one4_symmetric(mm, jj, kk, ll);
1415 
1416  auto t_viscous_stress = calculate_dot_deviator(t_dot_u);
1417  // t_viscous_stress(i, j) = t_dot_u(i, j);
1418  t_viscous_stress(i, j) *= alpha_u;
1419 
1420  auto t_viscous_stress_diff =
1421  calculate_diff_dot_deviator(t_dot_u, t_dot_u_diff);
1422  // t_viscous_stress_diff(i, j, k, l) = t_dot_u_diff(i, j, k, l);
1423  t_viscous_stress_diff(i, j, k, l) *= alpha_u;
1424 
1425  int rr = 0;
1426  for (; rr != row_nb_dofs / 6; ++rr) {
1427 
1428  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1429  auto t_m = get_ftensor4(K, 6 * rr, 0);
1430 
1431  for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1432 
1433  const double b = a * t_row_base_fun * t_col_base_fun;
1434 
1435  for (int ii : {0, 1, 2})
1436  for (int jj = 0; jj <= ii; ++jj)
1437  for (int kk : {0, 1, 2})
1438  for (int ll = 0; ll <= kk; ++ll)
1439  for (int mm : {0, 1, 2})
1440  for (int nn : {0, 1, 2})
1441  t_m(ii, jj, kk, ll) +=
1442  b * t_diff_u(mm, nn, ii, jj) *
1443  (t_stress_diff(mm, nn, kk, ll) - t_viscous_stress_diff(mm, nn, kk, ll));
1444 
1445  for (int ii : {0, 1, 2})
1446  for (int jj = 0; jj <= ii; ++jj) {
1447  int ss = 0;
1448  for (int kk : {0, 1, 2})
1449  for (int ll = 0; ll <= kk; ++ll, ++ss)
1450  for (int mm : {0, 1, 2})
1451  for (int nn : {0, 1, 2})
1452  t_m(ii, jj, kk, ll) +=
1453  b * (t_diff2_u[ss])(mm, nn, ii, jj) *
1454  (t_deltaP(mm, nn) - t_viscous_stress(mm, nn));
1455  }
1456 
1457  ++t_m;
1458  ++t_col_base_fun;
1459  }
1460 
1461  ++t_row_base_fun;
1462  }
1463  for (; rr != row_nb_base_functions; ++rr)
1464  ++t_row_base_fun;
1465 
1466  ++t_w;
1467  ++t_P_dh0;
1468  ++t_P_dh1;
1469  ++t_P_dh2;
1470  ++t_R;
1471  ++t_approx_P;
1472  ++t_P;
1473  ++t_dot_log_u;
1474  ++t_u;
1475  ++t_diff_u;
1476  for (auto ll : {0, 1, 2, 3, 4, 5})
1477  ++(t_diff2_u[ll]);
1478  }
1480 }
1481 
1482 MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
1483  EntData &col_data) {
1485  int row_nb_dofs = row_data.getIndices().size();
1486  int col_nb_dofs = col_data.getIndices().size();
1487  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1489 
1490  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1491 
1492  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1493 
1494  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1495 
1496  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1497 
1498  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1499 
1500  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1501 
1502  );
1503  };
1508 
1509  auto v = getVolume();
1510  auto t_w = getFTensor0IntegrationWeight();
1511  auto t_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
1512  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
1513  auto t_diff_u = FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>(
1514  &(*dataAtPts->diffStreachTensorAtPts)(0, 0),
1515  (*dataAtPts->diffStreachTensorAtPts).size2());
1516 
1517  int row_nb_base_functions = row_data.getN().size2();
1518  auto t_row_base_fun = row_data.getFTensor0N();
1519 
1520  int nb_integration_pts = row_data.getN().size1();
1521 
1522  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1523  double a = v * t_w;
1524 
1525  FTensor::Tensor3<double, 3, 3, 3> t_RT_domega_P;
1526  t_RT_domega_P(i, j, m) = t_diff_R(k, i, m) * t_P(k, j);
1528  t(i, j, k) = 0;
1529  for (int ii = 0; ii != 3; ++ii)
1530  for (int jj = 0; jj <= ii; ++jj)
1531  for (int mm = 0; mm != 3; ++mm)
1532  for (int nn = 0; nn != 3; ++nn)
1533  for (int kk = 0; kk != 3; ++kk)
1534  t(ii, jj, kk) +=
1535  t_diff_u(mm, nn, ii, jj) * t_RT_domega_P(mm, nn, kk);
1536 
1537  int rr = 0;
1538  for (; rr != row_nb_dofs / 6; ++rr) {
1539  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1540  auto t_m = get_ftensor3(K, 6 * rr, 0);
1541  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1542 
1543  double v = a * t_row_base_fun * t_col_base_fun;
1544  t_m(i, j, k) += v * t(i, j, k);
1545 
1546  ++t_m;
1547  ++t_col_base_fun;
1548  }
1549 
1550  ++t_row_base_fun;
1551  }
1552 
1553  for (; rr != row_nb_base_functions; ++rr)
1554  ++t_row_base_fun;
1555 
1556  ++t_w;
1557  ++t_P;
1558  ++t_diff_R;
1559  ++t_diff_u;
1560  }
1561 
1563 }
1564 
1565 MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
1566  EntData &col_data) {
1568  int nb_integration_pts = row_data.getN().size1();
1569  int row_nb_dofs = row_data.getIndices().size();
1570  int col_nb_dofs = col_data.getIndices().size();
1571  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1573 
1574  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1575 
1576  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1577 
1578  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1579 
1580  );
1581  };
1588 
1589  auto v = getVolume();
1590  auto t_w = getFTensor0IntegrationWeight();
1591  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1592 
1593  int row_nb_base_functions = row_data.getN().size2();
1594  auto t_row_base_fun = row_data.getFTensor0N();
1595 
1596  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1597  double a = v * t_w;
1598 
1599  int rr = 0;
1600  for (; rr != row_nb_dofs / 3; ++rr) {
1601  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1602  auto t_m = get_ftensor2(K, 3 * rr, 0);
1604  t_levi_h(i, j, k) =
1605  (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1606  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1607  FTensor::Tensor1<double, 3> t_pushed_base;
1608  t_m(k, i) += t_levi_h(i, j, k) * t_col_base_fun(j);
1609  ++t_m;
1610  ++t_col_base_fun;
1611  }
1612 
1613  ++t_row_base_fun;
1614  }
1615 
1616  for (; rr != row_nb_base_functions; ++rr)
1617  ++t_row_base_fun;
1618  ++t_w;
1619  ++t_R;
1620  }
1622 }
1623 
1624 MoFEMErrorCode OpSpatialRotation_domega_dBubble::integrate(EntData &row_data,
1625  EntData &col_data) {
1627  int nb_integration_pts = row_data.getN().size1();
1628  int row_nb_dofs = row_data.getIndices().size();
1629  int col_nb_dofs = col_data.getIndices().size();
1630  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1632  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
1633  };
1639 
1640  auto v = getVolume();
1641  auto t_w = getFTensor0IntegrationWeight();
1642  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1643 
1644  int row_nb_base_functions = row_data.getN().size2();
1645  auto t_row_base_fun = row_data.getFTensor0N();
1646 
1647  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1648  double a = v * t_w;
1649 
1650  int rr = 0;
1651  for (; rr != row_nb_dofs / 3; ++rr) {
1652 
1653  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1654  auto t_m = get_ftensor1(K, 3 * rr, 0);
1655 
1657  t_levi_h(i, j, k) =
1658  (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1659 
1660  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1661  FTensor::Tensor2<double, 3, 3> t_pushed_base;
1662  t_m(k) += t_levi_h(i, j, k) * t_col_base_fun(i, j);
1663  ++t_m;
1664  ++t_col_base_fun;
1665  }
1666 
1667  ++t_row_base_fun;
1668  }
1669 
1670  for (; rr != row_nb_base_functions; ++rr)
1671  ++t_row_base_fun;
1672  ++t_w;
1673  ++t_R;
1674  }
1676 }
1677 
1678 MoFEMErrorCode OpSpatialRotation_domega_domega::integrate(EntData &row_data,
1679  EntData &col_data) {
1681  int nb_integration_pts = row_data.getN().size1();
1682  int row_nb_dofs = row_data.getIndices().size();
1683  int col_nb_dofs = col_data.getIndices().size();
1684  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1686 
1687  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1688 
1689  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1690 
1691  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1692 
1693  );
1694  };
1701 
1702  auto v = getVolume();
1703  auto t_w = getFTensor0IntegrationWeight();
1704  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
1705  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1706  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
1707 
1708  int row_nb_base_functions = row_data.getN().size2();
1709  auto t_row_base_fun = row_data.getFTensor0N();
1710 
1711  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1712  double a = v * t_w;
1713 
1715  t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
1716 
1717  FTensor::Tensor2<double, 3, 3> t_leviPRT_domega;
1718  t_leviPRT_domega(n, l) =
1719  levi_civita(i, j, n) * (t_approx_P(i, k) * t_diff_R(j, k, l));
1720 
1721  int rr = 0;
1722  for (; rr != row_nb_dofs / 3; ++rr) {
1723 
1724  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1725  auto t_m = get_ftensor2(K, 3 * rr, 0);
1726  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1727 
1728  const double b = a * t_row_base_fun * t_col_base_fun;
1729 
1730  t_m(k, m) += b * t_leviPRT_domega(k, m);
1731 
1732  ++t_m;
1733  ++t_col_base_fun;
1734  }
1735 
1736  ++t_row_base_fun;
1737  }
1738 
1739  for (; rr != row_nb_base_functions; ++rr)
1740  ++t_row_base_fun;
1741  ++t_w;
1742  ++t_approx_P;
1743  ++t_R;
1744  ++t_diff_R;
1745  }
1747 }
1748 
1749 MoFEMErrorCode OpSpatialRotation_domega_du::integrate(EntData &row_data,
1750  EntData &col_data) {
1752  // int nb_integration_pts = row_data.getN().size1();
1753  // int row_nb_dofs = row_data.getIndices().size();
1754  // int col_nb_dofs = col_data.getIndices().size();
1755 
1756  // auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1757  // return FTensor::Christof<FTensor::PackPtr<double *, 6>, 3, 3>(
1758 
1759  // &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c +
1760  // 3), &m(r + 0, c + 4), &m(r + 0, c + 5),
1761 
1762  // &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c +
1763  // 3), &m(r + 1, c + 4), &m(r + 1, c + 5),
1764 
1765  // &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c +
1766  // 3), &m(r + 2, c + 4), &m(r + 2, c + 5)
1767 
1768  // );
1769  // };
1770 
1771  // FTensor::Index<'i', 3> i;
1772  // FTensor::Index<'j', 3> j;
1773  // FTensor::Index<'k', 3> k;
1774  // FTensor::Index<'l', 3> l;
1775  // FTensor::Index<'m', 3> m;
1776  // FTensor::Index<'n', 3> n;
1777 
1778  // auto v = getVolume();
1779  // auto t_w = getFTensor0IntegrationWeight();
1780  // auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
1781  // auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1782  // auto t_u =
1783  // getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts)); auto
1784 
1785  // int row_nb_base_functions = row_data.getN().size2();
1786  // auto t_row_base_fun = row_data.getFTensor0N();
1787 
1788  // for (int gg = 0; gg != nb_integration_pts; ++gg) {
1789  // double a = v * t_w;
1790 
1791  // FTensor::Tensor4<double, 3, 3, 3, 3> t_PhT_du;
1792  // t_PhT_du(i, j, k, l) = t_approx_P(i, l) * t_R(j, k);
1793 
1794  // FTensor::Tensor3<double, 3, 3, 3> t_levi;
1795  // t_levi(i, j, k) = levi_civita(i, j, k);
1796  // FTensor::Tensor3<double, 3, 3, 3> t_leviPhT_du;
1797  // t_leviPhT_du(n, k, l) = 0;
1798  // for (int ii : {0, 1, 2})
1799  // for (int jj : {0, 1, 2})
1800  // for (int nn : {0, 1, 2})
1801  // for (int kk : {0, 1, 2})
1802  // for (int ll : {0, 1, 2})
1803  // t_leviPhT_du(nn, kk, ll) += t_levi(ii, jj, nn) *
1804  // t_PhT_du(ii, jj, kk, ll);
1805 
1806  // int rr = 0;
1807  // for (; rr != row_nb_dofs / 3; ++rr) {
1808 
1809  // auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1810  // auto t_m = get_ftensor3(K, 3 * rr, 0);
1811  // for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1812 
1813  // const double b = a * t_row_base_fun * t_col_base_fun;
1814 
1815  // for (int ii : {0, 1, 2})
1816  // for (int jj : {0, 1, 2})
1817  // for (int kk : {0, 1, 2})
1818  // t_m(ii, jj, kk) += b * t_leviPhT_du(ii, jj, kk);
1819 
1820  // ++t_m;
1821  // ++t_col_base_fun;
1822  // }
1823 
1824  // ++t_row_base_fun;
1825  // }
1826 
1827  // for (; rr != row_nb_base_functions; ++rr)
1828  // ++t_row_base_fun;
1829  // ++t_w;
1830  // ++t_approx_P;
1831  // ++t_R;
1832  // ++t_u;
1833  // }
1835 }
1836 
1837 MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
1838  EntData &col_data) {
1840 
1841  int nb_integration_pts = row_data.getN().size1();
1842  int row_nb_dofs = row_data.getIndices().size();
1843  int col_nb_dofs = col_data.getIndices().size();
1844 
1845  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1847 
1848  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1849 
1850  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1851 
1852  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1853 
1854  );
1855  };
1856 
1861 
1862  auto v = getVolume();
1863  auto t_w = getFTensor0IntegrationWeight();
1864  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
1865  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
1866  int row_nb_base_functions = row_data.getN().size2() / 3;
1867  auto t_row_base_fun = row_data.getFTensor1N<3>();
1868 
1869  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1870  double a = v * t_w;
1871 
1872  int rr = 0;
1873  for (; rr != row_nb_dofs / 3; ++rr) {
1874 
1876  t_PRT(i, k) = t_row_base_fun(j) * (t_diff_R(i, l, k) * t_u(l, j));
1877 
1878  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1879  auto t_m = get_ftensor2(K, 3 * rr, 0);
1880  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1881  t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
1882  ++t_m;
1883  ++t_col_base_fun;
1884  }
1885 
1886  ++t_row_base_fun;
1887  }
1888 
1889  for (; rr != row_nb_base_functions; ++rr)
1890  ++t_row_base_fun;
1891  ++t_w;
1892  ++t_diff_R;
1893  ++t_u;
1894  }
1896 }
1897 
1899 OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
1900  EntData &col_data) {
1902 
1903  int nb_integration_pts = row_data.getN().size1();
1904  int row_nb_dofs = row_data.getIndices().size();
1905  int col_nb_dofs = col_data.getIndices().size();
1906 
1907  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1909  &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
1910  };
1911 
1916 
1917  auto v = getVolume();
1918  auto t_w = getFTensor0IntegrationWeight();
1919  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
1920  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
1921  int row_nb_base_functions = row_data.getN().size2() / 9;
1922  auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
1923 
1924  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1925  double a = v * t_w;
1926 
1927  int rr = 0;
1928  for (; rr != row_nb_dofs; ++rr) {
1929 
1931  t_PRT(k) = t_row_base_fun(i, j) * (t_diff_R(i, l, k) * t_u(l, j));
1932 
1933  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1934  auto t_m = get_ftensor2(K, rr, 0);
1935  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1936  t_m(j) += (a * t_col_base_fun) * t_PRT(j);
1937  ++t_m;
1938  ++t_col_base_fun;
1939  }
1940 
1941  ++t_row_base_fun;
1942  }
1943 
1944  for (; rr != row_nb_base_functions; ++rr)
1945  ++t_row_base_fun;
1946 
1947  ++t_w;
1948  ++t_diff_R;
1949  ++t_u;
1950  }
1952 }
1953 
1954 MoFEMErrorCode OpSpatialSchurBegin::assemble(int row_side, EntityType row_type,
1955  EntData &data) {
1957  if (row_type != MBTET)
1959  auto &bmc = dataAtPts->blockMatContainor;
1960  for (auto &bit : bmc)
1961  bit.unSetAtElement();
1962  // bmc.clear();
1964 }
1965 
1966 MoFEMErrorCode OpSpatialPreconditionMass::integrate(EntData &row_data) {
1968  const int nb_integration_pts = row_data.getN().size1();
1969  const int row_nb_dofs = row_data.getIndices().size();
1970  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1972  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
1973 
1974  );
1975  };
1977 
1978  auto v = getVolume();
1979  auto t_w = getFTensor0IntegrationWeight();
1980  auto &m = *mPtr;
1981 
1982  m.resize(row_nb_dofs, row_nb_dofs, false);
1983  m.clear();
1984 
1985  int row_nb_base_functions = row_data.getN().size2();
1986  auto t_row_base_fun = row_data.getFTensor0N();
1987 
1988  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1989  double a = v * t_w;
1990 
1991  int rr = 0;
1992  for (; rr != row_nb_dofs / 3; ++rr) {
1993 
1994  auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
1995  auto t_m = get_ftensor1(m, 3 * rr, 0);
1996  for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1997  const double b = a * t_row_base_fun * t_col_base_fun;
1998  t_m(i) += b;
1999  ++t_m;
2000  ++t_col_base_fun;
2001  }
2002 
2003  ++t_row_base_fun;
2004  }
2005 
2006  for (; rr != row_nb_base_functions; ++rr)
2007  ++t_row_base_fun;
2008 
2009  ++t_w;
2010  }
2011 
2013 }
2014 
2015 MoFEMErrorCode OpSpatialSchurEnd::assemble(int row_side, EntityType row_type,
2016  EntData &data) {
2018  if (row_type != MBTET)
2020 
2021  if (auto ep_fe_vol_ptr =
2022  dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
2023  getFEMethod())) {
2024 
2025  AO aoSuu = PETSC_NULL;
2026  AO aoSBB = PETSC_NULL;
2027  AO aoSOO = PETSC_NULL;
2028  AO aoSww = PETSC_NULL;
2029 
2030  Mat Suu = PETSC_NULL;
2031  Mat SBB = PETSC_NULL;
2032  Mat SOO = PETSC_NULL;
2033  Mat Sww = PETSC_NULL;
2034 
2035  auto set_data = [&](auto fe) {
2036  aoSuu = fe->aoSuu;
2037  aoSBB = fe->aoSBubble;
2038  aoSOO = fe->aoSOmega;
2039  aoSww = fe->aoSw;
2040 
2041  Suu = fe->Suu;
2042  SBB = fe->SBubble;
2043  SOO = fe->SOmega;
2044  Sww = fe->Sw;
2045  };
2046  set_data(ep_fe_vol_ptr);
2047 
2048  if (Suu) {
2049 
2050  auto find_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2051  auto &row_name, auto &col_name, auto row_side,
2052  auto col_side, auto row_type, auto col_type) {
2053  return add_bmc.get<0>().find(boost::make_tuple(
2054  row_name, col_name, row_type, col_type, row_side, col_side));
2055  };
2056 
2057  auto set_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2058  auto &row_name, auto &col_name, auto row_side,
2059  auto col_side, auto row_type, auto col_type,
2060  const auto &m, const auto &row_ind,
2061  const auto &col_ind) {
2062  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2063  row_type, col_type);
2064  if (it != add_bmc.get<0>().end()) {
2065  it->setMat(m);
2066  it->setInd(row_ind, col_ind);
2067  it->setSetAtElement();
2068  return it;
2069  } else {
2070  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2071  row_name, col_name, row_type, col_type, row_side, col_side, m,
2072  row_ind, col_ind));
2073  return p_it.first;
2074  }
2075  };
2076 
2077  auto add_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2078  auto &row_name, auto &col_name, auto row_side,
2079  auto col_side, auto row_type, auto col_type,
2080  const auto &m, const auto &row_ind,
2081  const auto &col_ind) {
2082  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2083  row_type, col_type);
2084  if (it != add_bmc.get<0>().end()) {
2085  it->addMat(m);
2086  it->setInd(row_ind, col_ind);
2087  it->setSetAtElement();
2088  return it;
2089  } else {
2090  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2091  row_name, col_name, row_type, col_type, row_side, col_side, m,
2092  row_ind, col_ind));
2093  return p_it.first;
2094  }
2095  };
2096 
2097  auto assemble_block = [&](auto &bit, Mat S) {
2099  const VectorInt &rind = bit.rowInd;
2100  const VectorInt &cind = bit.colInd;
2101  const MatrixDouble &m = bit.M;
2102  CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
2103  &*cind.begin(), &*m.data().begin(), ADD_VALUES);
2104 
2106  };
2107 
2108  auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
2110  const int nb = m.size1();
2111 
2112  inv.resize(nb, nb, false);
2113  inv.clear();
2114  for (int cc = 0; cc != nb; ++cc)
2115  inv(cc, cc) = -1;
2116 
2117  iPIV.resize(nb, false);
2118  lapackWork.resize(nb * nb, false);
2119  const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
2120  &*iPIV.begin(), &*inv.data().begin(), nb,
2121  &*lapackWork.begin(), nb * nb);
2122  // if (info != 0)
2123  // SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2124  // "Can not invert matrix info = %d", info);
2125 
2127  };
2128 
2129  auto invert_symm_schur = [&](DataAtIntegrationPts::BlockMatContainor &bmc,
2130  std::string field, auto &inv) {
2132 
2133  auto bit =
2134  bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
2135  if (bit != bmc.get<1>().end()) {
2136 
2137  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2138  CHKERR invert_symm_mat(m, inv);
2139 
2140  } else
2141  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2142  "%s matrix not found", field.c_str());
2143 
2145  };
2146 
2147  auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
2149 
2150  const int nb = m.size1();
2151 
2152  MatrixDouble trans_m = trans(m);
2153  MatrixDouble trans_inv;
2154  trans_inv.resize(nb, nb, false);
2155  trans_inv.clear();
2156  for (int c = 0; c != nb; ++c)
2157  trans_inv(c, c) = -1;
2158 
2159  iPIV.resize(nb, false);
2160  const auto info =
2161  lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb, &*iPIV.begin(),
2162  &*trans_inv.data().begin(), nb);
2163  if (info != 0)
2164  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2165  "Can not invert matrix info = %d", info);
2166 
2167  inv.resize(nb, nb, false);
2168  noalias(inv) = trans(trans_inv);
2169 
2171  };
2172 
2173  auto invert_nonsymm_schur =
2174  [&](DataAtIntegrationPts::BlockMatContainor &bmc, std::string field,
2175  auto &inv, const bool debug = false) {
2177 
2178  auto bit = bmc.get<1>().find(
2179  boost::make_tuple(field, field, MBTET, MBTET));
2180  if (bit != bmc.get<1>().end()) {
2181 
2182  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2183  CHKERR invert_nonsymm_mat(m, inv);
2184 
2185  if (debug) {
2186  std::cerr << prod(m, inv) << endl;
2187  std::cerr << endl;
2188  }
2189 
2190  } else
2191  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2192  "%s matrix not found", field.c_str());
2193 
2195  };
2196 
2197  auto create_block_schur =
2200  std::string field, AO ao, MatrixDouble &inv) {
2202 
2203  for (auto &bit : add_bmc) {
2204  bit.unSetAtElement();
2205  bit.clearMat();
2206  }
2207 
2208  for (auto &bit : bmc) {
2209  if (bit.setAtElement && bit.rowField != field &&
2210  bit.colField != field) {
2211  VectorInt rind = bit.rowInd;
2212  VectorInt cind = bit.colInd;
2213  const MatrixDouble &m = bit.M;
2214  if (ao) {
2215  CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
2216  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2217  }
2218  auto it = set_block(add_bmc, bit.rowField, bit.colField,
2219  bit.rowSide, bit.colSide, bit.rowType,
2220  bit.colType, m, rind, cind);
2221  }
2222  }
2223 
2224  for (auto &bit_col : bmc) {
2225  if (bit_col.setAtElement && bit_col.rowField == field &&
2226  bit_col.colField != field) {
2227  const MatrixDouble &cm = bit_col.M;
2228  VectorInt cind = bit_col.colInd;
2229  invMat.resize(inv.size1(), cm.size2(), false);
2230  noalias(invMat) = prod(inv, cm);
2231  if (ao)
2232  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2233  for (auto &bit_row : bmc) {
2234  if (bit_row.setAtElement && bit_row.rowField != field &&
2235  bit_row.colField == field) {
2236  const MatrixDouble &rm = bit_row.M;
2237  VectorInt rind = bit_row.rowInd;
2238  K.resize(rind.size(), cind.size(), false);
2239  noalias(K) = prod(rm, invMat);
2240  if (ao)
2241  CHKERR AOApplicationToPetsc(ao, rind.size(),
2242  &*rind.begin());
2243  auto it = add_block(add_bmc, bit_row.rowField,
2244  bit_col.colField, bit_row.rowSide,
2245  bit_col.colSide, bit_row.rowType,
2246  bit_col.colType, K, rind, cind);
2247  }
2248  }
2249  }
2250  }
2251 
2253  };
2254 
2255  auto assemble_schur =
2256  [&](DataAtIntegrationPts::BlockMatContainor &add_bmc, Mat S,
2257  bool debug = false) {
2259  for (auto &bit : add_bmc) {
2260  if (bit.setAtElement)
2261  CHKERR assemble_block(bit, S);
2262  }
2263  if (debug) {
2264  for (auto &bit : add_bmc) {
2265  if (bit.setAtElement) {
2266  std::cerr << "assemble: " << bit.rowField << " "
2267  << bit.colField << endl;
2268  std::cerr << bit.M << endl;
2269  }
2270  }
2271  std::cerr << std::endl;
2272  }
2274  };
2275 
2276  auto precondition_schur =
2279  const std::string field, const MatrixDouble &diag_mat,
2280  const double eps) {
2282 
2283  for (auto &bit : add_bmc) {
2284  bit.unSetAtElement();
2285  bit.clearMat();
2286  }
2287 
2288  for (auto &bit : bmc) {
2289  if (bit.setAtElement) {
2290  if (bit.rowField != field || bit.colField != field)
2291  auto it =
2292  set_block(add_bmc, bit.rowField, bit.colField,
2293  bit.rowSide, bit.colSide, bit.rowType,
2294  bit.colType, bit.M, bit.rowInd, bit.colInd);
2295  }
2296  }
2297 
2298  auto bit = bmc.get<1>().find(
2299  boost::make_tuple(field, field, MBTET, MBTET));
2300  if (bit->setAtElement && bit != bmc.get<1>().end()) {
2301  auto it =
2302  set_block(add_bmc, bit->rowField, bit->colField, bit->rowSide,
2303  bit->colSide, bit->rowType, bit->colType, bit->M,
2304  bit->rowInd, bit->colInd);
2305  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2306  m += eps * diag_mat;
2307  } else {
2308  auto row_it = bmc.get<3>().lower_bound(field);
2309  for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
2310  if (row_it->setAtElement) {
2311  auto it = set_block(add_bmc, field, field, 0, 0, MBTET, MBTET,
2312  diag_mat, row_it->rowInd, row_it->rowInd);
2313  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2314  m *= eps;
2315  break;
2316  }
2317  }
2318  if (row_it == bmc.get<3>().end())
2319  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2320  "row field not found %s", field.c_str());
2321  }
2322 
2324  };
2325 
2326  CHKERR invert_symm_schur(dataAtPts->blockMatContainor, "u",
2327  invBlockMat["uu"]);
2328  CHKERR create_block_schur(dataAtPts->blockMatContainor, blockMat["uu"],
2329  "u", aoSuu, invBlockMat["uu"]);
2330  CHKERR assemble_schur(blockMat["uu"], Suu);
2331 
2332  if (SBB) {
2333  CHKERR invert_symm_schur(blockMat["uu"], "BUBBLE", invBlockMat["BB"]);
2334  CHKERR create_block_schur(blockMat["uu"], blockMat["BB"], "BUBBLE",
2335  aoSBB, invBlockMat["BB"]);
2336  CHKERR precondition_schur(blockMat["BB"], blockMat["precBBOO"], "omega",
2337  *dataAtPts->ooMatPtr, eps);
2338  CHKERR assemble_schur(blockMat["precBBOO"], SBB);
2339 
2340  if (SOO) {
2341  CHKERR invert_nonsymm_schur(blockMat["precBBOO"], "omega",
2342  invBlockMat["OO"]);
2343  CHKERR create_block_schur(blockMat["precBBOO"], blockMat["OO"],
2344  "omega", aoSOO, invBlockMat["OO"]);
2345  if (dataAtPts->wwMatPtr) {
2346  CHKERR precondition_schur(blockMat["OO"], blockMat["precOOww"], "w",
2347  *dataAtPts->wwMatPtr, -eps);
2348  } else {
2349  blockMat["precOOww"] = blockMat["OO"];
2350  }
2351  CHKERR assemble_schur(blockMat["precOOww"], SOO);
2352 
2353  if (Sww) {
2354  CHKERR invert_symm_schur(blockMat["precOOww"], "w",
2355  invBlockMat["ww"]);
2356  CHKERR create_block_schur(blockMat["precOOww"], blockMat["ww"], "w",
2357  aoSww, invBlockMat["ww"]);
2358  CHKERR assemble_schur(blockMat["ww"], Sww);
2359  }
2360  }
2361  }
2362  }
2363  }
2364 
2366 }
2367 
2368 MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
2369  EntData &data) {
2371  if (type != MBVERTEX)
2373 
2374  auto create_tag = [this](const std::string tag_name, const int size) {
2375  double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2376  Tag th;
2377  CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
2378  th, MB_TAG_CREAT | MB_TAG_SPARSE,
2379  def_VAL);
2380  return th;
2381  };
2382 
2383  Tag th_w = create_tag("SpatialDisplacement", 3);
2384  Tag th_omega = create_tag("Omega", 3);
2385  Tag th_approxP = create_tag("Piola1Stress", 9);
2386  Tag th_sigma = create_tag("CauchyStress", 9);
2387  Tag th_log_u = create_tag("LogSpatialStreach", 9);
2388  Tag th_u = create_tag("SpatialStreach", 9);
2389  Tag th_h = create_tag("h", 9);
2390  Tag th_X = create_tag("X", 3);
2391  Tag th_detF = create_tag("detF", 1);
2392  Tag th_angular_momentum = create_tag("AngularMomentum", 3);
2393 
2394  int nb_gauss_pts = data.getN().size1();
2395  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
2396  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2397  "Nb. of integration points is not equal to number points on "
2398  "post-processing mesh");
2399  }
2400 
2401  auto t_w = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
2402  auto t_omega = getFTensor1FromMat<3>(*(dataAtPts->rotAxisAtPts));
2403  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
2404  auto t_log_u =
2405  getFTensor2SymmetricFromMat<3>(*(dataAtPts->logStreachTensorAtPts));
2406  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
2407  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
2408  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
2409  auto t_coords = getFTensor1CoordsAtGaussPts();
2410 
2415 
2416  // vectors
2417  VectorDouble3 v(3);
2418  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
2419  auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
2421  t_v(i) = t_d(i);
2422  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2423  &*v.data().begin());
2425  };
2426 
2427  MatrixDouble3by3 m(3, 3);
2429  &m(0, 0), &m(0, 1), &m(0, 2),
2430 
2431  &m(1, 0), &m(1, 1), &m(1, 2),
2432 
2433  &m(2, 0), &m(2, 1), &m(2, 2));
2434 
2435  auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
2437  t_m(i, j) = t_d(i, j);
2438  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2439  &*m.data().begin());
2441  };
2442 
2443  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2444 
2445  // vetors
2446  CHKERR save_vec_tag(th_w, t_w, gg);
2447  CHKERR save_vec_tag(th_X, t_coords, gg);
2448  CHKERR save_vec_tag(th_omega, t_omega, gg);
2449 
2450  // tensors
2451  CHKERR save_mat_tag(th_h, t_h, gg);
2452 
2453  FTensor::Tensor2<double, 3, 3> t_log_u_tmp;
2454  for (int ii = 0; ii != 3; ++ii)
2455  for (int jj = 0; jj != 3; ++jj)
2456  t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
2457 
2458  CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
2459 
2461  for (int ii = 0; ii != 3; ++ii)
2462  for (int jj = 0; jj != 3; ++jj)
2463  t_u_tmp(ii, jj) = t_u(ii, jj);
2464 
2465  CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2466  CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2467 
2468  const double jac = dEterminant(t_h);
2470  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
2471  CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2472  CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
2473 
2475  t_PhT(i, k) = t_approx_P(i, j) * t_R(k, j);
2476  FTensor::Tensor1<double, 3> t_leviPRT;
2477  t_leviPRT(k) = levi_civita(i, l, k) * t_PhT(i, l);
2478 
2479  CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
2480  &t_leviPRT(0));
2481 
2482  ++t_w;
2483  ++t_h;
2484  ++t_log_u;
2485  ++t_u;
2486  ++t_omega;
2487  ++t_R;
2488  ++t_approx_P;
2489  ++t_coords;
2490  }
2491 
2493 }
2494 
2496 OpCalculateStrainEnergy::doWork(int side, EntityType type,
2499  int nb_integration_pts = data.getN().size1();
2500  auto v = getVolume();
2501  auto t_w = getFTensor0IntegrationWeight();
2502  auto t_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
2503  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
2504 
2507 
2508  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2509  const double a = t_w * v;
2510  (*energy) += a * t_P(i, J) * t_h(i, J);
2511  ++t_w;
2512  ++t_P;
2513  ++t_h;
2514  }
2515 
2517 }
2518 
2519 } // namespace EshelbianPlasticity
static MoFEMErrorCode matrix_exponent_derivative(MatrixDouble &log_u, MatrixDouble &diff_u, std::array< MatrixDouble, 6 > &diff2_u)
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:506
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
Fully Antisymmetric Levi-Civita Tensor.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:71
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)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
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 double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:452
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
static const bool debug
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > calculate_dot_deviator(FTensor::Tensor2< T, 3, 3 > &t)
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:97
#define CHKERR
Inline error check.
Definition: definitions.h:601
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:71
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:87
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)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
static FTensor::Tensor4< typename TensorTypeExtractor< T >::Type, 3, 3, 3, 3 > calculate_diff_dot_deviator(FTensor::Tensor2< T, 3, 3 > &t, FTensor::Tensor4< T, 3, 3, 3, 3 > &t_diff)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
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
const int N
Definition: speed_test.cpp:3
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.
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
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.
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
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