v0.14.0
EshelbianOperators.cpp
Go to the documentation of this file.
1 /**
2  * \file EshelbianOperators.cpp
3  * \example EshelbianOperators.cpp
4  *
5  * \brief Implementation of operators
6  */
7 
8 #include <MoFEM.hpp>
9 using namespace MoFEM;
10 
11 #include <EshelbianPlasticity.hpp>
12 
13 #include <boost/math/constants/constants.hpp>
14 
15 #include <EshelbianAux.hpp>
16 
17 #include <lapack_wrap.h>
18 
19 #include <MatrixFunction.hpp>
20 
21 namespace EshelbianPlasticity {
22 
23 template <typename T>
25  RotSelector rotSelector = LARGE_ROT) {
26 
27  using D = typename TensorTypeExtractor<T>::Type;
28 
33 
34  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
35  t_R(i, j) = t_kd(i, j);
36 
38  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
39 
40  if (rotSelector == SMALL_ROT) {
41  t_R(i, j) += t_Omega(i, j);
42  return t_R;
43  }
44 
45  const auto angle =
46  sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
47  const auto a = sin(angle) / angle;
48  t_R(i, j) += a * t_Omega(i, j);
49  if (rotSelector == MODERATE_ROT)
50  return t_R;
51 
52  const auto ss_2 = sin(angle / 2.) / angle;
53  const auto b = 2. * ss_2 * ss_2;
54  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
55 
56  return t_R;
57 }
58 
59 template <typename T>
61  RotSelector rotSelector = LARGE_ROT) {
62 
63  using D = typename TensorTypeExtractor<T>::Type;
64 
69 
71 
72  const auto angle =
73  sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
74  if (rotSelector == SMALL_ROT) {
75  t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, j, k);
76  return t_diff_R;
77  }
78 
79  const auto ss = sin(angle);
80  const auto a = ss / angle;
81  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
82 
84  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
85 
86  const auto angle2 = angle * angle;
87  const auto cc = cos(angle);
88  const auto diff_a = (angle * cc - ss) / angle2;
89  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
90  if (rotSelector == MODERATE_ROT)
91  return t_diff_R;
92 
93  const auto ss_2 = sin(angle / 2.);
94  const auto cc_2 = cos(angle / 2.);
95  const auto b = 2. * ss_2 * ss_2 / angle2;
96  t_diff_R(i, j, k) +=
97  b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
98  FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
99  const auto diff_b =
100  (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
101  t_diff_R(i, j, k) +=
102  diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
103 
104  return t_diff_R;
105 }
106 
107 template <typename T>
109  RotSelector rotSelector = LARGE_ROT) {
110 
111  using D = typename TensorTypeExtractor<T>::Type;
112 
115 
116  constexpr double eps = 1e-10;
117  for (int l = 0; l != 3; ++l) {
119  t_omega_c(i) = t_omega(i);
120  t_omega_c(l) += std::complex<double>(0, eps);
121  auto t_diff_R_c = get_diff_rotation_form_vector(t_omega_c, rotSelector);
122  for (int i = 0; i != 3; ++i) {
123  for (int j = 0; j != 3; ++j) {
124  for (int k = 0; k != 3; ++k) {
125  t_diff2_R(i, j, k, l) = t_diff_R_c(i, j, k).imag() / eps;
126  }
127  }
128  }
129  }
130 
131  return t_diff2_R;
132 }
133 
134 struct isEq {
135  static inline auto check(const double &a, const double &b) {
136  constexpr double eps = std::numeric_limits<float>::epsilon();
137  return std::abs(a - b) / absMax < eps;
138  }
139  static double absMax;
140 };
141 
142 double isEq::absMax = 1;
143 
144 inline auto is_eq(const double &a, const double &b) {
145  return isEq::check(a, b);
146 };
147 
148 template <int DIM> inline auto get_uniq_nb(double *ptr) {
149  std::array<double, DIM> tmp;
150  std::copy(ptr, &ptr[DIM], tmp.begin());
151  std::sort(tmp.begin(), tmp.end());
152  isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
153  constexpr double eps = std::numeric_limits<float>::epsilon();
154  isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
155  return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
156 }
157 
158 template <int DIM>
161 
162  isEq::absMax =
163  std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
164 
165  int i = 0, j = 1, k = 2;
166 
167  if (is_eq(eig(0), eig(1))) {
168  i = 0;
169  j = 2;
170  k = 1;
171  } else if (is_eq(eig(0), eig(2))) {
172  i = 0;
173  j = 1;
174  k = 2;
175  } else if (is_eq(eig(1), eig(2))) {
176  i = 1;
177  j = 0;
178  k = 2;
179  }
180 
181  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
182  eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
183 
184  eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
185 
186  eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
187 
188  FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
189 
190  {
193  eig(i) = eig_c(i);
194  eigen_vec(i, j) = eigen_vec_c(i, j);
195  }
196 }
197 
198 
199 
201  EntData &data) {
206 
207  int nb_integration_pts = getGaussPts().size2();
208 
209  auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
210  auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
211  auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
212 
213  dataAtPts->SigmaAtPts.resize(SPACE_DIM * SPACE_DIM, nb_integration_pts,
214  false);
215  auto t_eshelby_stress =
216  getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->SigmaAtPts);
217 
219 
220  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
221  t_eshelby_stress(i, j) = t_energy * t_kd(i, j) - t_F(m, i) * t_P(m, j);
222  ++t_energy;
223  ++t_P;
224  ++t_F;
225  ++t_eshelby_stress;
226  }
227 
229 }
230 
231 MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
232  EntityType type,
233  EntData &data) {
235 
237  auto t_L = symm_L_tensor();
238 
239  int nb_integration_pts = getGaussPts().size2();
246 
247  dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
248  dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts, false);
249  dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts, false);
250  dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts, false);
251  dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts, false);
252  dataAtPts->leviKirchhoffOmegaAtPts.resize(9, nb_integration_pts, false);
253 
254  dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts, false);
255  dataAtPts->adjointPdUAtPts.resize(size_symm, nb_integration_pts, false);
256  dataAtPts->adjointPdUdPAtPts.resize(9 * size_symm, nb_integration_pts, false);
257  dataAtPts->adjointPdUdOmegaAtPts.resize(3 * size_symm, nb_integration_pts,
258  false);
259  dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
260  dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
261  dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts, false);
262  dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts, false);
263  dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
264  dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
265  dataAtPts->nbUniq.resize(nb_integration_pts, false);
266 
267  dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts, false);
268  dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
269 
270  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
271  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
272  auto t_h_dlog_u =
273  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
274  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
275  auto t_levi_kirchoff_dP =
276  getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
277  auto t_levi_kirchoff_domega =
278  getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
279  auto t_approx_P_adjont_dstretch =
280  getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
281  auto t_approx_P_adjont_log_du =
282  getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
283  auto t_approx_P_adjont_log_du_dP =
284  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
285  auto t_approx_P_adjont_log_du_domega =
286  getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
287  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
288  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
289  auto t_diff_R = getFTensor3FromMat<3, 3, 3>(dataAtPts->diffRotMatAtPts);
290  auto t_log_u =
291  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
292  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
293  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
294  auto t_diff_u =
295  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
296  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
297  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
298 
299  auto &nbUniq = dataAtPts->nbUniq;
300  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
301  auto t_log_stretch_total =
302  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
303  auto t_log_u2_h1 =
304  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
306 
307  for (int gg = 0; gg != nb_integration_pts; ++gg) {
308 
314  FTensor::Tensor2<double, 3, 3> t_approx_P_intermidiate;
315 
316  switch (EshelbianCore::gradApproximator) {
317  case NO_H1_CONFIGURATION:
318  t_h1(i, j) = t_kd(i, j);
319  break;
320  case LARGE_ROT:
321  case MODERATE_ROT:
322  t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
323  break;
324  case SMALL_ROT:
325  t_h1(i, j) = t_kd(i, j);
326  break;
327  };
328 
329  auto calculate_rotation = [&]() {
330  auto t0_diff =
331  get_diff_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
332  auto t0 = get_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
333  t_diff_R(i, j, k) = t0_diff(i, j, k);
334  t_R(i, j) = t0(i, j);
335  };
336 
337  auto calculate_stretch = [&]() {
338  if (EshelbianCore::stretchSelector == StretchSelector::LINEAR) {
339  constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
340  t_u(i, j) = t_log_u(i, j) + t_kd_sym(i, j);
341  t_diff_u(i, j, k, l) = (t_kd_sym(i, k) ^ t_kd_sym(j, l)) / 4.;
342  t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
343  } else {
344  eigen_vec(i, j) = t_log_u(i, j);
345  CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
346  // rare case when two eigen values are equal
347  nbUniq[gg] = get_uniq_nb<3>(&eig(0));
348  if (nbUniq[gg] < 3) {
349  sort_eigen_vals<3>(eig, eigen_vec);
350  }
351  t_eigen_vals(i) = eig(i);
352  t_eigen_vecs(i, j) = eigen_vec(i, j);
353  auto t_u_tmp =
354  EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f);
355  t_u(i, j) = t_u_tmp(i, j);
356  auto t_diff_u_tmp = EigenMatrix::getDiffMat(
357  t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
358  nbUniq[gg]);
359  t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
360  t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
361  }
362  };
363 
364  calculate_rotation();
365 
366  if (!EshelbianCore::noStretch) {
367 
368  calculate_stretch();
369 
370  switch (EshelbianCore::gradApproximator) {
371  case NO_H1_CONFIGURATION:
372  case LARGE_ROT:
373 
374  // rotation can be large, moderate or small
375  t_Ru(i, m) = t_R(i, l) * t_u(l, m);
376  t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
377  t_h_domega(i, j, k) = (t_diff_R(i, l, k) * t_u(l, m)) * t_h1(m, j);
378  t_h_dlog_u(i, j, L) = (t_R(i, l) * t_Ldiff_u(l, m, L)) * t_h1(m, j);
379 
380  // conservation of angular momentum at current configuration
381  t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
382  t_approx_P_adjont_dstretch(l, m) =
383  t_approx_P_intermidiate(i, m) * t_R(i, l);
384 
385  t_levi_kirchoff(k) =
386  levi_civita(l, m, k) * (t_approx_P_adjont_dstretch(l, m));
387  t_levi_kirchoff_dP(i, j, k) =
388  (levi_civita(l, m, k) * t_R(i, l)) * t_h1(m, j);
389  t_levi_kirchoff_domega(k, n) =
390  levi_civita(l, m, k) *
391  (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
392 
393  t_approx_P_adjont_log_du(L) =
394  t_Ldiff_u(l, m, L) * t_approx_P_adjont_dstretch(l, m);
395  t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
396  t_approx_P_adjont_log_du_domega(n, L) =
397  t_Ldiff_u(l, m, L) *
398  (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
399 
400  break;
401  case MODERATE_ROT:
402 
403  {
404 
405  // assume small current rotation
407  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
408  t_Ru(i, m) = t_Omega(i, m) + t_u(i, m);
409  t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
410  t_h_domega(i, j, k) =
411  FTensor::levi_civita<double>(i, m, k) * t_h1(m, j);
412  t_h_dlog_u(i, j, L) = t_Ldiff_u(i, m, L) * t_h1(m, j);
413 
414  // conservation of angular momentum at intermediate configuration
415  t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
416  t_approx_P_adjont_dstretch(i, m) = t_approx_P_intermidiate(i, m);
417 
418  t_levi_kirchoff(k) =
419  levi_civita(i, m, k) * (t_approx_P_adjont_dstretch(i, m));
420  t_levi_kirchoff_dP(i, j, k) = levi_civita(i, m, k) * t_h1(m, j);
421  t_levi_kirchoff_domega(k, n) = 0;
422 
423  t_approx_P_adjont_log_du(L) =
424  t_Ldiff_u(i, m, L) * t_approx_P_adjont_dstretch(i, m);
425  t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
426  t_approx_P_adjont_log_du_domega(n, L) = 0;
427 
428  }
429 
430  break;
431 
432  case SMALL_ROT:
433 
434  {
435 
436  // assume small current rotation
438  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
439  t_h(i, j) = t_Omega(i, j) + t_u(i, j);
440  t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, j, k);
441  t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
442 
443  // conservation of angular momentum at reference condition
444  t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
445  t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
446  t_levi_kirchoff_domega(k, l) = 0;
447 
448  t_approx_P_adjont_dstretch(i, j) = t_approx_P(i, j);
449  t_approx_P_adjont_log_du(L) =
450  t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
451  t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
452  t_approx_P_adjont_log_du_domega(m, L) = 0;
453 
454  }
455 
456  break;
457  }
458 
460  t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
461  eigen_vec(i, j) = t_C_h1(i, j);
462  CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
463  switch (EshelbianCore::stretchSelector) {
465  break;
467  // for (int ii = 0; ii != 3; ++ii) {
468  // eig(ii) = std::max(
469  // eig(ii), 2 * EshelbianCore::inv_f(EshelbianCore::f(
470  // std::numeric_limits<double>::min_exponent)));
471  // }
472  break;
473  }
474 
475  switch (EshelbianCore::gradApproximator) {
476  case NO_H1_CONFIGURATION:
477  t_log_u2_h1(i, j) = 0;
478  t_log_stretch_total(i, j) = t_log_u(i, j);
479  break;
480  case LARGE_ROT:
481  case MODERATE_ROT: {
482  auto t_log_u2_h1_tmp =
483  EigenMatrix::getMat(eig, eigen_vec, EshelbianCore::inv_f);
484  t_log_u2_h1(i, j) = t_log_u2_h1_tmp(i, j);
485  t_log_stretch_total(i, j) = t_log_u2_h1_tmp(i, j) / 2 + t_log_u(i, j);
486  } break;
487  case SMALL_ROT:
488  t_log_u2_h1(i, j) = 0;
489  t_log_stretch_total(i, j) = t_log_u(i, j);
490  break;
491  };
492 
493  } else {
494 
496 
498  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
499  t_u(i, j) = t_kd(i, j) + t_log_u(i, j);
500  t_h(i, j) = t_Omega(i, j) + t_u(i, j);
501  t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, j, k);
502 
503  // conservation of angular momentum at reference condition
504  t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
505  t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
506  t_levi_kirchoff_domega(k, l) = 0;
507 
508  t_log_stretch_total(i, j) = t_log_u(i, j);
509  }
510 
511  ++t_h;
512  ++t_h_domega;
513  ++t_h_dlog_u;
514  ++t_levi_kirchoff;
515  ++t_levi_kirchoff_dP;
516  ++t_levi_kirchoff_domega;
517  ++t_approx_P_adjont_dstretch;
518  ++t_approx_P_adjont_log_du;
519  ++t_approx_P_adjont_log_du_dP;
520  ++t_approx_P_adjont_log_du_domega;
521  ++t_approx_P;
522  ++t_R;
523  ++t_diff_R;
524  ++t_log_u;
525  ++t_u;
526  ++t_diff_u;
527  ++t_eigen_vals;
528  ++t_eigen_vecs;
529  ++t_omega;
530  ++t_grad_h1;
531  ++t_log_u2_h1;
532  ++t_log_stretch_total;
533  }
534 
536 }
537 
538 MoFEMErrorCode OpSpatialEquilibrium::integrate(EntData &data) {
540  int nb_dofs = data.getIndices().size();
541  int nb_integration_pts = data.getN().size1();
542  auto v = getVolume();
543  auto t_w = getFTensor0IntegrationWeight();
544  auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
545  auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
546  if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
547  dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
548  dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
549  dataAtPts->wL2DotDotAtPts.clear();
550  }
551  auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
552 
553  int nb_base_functions = data.getN().size2();
554  auto t_row_base_fun = data.getFTensor0N();
555 
557  auto get_ftensor1 = [](auto &v) {
559  &v[2]);
560  };
561 
562  for (int gg = 0; gg != nb_integration_pts; ++gg) {
563  double a = v * t_w;
564  auto t_nf = get_ftensor1(nF);
565  int bb = 0;
566  for (; bb != nb_dofs / 3; ++bb) {
567  t_nf(i) += a * t_row_base_fun * t_div_P(i);
568  t_nf(i) -= a * t_row_base_fun * alphaW * t_s_dot_w(i);
569  t_nf(i) -= a * t_row_base_fun * alphaRho * t_s_dot_dot_w(i);
570  ++t_nf;
571  ++t_row_base_fun;
572  }
573  for (; bb != nb_base_functions; ++bb)
574  ++t_row_base_fun;
575  ++t_w;
576  ++t_div_P;
577  ++t_s_dot_w;
578  ++t_s_dot_dot_w;
579  }
580 
582 }
583 
584 MoFEMErrorCode OpSpatialRotation::integrate(EntData &data) {
586  int nb_dofs = data.getIndices().size();
587  int nb_integration_pts = getGaussPts().size2();
588  auto v = getVolume();
589  auto t_w = getFTensor0IntegrationWeight();
590  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
591  int nb_base_functions = data.getN().size2();
592  auto t_row_base_fun = data.getFTensor0N();
596  auto get_ftensor1 = [](auto &v) {
598  &v[2]);
599  };
600  for (int gg = 0; gg != nb_integration_pts; ++gg) {
601  double a = v * t_w;
602  auto t_nf = get_ftensor1(nF);
603  int bb = 0;
604  for (; bb != nb_dofs / 3; ++bb) {
605  t_nf(k) += (a * t_row_base_fun) * t_levi_kirchoff(k);
606  ++t_nf;
607  ++t_row_base_fun;
608  }
609  for (; bb != nb_base_functions; ++bb) {
610  ++t_row_base_fun;
611  }
612  ++t_w;
613  ++t_levi_kirchoff;
614  }
616 }
617 
618 MoFEMErrorCode OpSpatialConsistencyP::integrate(EntData &data) {
620  int nb_dofs = data.getIndices().size();
621  int nb_integration_pts = data.getN().size1();
622  auto v = getVolume();
623  auto t_w = getFTensor0IntegrationWeight();
624  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
625 
626  int nb_base_functions = data.getN().size2() / 3;
627  auto t_row_base_fun = data.getFTensor1N<3>();
632 
633  auto get_ftensor1 = [](auto &v) {
635  &v[2]);
636  };
637 
638  for (int gg = 0; gg != nb_integration_pts; ++gg) {
639  double a = v * t_w;
640  auto t_nf = get_ftensor1(nF);
641 
642  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
644  t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
645 
646  int bb = 0;
647  for (; bb != nb_dofs / 3; ++bb) {
648  t_nf(i) += a * t_row_base_fun(j) * t_residuum(i, j);
649  ++t_nf;
650  ++t_row_base_fun;
651  }
652 
653  for (; bb != nb_base_functions; ++bb)
654  ++t_row_base_fun;
655 
656  ++t_w;
657  ++t_h;
658  }
659 
661 }
662 
663 MoFEMErrorCode OpSpatialConsistencyBubble::integrate(EntData &data) {
665  int nb_dofs = data.getIndices().size();
666  int nb_integration_pts = data.getN().size1();
667  auto v = getVolume();
668  auto t_w = getFTensor0IntegrationWeight();
669  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
670 
671  int nb_base_functions = data.getN().size2() / 9;
672  auto t_row_base_fun = data.getFTensor2N<3, 3>();
677 
678  auto get_ftensor0 = [](auto &v) {
680  };
681 
682  for (int gg = 0; gg != nb_integration_pts; ++gg) {
683  double a = v * t_w;
684  auto t_nf = get_ftensor0(nF);
685 
686  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
688  t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
689 
690  int bb = 0;
691  for (; bb != nb_dofs; ++bb) {
692  t_nf += a * t_row_base_fun(i, j) * t_residuum(i, j);
693  ++t_nf;
694  ++t_row_base_fun;
695  }
696  for (; bb != nb_base_functions; ++bb) {
697  ++t_row_base_fun;
698  }
699  ++t_w;
700  ++t_h;
701  }
702 
704 }
705 
706 MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
708  int nb_dofs = data.getIndices().size();
709  int nb_integration_pts = data.getN().size1();
710  auto v = getVolume();
711  auto t_w = getFTensor0IntegrationWeight();
712  auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
713  int nb_base_functions = data.getN().size2() / 3;
714  auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
716  auto get_ftensor1 = [](auto &v) {
718  &v[2]);
719  };
720 
721  for (int gg = 0; gg != nb_integration_pts; ++gg) {
722  double a = v * t_w;
723  auto t_nf = get_ftensor1(nF);
724  int bb = 0;
725  for (; bb != nb_dofs / 3; ++bb) {
726  double div_row_base = t_row_diff_base_fun(i, i);
727  t_nf(i) += a * div_row_base * t_w_l2(i);
728  ++t_nf;
729  ++t_row_diff_base_fun;
730  }
731  for (; bb != nb_base_functions; ++bb) {
732  ++t_row_diff_base_fun;
733  }
734  ++t_w;
735  ++t_w_l2;
736  }
737 
739 }
740 
741 template <AssemblyType A>
744  // get entity of face
745  EntityHandle fe_ent = OP::getFEEntityHandle();
746  // iterate over all boundary data
747  for (auto &bc : (*bcDispPtr)) {
748  // check if finite element entity is part of boundary condition
749  if (bc.faces.find(fe_ent) != bc.faces.end()) {
750  int nb_dofs = data.getIndices().size();
751  int nb_integration_pts = OP::getGaussPts().size2();
752  auto t_normal = OP::getFTensor1NormalsAtGaussPts();
753  auto t_w = OP::getFTensor0IntegrationWeight();
754  int nb_base_functions = data.getN().size2() / 3;
755  auto t_row_base_fun = data.getFTensor1N<3>();
756 
759 
760  double scale = 1;
761  for (auto &sm : scalingMethodsVec) {
762  scale *= sm->getScale(OP::getFEMethod()->ts_t);
763  }
764 
765  // get bc data
766  FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
767  t_bc_disp(i) *= scale;
768 
769  for (int gg = 0; gg != nb_integration_pts; ++gg) {
770  auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
771  int bb = 0;
772  for (; bb != nb_dofs / SPACE_DIM; ++bb) {
773  t_nf(i) -=
774  t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_disp(i) * 0.5;
775  ++t_nf;
776  ++t_row_base_fun;
777  }
778  for (; bb != nb_base_functions; ++bb)
779  ++t_row_base_fun;
780 
781  ++t_w;
782  ++t_normal;
783  }
784  }
785  }
787 }
788 
789 MoFEMErrorCode OpDispBc::iNtegrate(EntData &data) {
790  return OP::iNtegrate(data);
791 }
792 
793 template <AssemblyType A>
796 
797  FTENSOR_INDEX(3, i);
798  FTENSOR_INDEX(3, j);
799  FTENSOR_INDEX(3, k);
800 
801  // get entity of face
802  EntityHandle fe_ent = OP::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 = OP::getGaussPts().size2();
809  auto t_normal = OP::getFTensor1NormalsAtGaussPts();
810  auto t_w = OP::getFTensor0IntegrationWeight();
811 
812  int nb_base_functions = data.getN().size2() / 3;
813  auto t_row_base_fun = data.getFTensor1N<3>();
814 
815  auto get_ftensor1 = [](auto &v) {
817  &v[2]);
818  };
819 
820  // get bc data
821  FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
822 
823  double theta = bc.theta;
824  theta *= OP::getFEMethod()->ts_t;
825 
827  const double a = sqrt(t_normal(i) * t_normal(i));
828  t_omega(i) = t_normal(i) * (theta / a);
829  auto t_R = get_rotation_form_vector(t_omega, LARGE_ROT);
830  auto t_coords = OP::getFTensor1CoordsAtGaussPts();
831 
832  for (int gg = 0; gg != nb_integration_pts; ++gg) {
834  t_delta(i) = t_center(i) - t_coords(i);
836  t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
837 
838  auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
839  int bb = 0;
840  for (; bb != nb_dofs / 3; ++bb) {
841  t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
842  ++t_nf;
843  ++t_row_base_fun;
844  }
845  for (; bb != nb_base_functions; ++bb)
846  ++t_row_base_fun;
847 
848  ++t_w;
849  ++t_normal;
850  ++t_coords;
851  }
852  }
853  }
855 }
856 
857 MoFEMErrorCode OpRotationBc::iNtegrate(EntData &data) {
858  return OP::iNtegrate(data);
859 }
860 
861 MoFEMErrorCode OpBrokenTractionBc::iNtegrate(EntData &data) {
863 
864  FTENSOR_INDEX(3, i);
865 
866  int nb_dofs = data.getFieldData().size();
867  int nb_integration_pts = getGaussPts().size2();
868  int nb_base_functions = data.getN().size2();
869  double ts_t = getFEMethod()->ts_t;
870 
871 #ifndef NDEBUG
872  if (this->locF.size() != nb_dofs)
873  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
874  "Size of locF %d != nb_dofs %d", this->locF.size(), nb_dofs);
875 #endif // NDEBUG
876 
877  auto integrate_rhs = [&](auto &bc) {
879 
880  auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
881  auto t_row_base = data.getFTensor0N();
882  auto t_w = getFTensor0IntegrationWeight();
883 
884  for (int gg = 0; gg != nb_integration_pts; ++gg) {
885  auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
886  int rr = 0;
887  for (; rr != nb_dofs / SPACE_DIM; ++rr) {
888  t_f(i) -= ts_t * t_w * t_row_base * t_val(i);
889  ++t_row_base;
890  ++t_f;
891  }
892  for (; rr != nb_base_functions; ++rr)
893  ++t_row_base;
894  ++t_w;
895  }
896 
897  this->locF *= getMeasure();
899  };
900 
901  auto integrate_rhs_cook = [&](auto &bc) {
903 
904  auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
905  auto t_row_base = data.getFTensor0N();
906  auto t_w = getFTensor0IntegrationWeight();
907  auto t_coords = getFTensor1CoordsAtGaussPts();
908 
909  for (int gg = 0; gg != nb_integration_pts; ++gg) {
910 
911  auto calc_tau = [](double y) {
912  y -= 44;
913  y /= (60 - 44);
914  return -y * (y - 1) / 0.25;
915  };
916 
917  const auto tau = calc_tau(t_coords(1));
918  auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
919  int rr = 0;
920  for (; rr != nb_dofs / SPACE_DIM; ++rr) {
921  t_f(i) -= ts_t * t_w * t_row_base * tau * t_val(i);
922  ++t_row_base;
923  ++t_f;
924  }
925 
926  for (; rr != nb_base_functions; ++rr)
927  ++t_row_base;
928  ++t_w;
929  ++t_coords;
930  }
931 
932  // Multiply by 2 since we integrate on triangle, and hybrid constrain is
933  // multiplied by norm.
934  this->locF *= 2. * getMeasure();
935 
937  };
938 
939  // get entity of face
940  EntityHandle fe_ent = getFEEntityHandle();
941  for (auto &bc : *(bcData)) {
942  if (bc.faces.find(fe_ent) != bc.faces.end()) {
943 
944  int nb_dofs = data.getFieldData().size();
945  if (nb_dofs) {
946 
947  if (std::regex_match(bc.blockName, std::regex(".*COOK.*")))
948  CHKERR integrate_rhs_cook(bc);
949  else
950  CHKERR integrate_rhs(bc);
951  }
952  }
953  }
954 
956 }
957 
958 MoFEMErrorCode OpSpatialEquilibrium_dw_dP::integrate(EntData &row_data,
959  EntData &col_data) {
961  int nb_integration_pts = row_data.getN().size1();
962  int row_nb_dofs = row_data.getIndices().size();
963  int col_nb_dofs = col_data.getIndices().size();
964  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
966  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
967  };
969  auto v = getVolume();
970  auto t_w = getFTensor0IntegrationWeight();
971  int row_nb_base_functions = row_data.getN().size2();
972  auto t_row_base_fun = row_data.getFTensor0N();
973  for (int gg = 0; gg != nb_integration_pts; ++gg) {
974  double a = v * t_w;
975  int rr = 0;
976  for (; rr != row_nb_dofs / 3; ++rr) {
977  auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
978  auto t_m = get_ftensor1(K, 3 * rr, 0);
979  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
980  double div_col_base = t_col_diff_base_fun(i, i);
981  t_m(i) += a * t_row_base_fun * div_col_base;
982  ++t_m;
983  ++t_col_diff_base_fun;
984  }
985  ++t_row_base_fun;
986  }
987  for (; rr != row_nb_base_functions; ++rr)
988  ++t_row_base_fun;
989  ++t_w;
990  }
992 }
993 
994 MoFEMErrorCode OpSpatialEquilibrium_dw_dw::integrate(EntData &row_data,
995  EntData &col_data) {
997 
998  if (alphaW < std::numeric_limits<double>::epsilon() &&
999  alphaRho < std::numeric_limits<double>::epsilon())
1001 
1002  const int nb_integration_pts = row_data.getN().size1();
1003  const int row_nb_dofs = row_data.getIndices().size();
1004  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1006  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
1007 
1008  );
1009  };
1011 
1012  auto v = getVolume();
1013  auto t_w = getFTensor0IntegrationWeight();
1014 
1015  int row_nb_base_functions = row_data.getN().size2();
1016  auto t_row_base_fun = row_data.getFTensor0N();
1017 
1018  double ts_scale = alphaW * getTSa();
1019  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
1020  ts_scale += alphaRho * getTSaa();
1021 
1022  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1023  double a = v * t_w * ts_scale;
1024 
1025  int rr = 0;
1026  for (; rr != row_nb_dofs / 3; ++rr) {
1027 
1028  auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
1029  auto t_m = get_ftensor1(K, 3 * rr, 0);
1030  for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1031  const double b = a * t_row_base_fun * t_col_base_fun;
1032  t_m(i) -= b;
1033  ++t_m;
1034  ++t_col_base_fun;
1035  }
1036 
1037  ++t_row_base_fun;
1038  }
1039 
1040  for (; rr != row_nb_base_functions; ++rr)
1041  ++t_row_base_fun;
1042 
1043  ++t_w;
1044  }
1045 
1047 }
1048 
1049 MoFEMErrorCode OpSpatialPhysical_du_dP::integrate(EntData &row_data,
1050  EntData &col_data) {
1052 
1054 
1055  int nb_integration_pts = row_data.getN().size1();
1056  int row_nb_dofs = row_data.getIndices().size();
1057  int col_nb_dofs = col_data.getIndices().size();
1058  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1060 
1061  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1062 
1063  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1064 
1065  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1066 
1067  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1068 
1069  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1070 
1071  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
1072  };
1073 
1076 
1077  auto v = getVolume();
1078  auto t_w = getFTensor0IntegrationWeight();
1079  auto t_approx_P_adjont_log_du_dP =
1080  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1081  int row_nb_base_functions = row_data.getN().size2();
1082  auto t_row_base_fun = row_data.getFTensor0N();
1083  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1084  double a = v * t_w;
1085  int rr = 0;
1086  for (; rr != row_nb_dofs / 6; ++rr) {
1087 
1088  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1089  auto t_m = get_ftensor3(K, 6 * rr, 0);
1090 
1091  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1092  t_m(L, i) +=
1093  a * (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(j)) *
1094  t_row_base_fun;
1095  ++t_col_base_fun;
1096  ++t_m;
1097  }
1098 
1099  ++t_row_base_fun;
1100  }
1101  for (; rr != row_nb_base_functions; ++rr)
1102  ++t_row_base_fun;
1103  ++t_w;
1104  ++t_approx_P_adjont_log_du_dP;
1105  }
1107 }
1108 
1109 MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
1110  EntData &col_data) {
1112 
1114 
1115  int nb_integration_pts = row_data.getN().size1();
1116  int row_nb_dofs = row_data.getIndices().size();
1117  int col_nb_dofs = col_data.getIndices().size();
1118  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1120  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
1121  &m(r + 5, c));
1122  };
1123 
1126 
1127  auto v = getVolume();
1128  auto t_w = getFTensor0IntegrationWeight();
1129  auto t_approx_P_adjont_log_du_dP =
1130  getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1131 
1132  int row_nb_base_functions = row_data.getN().size2();
1133  auto t_row_base_fun = row_data.getFTensor0N();
1134  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1135  double a = v * t_w;
1136  int rr = 0;
1137  for (; rr != row_nb_dofs / 6; ++rr) {
1138  auto t_m = get_ftensor2(K, 6 * rr, 0);
1139  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1140  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1141  t_m(L) +=
1142  a * (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(i, j)) *
1143  t_row_base_fun;
1144  ++t_m;
1145  ++t_col_base_fun;
1146  }
1147  ++t_row_base_fun;
1148  }
1149  for (; rr != row_nb_base_functions; ++rr)
1150  ++t_row_base_fun;
1151  ++t_w;
1152  ++t_approx_P_adjont_log_du_dP;
1153  }
1155 }
1156 
1157 MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
1158  EntData &col_data) {
1160 
1162  auto t_L = symm_L_tensor();
1163 
1164  int row_nb_dofs = row_data.getIndices().size();
1165  int col_nb_dofs = col_data.getIndices().size();
1166  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1168 
1169  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1170 
1171  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1172 
1173  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1174 
1175  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1176 
1177  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1178 
1179  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1180 
1181  );
1182  };
1188 
1189  auto v = getVolume();
1190  auto t_w = getFTensor0IntegrationWeight();
1191  auto t_approx_P_adjont_log_du_domega =
1192  getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
1193 
1194  int row_nb_base_functions = row_data.getN().size2();
1195  auto t_row_base_fun = row_data.getFTensor0N();
1196 
1197  int nb_integration_pts = row_data.getN().size1();
1198 
1199  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1200  double a = v * t_w;
1201 
1202  int rr = 0;
1203  for (; rr != row_nb_dofs / 6; ++rr) {
1204  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1205  auto t_m = get_ftensor3(K, 6 * rr, 0);
1206  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1207  double v = a * t_row_base_fun * t_col_base_fun;
1208  t_m(L, k) += v * t_approx_P_adjont_log_du_domega(k, L);
1209  ++t_m;
1210  ++t_col_base_fun;
1211  }
1212  ++t_row_base_fun;
1213  }
1214 
1215  for (; rr != row_nb_base_functions; ++rr)
1216  ++t_row_base_fun;
1217 
1218  ++t_w;
1219  ++t_approx_P_adjont_log_du_domega;
1220  }
1221 
1223 }
1224 
1225 MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
1226  EntData &col_data) {
1228  int nb_integration_pts = getGaussPts().size2();
1229  int row_nb_dofs = row_data.getIndices().size();
1230  int col_nb_dofs = col_data.getIndices().size();
1231  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1233 
1234  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1235 
1236  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1237 
1238  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1239 
1240  );
1241  };
1248  auto v = getVolume();
1249  auto t_w = getFTensor0IntegrationWeight();
1250  auto t_levi_kirchoff_dP =
1251  getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
1252 
1253  int row_nb_base_functions = row_data.getN().size2();
1254  auto t_row_base_fun = row_data.getFTensor0N();
1255 
1256  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1257  double a = v * t_w;
1258  int rr = 0;
1259  for (; rr != row_nb_dofs / 3; ++rr) {
1260  double b = a * t_row_base_fun;
1261  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1262  auto t_m = get_ftensor2(K, 3 * rr, 0);
1263  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1264  t_m(k, i) += b * (t_levi_kirchoff_dP(i, l, k) * t_col_base_fun(l));
1265  ++t_m;
1266  ++t_col_base_fun;
1267  }
1268  ++t_row_base_fun;
1269  }
1270  for (; rr != row_nb_base_functions; ++rr) {
1271  ++t_row_base_fun;
1272  }
1273 
1274  ++t_w;
1275  ++t_levi_kirchoff_dP;
1276  }
1278 }
1279 
1280 MoFEMErrorCode OpSpatialRotation_domega_dBubble::integrate(EntData &row_data,
1281  EntData &col_data) {
1283  int nb_integration_pts = getGaussPts().size2();
1284  int row_nb_dofs = row_data.getIndices().size();
1285  int col_nb_dofs = col_data.getIndices().size();
1286 
1287  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1289  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
1290  };
1291 
1297 
1298  auto v = getVolume();
1299  auto t_w = getFTensor0IntegrationWeight();
1300  auto t_levi_kirchoff_dP =
1301  getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
1302 
1303  int row_nb_base_functions = row_data.getN().size2();
1304  auto t_row_base_fun = row_data.getFTensor0N();
1305 
1306  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1307  double a = v * t_w;
1308  int rr = 0;
1309  for (; rr != row_nb_dofs / 3; ++rr) {
1310  double b = a * t_row_base_fun;
1311  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1312  auto t_m = get_ftensor1(K, 3 * rr, 0);
1313  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1314  t_m(k) += b * (t_levi_kirchoff_dP(i, j, k) * t_col_base_fun(i, j));
1315  ++t_m;
1316  ++t_col_base_fun;
1317  }
1318  ++t_row_base_fun;
1319  }
1320 
1321  for (; rr != row_nb_base_functions; ++rr) {
1322  ++t_row_base_fun;
1323  }
1324  ++t_w;
1325  ++t_levi_kirchoff_dP;
1326  }
1328 }
1329 
1330 MoFEMErrorCode OpSpatialRotation_domega_domega::integrate(EntData &row_data,
1331  EntData &col_data) {
1333  int nb_integration_pts = getGaussPts().size2();
1334  int row_nb_dofs = row_data.getIndices().size();
1335  int col_nb_dofs = col_data.getIndices().size();
1336  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1338 
1339  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1340 
1341  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1342 
1343  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1344 
1345  );
1346  };
1353 
1354  auto v = getVolume();
1355  auto t_w = getFTensor0IntegrationWeight();
1356  auto t_levi_kirchoff_domega =
1357  getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
1358  int row_nb_base_functions = row_data.getN().size2();
1359  auto t_row_base_fun = row_data.getFTensor0N();
1360  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1361  double a = v * t_w;
1362  int rr = 0;
1363  for (; rr != row_nb_dofs / 3; ++rr) {
1364  auto t_m = get_ftensor2(K, 3 * rr, 0);
1365  const double b = a * t_row_base_fun;
1366  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1367  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1368  t_m(k, l) += (b * t_col_base_fun) * t_levi_kirchoff_domega(k, l);
1369  ++t_m;
1370  ++t_col_base_fun;
1371  }
1372  ++t_row_base_fun;
1373  }
1374  for (; rr != row_nb_base_functions; ++rr) {
1375  ++t_row_base_fun;
1376  }
1377  ++t_w;
1378  ++t_levi_kirchoff_domega;
1379  }
1381 }
1382 
1383 MoFEMErrorCode OpSpatialConsistency_dP_dP::integrate(EntData &row_data,
1384  EntData &col_data) {
1386  if (dataAtPts->matInvD.size1() == size_symm &&
1387  dataAtPts->matInvD.size2() == size_symm) {
1388  return integrateImpl<0>(row_data, col_data);
1389  } else {
1390  return integrateImpl<size_symm>(row_data, col_data);
1391  }
1393 };
1394 
1395 template <int S>
1396 MoFEMErrorCode OpSpatialConsistency_dP_dP::integrateImpl(EntData &row_data,
1397  EntData &col_data) {
1399 
1400  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1402 
1403  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1404 
1405  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1406 
1407  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1408 
1409  );
1410  };
1411 
1412  int nb_integration_pts = getGaussPts().size2();
1413  int row_nb_dofs = row_data.getIndices().size();
1414  int col_nb_dofs = col_data.getIndices().size();
1415 
1416  auto v = getVolume();
1417  auto t_w = getFTensor0IntegrationWeight();
1418  int row_nb_base_functions = row_data.getN().size2() / 3;
1419 
1424 
1425  auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
1426  &*dataAtPts->matInvD.data().begin());
1427 
1428  auto t_row_base = row_data.getFTensor1N<3>();
1429  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1430  double a = v * t_w;
1431 
1432  int rr = 0;
1433  for (; rr != row_nb_dofs / 3; ++rr) {
1434  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1435  auto t_m = get_ftensor2(K, 3 * rr, 0);
1436  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1437  t_m(i, k) += a * t_row_base(j) * (t_inv_D(i, j, k, l) * t_col_base(l));
1438  ++t_m;
1439  ++t_col_base;
1440  }
1441 
1442  ++t_row_base;
1443  }
1444 
1445  for (; rr != row_nb_base_functions; ++rr)
1446  ++t_row_base;
1447 
1448  ++t_w;
1449  ++t_inv_D;
1450  }
1452 }
1453 
1455 OpSpatialConsistency_dBubble_dBubble::integrate(EntData &row_data,
1456  EntData &col_data) {
1458  if (dataAtPts->matInvD.size1() == size_symm &&
1459  dataAtPts->matInvD.size2() == size_symm) {
1460  return integrateImpl<0>(row_data, col_data);
1461  } else {
1462  return integrateImpl<size_symm>(row_data, col_data);
1463  }
1465 };
1466 
1467 template <int S>
1469 OpSpatialConsistency_dBubble_dBubble::integrateImpl(EntData &row_data,
1470  EntData &col_data) {
1472 
1473  int nb_integration_pts = getGaussPts().size2();
1474  int row_nb_dofs = row_data.getIndices().size();
1475  int col_nb_dofs = col_data.getIndices().size();
1476 
1477  auto v = getVolume();
1478  auto t_w = getFTensor0IntegrationWeight();
1479  int row_nb_base_functions = row_data.getN().size2() / 9;
1480 
1485 
1486  auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
1487  &*dataAtPts->matInvD.data().begin());
1488 
1489  auto t_row_base = row_data.getFTensor2N<3, 3>();
1490  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1491  double a = v * t_w;
1492 
1493  int rr = 0;
1494  for (; rr != row_nb_dofs; ++rr) {
1495  auto t_col_base = col_data.getFTensor2N<3, 3>(gg, 0);
1496  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1497  K(rr, cc) +=
1498  a * (t_row_base(i, j) * (t_inv_D(i, j, k, l) * t_col_base(k, l)));
1499  ++t_col_base;
1500  }
1501 
1502  ++t_row_base;
1503  }
1504 
1505  for (; rr != row_nb_base_functions; ++rr)
1506  ++t_row_base;
1507  ++t_w;
1508  ++t_inv_D;
1509  }
1511 }
1512 
1513 MoFEMErrorCode OpSpatialConsistency_dBubble_dP::integrate(EntData &row_data,
1514  EntData &col_data) {
1516  if (dataAtPts->matInvD.size1() == size_symm &&
1517  dataAtPts->matInvD.size2() == size_symm) {
1518  return integrateImpl<0>(row_data, col_data);
1519  } else {
1520  return integrateImpl<size_symm>(row_data, col_data);
1521  }
1523 };
1524 
1525 template <int S>
1526 MoFEMErrorCode OpSpatialConsistency_dBubble_dP::integrateImpl(EntData &row_data,
1527  EntData &col_data) {
1529 
1530  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1532 
1533  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2)
1534 
1535  );
1536  };
1537 
1538  int nb_integration_pts = getGaussPts().size2();
1539  int row_nb_dofs = row_data.getIndices().size();
1540  int col_nb_dofs = col_data.getIndices().size();
1541 
1542  auto v = getVolume();
1543  auto t_w = getFTensor0IntegrationWeight();
1544  int row_nb_base_functions = row_data.getN().size2() / 9;
1545 
1552 
1553  auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
1554  &*dataAtPts->matInvD.data().begin());
1555 
1556  auto t_row_base = row_data.getFTensor2N<3, 3>();
1557  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1558  double a = v * t_w;
1559 
1560  auto t_m = get_ftensor1(K, 0, 0);
1561 
1562  int rr = 0;
1563  for (; rr != row_nb_dofs; ++rr) {
1564  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1565  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1566  t_m(k) += a * (t_row_base(i, j) * t_inv_D(i, j, k, l)) * t_col_base(l);
1567  ++t_col_base;
1568  ++t_m;
1569  }
1570 
1571  ++t_row_base;
1572  }
1573 
1574  for (; rr != row_nb_base_functions; ++rr)
1575  ++t_row_base;
1576  ++t_w;
1577  ++t_inv_D;
1578  }
1580 }
1581 
1582 MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
1583  EntData &col_data) {
1585 
1586  int nb_integration_pts = row_data.getN().size1();
1587  int row_nb_dofs = row_data.getIndices().size();
1588  int col_nb_dofs = col_data.getIndices().size();
1589 
1590  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1592 
1593  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1594 
1595  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1596 
1597  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1598 
1599  );
1600  };
1601 
1608 
1609  auto v = getVolume();
1610  auto t_w = getFTensor0IntegrationWeight();
1611  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
1612  int row_nb_base_functions = row_data.getN().size2() / 3;
1613  auto t_row_base_fun = row_data.getFTensor1N<3>();
1614 
1615  const double ts_a = getTSa();
1616 
1617  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1618  double a = v * t_w;
1619 
1620  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
1621 
1622  int rr = 0;
1623  for (; rr != row_nb_dofs / 3; ++rr) {
1624 
1626  t_PRT(i, k) = t_row_base_fun(j) * t_h_domega(i, j, k);
1627 
1628  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1629  auto t_m = get_ftensor2(K, 3 * rr, 0);
1630  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1631  t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
1632  ++t_m;
1633  ++t_col_base_fun;
1634  }
1635 
1636  ++t_row_base_fun;
1637  }
1638 
1639  for (; rr != row_nb_base_functions; ++rr)
1640  ++t_row_base_fun;
1641  ++t_w;
1642  ++t_h_domega;
1643  }
1645 }
1646 
1648 OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
1649  EntData &col_data) {
1651 
1652  int nb_integration_pts = row_data.getN().size1();
1653  int row_nb_dofs = row_data.getIndices().size();
1654  int col_nb_dofs = col_data.getIndices().size();
1655 
1656  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1658  &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
1659  };
1660 
1667 
1668  auto v = getVolume();
1669  auto t_w = getFTensor0IntegrationWeight();
1670  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
1671  int row_nb_base_functions = row_data.getN().size2() / 9;
1672  auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
1673 
1674  const double ts_a = getTSa();
1675 
1676  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1677  double a = v * t_w;
1678 
1679  int rr = 0;
1680  for (; rr != row_nb_dofs; ++rr) {
1681 
1683  t_PRT(k) = t_row_base_fun(i, j) * t_h_domega(i, j, k);
1684 
1685  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1686  auto t_m = get_ftensor2(K, rr, 0);
1687  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1688  t_m(j) += (a * t_col_base_fun) * t_PRT(j);
1689  ++t_m;
1690  ++t_col_base_fun;
1691  }
1692 
1693  ++t_row_base_fun;
1694  }
1695 
1696  for (; rr != row_nb_base_functions; ++rr)
1697  ++t_row_base_fun;
1698 
1699  ++t_w;
1700  ++t_h_domega;
1701  }
1703 }
1704 
1705 MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
1706  EntData &data) {
1708 
1709  if(tagSense != getSkeletonSense())
1711 
1712  auto create_tag = [this](const std::string tag_name, const int size) {
1713  double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1714  Tag th;
1715  CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
1716  th, MB_TAG_CREAT | MB_TAG_SPARSE,
1717  def_VAL);
1718  return th;
1719  };
1720 
1721  Tag th_w = create_tag("SpatialDisplacement", 3);
1722  Tag th_omega = create_tag("Omega", 3);
1723  Tag th_approxP = create_tag("Piola1Stress", 9);
1724  Tag th_calcSigma = create_tag("EshelbyStress", 9);
1725  Tag th_sigma = create_tag("CauchyStress", 9);
1726  Tag th_log_u = create_tag("LogSpatialStretch", 9);
1727  Tag th_u = create_tag("SpatialStretch", 9);
1728  Tag th_h = create_tag("h", 9);
1729  Tag th_X = create_tag("X", 3);
1730  Tag th_detF = create_tag("detF", 1);
1731  Tag th_angular_momentum = create_tag("AngularMomentum", 3);
1732 
1733  // Tag th_u_eig_vec = create_tag("SpatialStretchEigenVec", 9);
1734  // Tag th_u_eig_vals = create_tag("SpatialStretchEigenVals", 3);
1735  Tag th_traction = create_tag("traction", 3);
1736 
1737  Tag th_disp = create_tag("H1Displacement", 3);
1738  Tag th_disp_error = create_tag("DisplacementError", 1);
1739  Tag th_lambda_disp = create_tag("ContactDisplacement", 3);
1740 
1741  Tag th_energy = create_tag("Energy", 1);
1742 
1743  auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1744  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
1745  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1746  auto t_log_u =
1747  getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
1748  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1749  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1750  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1751  if (dataAtPts->SigmaAtPts.size2() != dataAtPts->approxPAtPts.size2()) {
1752  // that is for case that Eshelby stress is not calculated
1753  dataAtPts->SigmaAtPts.resize(dataAtPts->approxPAtPts.size1(),
1754  dataAtPts->approxPAtPts.size2(), false);
1755  dataAtPts->SigmaAtPts.clear();
1756  }
1757 
1758  auto t_calcSigma_P = getFTensor2FromMat<3, 3>(dataAtPts->SigmaAtPts);
1759  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
1760  auto t_coords = getFTensor1FromMat<3>(dataAtPts->XH1AtPts);
1761  auto t_normal = getFTensor1NormalsAtGaussPts();
1762  auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
1763  auto t_lambda_disp = getFTensor1FromMat<3>(dataAtPts->contactL2AtPts);
1764 
1765  if (dataAtPts->energyAtPts.size() == 0) {
1766  // that is for case that energy is not calculated
1767  dataAtPts->energyAtPts.resize(getGaussPts().size2());
1768  dataAtPts->energyAtPts.clear();
1769  }
1770  auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
1771 
1776 
1777  auto set_float_precision = [](const double x) {
1778  if (std::abs(x) < std::numeric_limits<float>::epsilon())
1779  return 0.;
1780  else
1781  return x;
1782  };
1783 
1784  // scalars
1785  auto save_scal_tag = [&](auto &th, auto v, const int gg) {
1787  v = set_float_precision(v);
1788  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, &v);
1790  };
1791 
1792  // vectors
1793  VectorDouble3 v(3);
1794  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
1795  auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
1797  t_v(i) = t_d(i);
1798  for (auto &a : v.data())
1799  a = set_float_precision(a);
1800  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1801  &*v.data().begin());
1803  };
1804 
1805  // tensors
1806 
1807  MatrixDouble3by3 m(3, 3);
1809  &m(0, 0), &m(0, 1), &m(0, 2),
1810 
1811  &m(1, 0), &m(1, 1), &m(1, 2),
1812 
1813  &m(2, 0), &m(2, 1), &m(2, 2));
1814 
1815  auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
1817  t_m(i, j) = t_d(i, j);
1818  for (auto &v : m.data())
1819  v = set_float_precision(v);
1820  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1821  &*m.data().begin());
1823  };
1824 
1825  const auto nb_gauss_pts = getGaussPts().size2();
1826  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1827 
1828  FTensor::Tensor1<double, 3> t_traction;
1829  t_traction(i) = t_approx_P(i, j) * t_normal(j) / t_normal.l2();
1830 
1831  // vectors
1832  CHKERR save_vec_tag(th_w, t_w, gg);
1833  CHKERR save_vec_tag(th_X, t_coords, gg);
1834  CHKERR save_vec_tag(th_omega, t_omega, gg);
1835  CHKERR save_vec_tag(th_traction, t_traction, gg);
1836 
1837  // tensors
1838  CHKERR save_mat_tag(th_h, t_h, gg);
1839 
1840  FTensor::Tensor2<double, 3, 3> t_log_u_tmp;
1841  for (int ii = 0; ii != 3; ++ii)
1842  for (int jj = 0; jj != 3; ++jj)
1843  t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
1844 
1845  CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
1846 
1848  for (int ii = 0; ii != 3; ++ii)
1849  for (int jj = 0; jj != 3; ++jj)
1850  t_u_tmp(ii, jj) = t_u(ii, jj);
1851 
1852  CHKERR save_mat_tag(th_u, t_u_tmp, gg);
1853  CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
1854  CHKERR save_mat_tag(th_calcSigma, t_calcSigma_P, gg);
1855  CHKERR save_vec_tag(th_disp, t_disp, gg);
1856  CHKERR save_vec_tag(th_lambda_disp, t_lambda_disp, gg);
1857 
1858  double u_error = sqrt((t_disp(i) - t_w(i)) * (t_disp(i) - t_w(i)));
1859  CHKERR save_scal_tag(th_disp_error, u_error, gg);
1860  CHKERR save_scal_tag(th_energy, t_energy, gg);
1861 
1862  const double jac = determinantTensor3by3(t_h);
1864  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
1865  CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
1866  CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
1867 
1869  t_levi(k) = t_levi_kirchoff(k);
1870  CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
1871  &t_levi(0));
1872 
1873  auto get_eiegn_vector_symmetric = [&](auto &t_u) {
1875 
1876  // for (int ii = 0; ii != 3; ++ii)
1877  // for (int jj = 0; jj != 3; ++jj)
1878  // t_m(ii, jj) = t_u(ii, jj);
1879 
1880  // VectorDouble3 eigen_values(3);
1881  // auto t_eigen_values = getFTensor1FromArray<3>(eigen_values);
1882  // if (computeEigenValuesSymmetric(t_m, t_eigen_values) != 0) {
1883  // };
1884 
1885  // CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
1886  // &*m.data().begin());
1887  // CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
1888  // &*eigen_values.data().begin());
1889 
1891  };
1892 
1893  CHKERR get_eiegn_vector_symmetric(t_u);
1894 
1895  ++t_w;
1896  ++t_h;
1897  ++t_log_u;
1898  ++t_u;
1899  ++t_omega;
1900  ++t_R;
1901  ++t_approx_P;
1902  ++t_calcSigma_P;
1903  ++t_levi_kirchoff;
1904  ++t_coords;
1905  ++t_normal;
1906  ++t_disp;
1907  ++t_lambda_disp;
1908  ++t_energy;
1909  }
1910 
1912 }
1913 
1915  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1916  std::vector<FieldSpace> spaces, std::string geom_field_name,
1917  boost::shared_ptr<Range> crack_front_edges_ptr) {
1919 
1920  {
1921 
1922  struct OpGetHONormalsOnFace
1923  : public MoFEM::OpGetHONormalsOnFace<SPACE_DIM> {
1924 
1926 
1927  OpGetHONormalsOnFace(std::string &field_name,
1928  boost::shared_ptr<Range> edges_ptr)
1929  : OP(field_name), edgesPtr(edges_ptr) {}
1930 
1931  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1932 
1933  auto ent = data.getFieldEntities().size()
1934  ? data.getFieldEntities()[0]->getEnt()
1935  : 0;
1936 
1937  if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
1938  return 0;
1939  } else {
1940  return OP::doWork(side, type, data);
1941  }
1942  };
1943 
1944  private:
1945  boost::shared_ptr<Range> edgesPtr;
1946  };
1947 
1948  auto jac = boost::make_shared<MatrixDouble>();
1949  auto det = boost::make_shared<VectorDouble>();
1950  pipeline.push_back(new OpGetHONormalsOnFace(
1951  geom_field_name, EshelbianCore::setSingularity
1952  ? crack_front_edges_ptr
1953  : boost::make_shared<Range>()));
1954  pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac));
1955  pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
1956  pipeline.push_back(new OpScaleBaseBySpaceInverseOfMeasure(
1957  L2, AINSWORTH_LEGENDRE_BASE, det));
1958  }
1959 
1960  CHKERR MoFEM::AddHOOps<2, 3, 3>::add(pipeline, spaces, geom_field_name);
1961 
1963 }
1964 
1966  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1967  std::vector<FieldSpace> spaces, std::string geom_field_name,
1968  boost::shared_ptr<Range> crack_front_edges_ptr) {
1970 
1971  {
1972  auto jac = boost::make_shared<MatrixDouble>();
1973  auto det = boost::make_shared<VectorDouble>();
1975  geom_field_name, jac));
1976  pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
1977  pipeline.push_back(
1979  }
1980 
1981  {
1982 
1984  : public MoFEM::OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM> {
1985 
1987 
1988  OpCalculateVectorFieldGradient(const std::string &field_name,
1989  boost::shared_ptr<MatrixDouble> jac,
1990  boost::shared_ptr<Range> edges_ptr)
1991  : OP(field_name, jac), edgesPtr(edges_ptr) {}
1992 
1993  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1994 
1995  auto ent = data.getFieldEntities().size()
1996  ? data.getFieldEntities()[0]->getEnt()
1997  : 0;
1998 
1999  if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
2000  return 0;
2001  } else {
2002  return OP::doWork(side, type, data);
2003  }
2004  };
2005 
2006  private:
2007  boost::shared_ptr<Range> edgesPtr;
2008  };
2009 
2010  auto jac = boost::make_shared<MatrixDouble>();
2011  auto det = boost::make_shared<VectorDouble>();
2012  pipeline.push_back(new OpCalculateVectorFieldGradient(
2013  geom_field_name, jac,
2014  EshelbianCore::setSingularity ? crack_front_edges_ptr
2015  : boost::make_shared<Range>()));
2016  pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
2017  pipeline.push_back(new OpScaleBaseBySpaceInverseOfMeasure(
2018  L2, AINSWORTH_LEGENDRE_BASE, det));
2019  }
2020 
2021  CHKERR MoFEM::AddHOOps<3, 3, 3>::add(pipeline, spaces, geom_field_name,
2022  nullptr, nullptr, nullptr);
2023 
2025 }
2026 
2027 } // namespace EshelbianPlasticity
EshelbianPlasticity::sort_eigen_vals
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
Definition: EshelbianOperators.cpp:159
EigenMatrix::getMat
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
Definition: MatrixFunction.cpp:53
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
EshelbianPlasticity::is_eq
auto is_eq(const double &a, const double &b)
Definition: EshelbianOperators.cpp:144
MoFEM::EntitiesFieldData::EntData::getFieldEntities
const VectorFieldEntities & getFieldEntities() const
get field entities
Definition: EntitiesFieldData.hpp:1277
EshelbianPlasticity::get_rotation_form_vector
auto get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:24
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:44
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
OpCalculateEshelbyStress::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::OpCalculateHOJacForFaceEmbeddedIn3DSpace
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
Definition: HODataOperators.hpp:265
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FTensor::Tensor1::l2
T l2() const
Definition: Tensor1_value.hpp:84
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
EshelbianPlasticity::AddHOOps
Definition: EshelbianPlasticity.hpp:1150
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
EshelbianPlasticity::isEq
Definition: EshelbianOperators.cpp:134
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1254
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
MoFEM.hpp
MoFEM::OpScaleBaseBySpaceInverseOfMeasure
Scale base functions by inverses of measure of element.
Definition: HODataOperators.hpp:390
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition: Templates.hpp:2011
EshelbianPlasticity::NO_H1_CONFIGURATION
@ NO_H1_CONFIGURATION
Definition: EshelbianPlasticity.hpp:43
scale
double scale
Definition: plastic.cpp:119
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
EshelbianPlasticity::isEq::absMax
static double absMax
Definition: EshelbianOperators.cpp:139
EshelbianPlasticity::get_uniq_nb
auto get_uniq_nb(double *ptr)
Definition: EshelbianOperators.cpp:148
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::EntitiesFieldData::EntData::getFTensor2DiffN
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Definition: EntitiesFieldData.cpp:660
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
a
constexpr double a
Definition: approx_sphere.cpp:30
convert.type
type
Definition: convert.py:64
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
lapack_wrap.h
MatrixFunction.hpp
MoFEM::EntitiesFieldData::EntData::getFTensor2N
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:762
EshelbianPlasticity::OpDispBcImpl
Definition: EshelbianPlasticity.hpp:612
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::OpGetHONormalsOnFace
Calculate normals at Gauss points of triangle element.
Definition: HODataOperators.hpp:281
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
convert.n
n
Definition: convert.py:82
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:43
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:43
FTensor::Tensor0
Definition: Tensor0.hpp:16
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:44
PlasticOps::symm_L_tensor
auto symm_L_tensor(FTensor::Number< DIM >)
Definition: PlasticOpsGeneric.hpp:21
EshelbianPlasticity::OpRotationBcImpl
Definition: EshelbianPlasticity.hpp:643
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
EshelbianPlasticity::get_diff2_rotation_form_vector
auto get_diff2_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:108
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::EntitiesFieldData::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: EntitiesFieldData.cpp:640
EshelbianPlasticity::RotSelector
RotSelector
Definition: EshelbianPlasticity.hpp:43
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EshelbianPlasticity::get_diff_rotation_form_vector
auto get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Definition: EshelbianOperators.cpp:60
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
MoFEM::Types::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
EshelbianAux.hpp
Auxilary functions for Eshelbian plasticity.
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
EshelbianPlasticity::isEq::check
static auto check(const double &a, const double &b)
Definition: EshelbianOperators.cpp:135
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:43
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EshelbianPlasticity::TensorTypeExtractor::Type
std::remove_pointer< T >::type Type
Definition: EshelbianAux.hpp:104
OP
EigenMatrix::getDiffMat
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
Definition: MatrixFunction.cpp:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1452