v0.10.0
EshelbianOperators.cpp
Go to the documentation of this file.
1 /**
2  * \file EshelbianOperators.cpp
3  * \example EshelbianOperators.cpp
4  *
5  * \brief Implementation of operators
6  */
7 
8 #include <MoFEM.hpp>
9 using namespace MoFEM;
10 
11 #include <BasicFiniteElements.hpp>
12 
13 #include <EshelbianPlasticity.hpp>
14 #include <boost/math/constants/constants.hpp>
15 
16 #include <lapack_wrap.h>
17 
18 #include <MatrixFunction.hpp>
19 
20 namespace EshelbianPlasticity {
21 
22 double f(const double v) { return exp(v); }
23 double d_f(const double v) { return exp(v); }
24 double dd_f(const double v) { return exp(v); }
25 
26 template <class T> struct TensorTypeExtractor {
28 };
29 template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
31 };
32 
33 template <typename T>
36  &t_omega) {
41  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
42  t_R(i, j) = t_kd(i, j);
43 
44  const typename TensorTypeExtractor<T>::Type angle =
45  sqrt(t_omega(i) * t_omega(i));
46 
47  constexpr double tol = 1e-18;
48  if (std::abs(angle) < tol) {
49  return t_R;
50  }
51 
53  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
54  const typename TensorTypeExtractor<T>::Type a = sin(angle) / angle;
55  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
56  const typename TensorTypeExtractor<T>::Type b =
57  2. * ss_2 * ss_2 / (angle * angle);
58  t_R(i, j) += a * t_Omega(i, j);
59  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
60 
61  return t_R;
62 }
63 
64 template <typename T>
67  &t_omega) {
68 
73 
74  const typename TensorTypeExtractor<T>::Type angle =
75  sqrt(t_omega(i) * t_omega(i));
76 
77  constexpr double tol = 1e-18;
78  if (std::abs(angle) < tol) {
80  auto t_R = get_rotation_form_vector(t_omega);
81  t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, l, k) * t_R(l, j);
82  return t_diff_R;
83  }
84 
87  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
88 
89  const typename TensorTypeExtractor<T>::Type angle2 = angle * angle;
90  const typename TensorTypeExtractor<T>::Type ss = sin(angle);
91  const typename TensorTypeExtractor<T>::Type a = ss / angle;
92  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
93  const typename TensorTypeExtractor<T>::Type b = 2. * ss_2 * ss_2 / angle2;
94  t_diff_R(i, j, k) = 0;
95  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
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 
100  const typename TensorTypeExtractor<T>::Type cc = cos(angle);
101  const typename TensorTypeExtractor<T>::Type cc_2 = cos(angle / 2.);
102  const typename TensorTypeExtractor<T>::Type diff_a =
103  (angle * cc - ss) / angle2;
104 
105  const typename TensorTypeExtractor<T>::Type diff_b =
106  (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
107  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
108  t_diff_R(i, j, k) +=
109  diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
110 
111  return t_diff_R;
112 }
113 
114 template <typename T1, typename T2, typename T3>
116  T3 &eig_vec) {
118  for (int ii = 0; ii != 3; ii++)
119  for (int jj = 0; jj != 3; jj++)
120  eig_vec(ii, jj) = X(ii, jj);
121 
122  int n = 3;
123  int lda = 3;
124  int lwork = (3 + 2) * 3;
125  std::array<double, (3 + 2) * 3> work;
126 
127  if (lapack_dsyev('V', 'U', n, &(eig_vec(0, 0)), lda, &eig(0), work.data(),
128  lwork) > 0)
129  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
130  "The algorithm failed to compute eigenvalues.");
132 }
133 
134 inline bool is_eq(const double &a, const double &b) {
135  constexpr double eps = 1e-8;
136  return abs(a - b) < eps;
137 }
138 
139 template <typename T> inline int get_uniq_nb(T &ec) {
140  return distance(&ec(0), unique(&ec(0), &ec(0) + 3, is_eq));
141 }
142 
143 template <typename T1, typename T2>
144 void sort_eigen_vals(T1 &eig, T2 &eigen_vec) {
147  if (is_eq(eig(0), eig(1))) {
148  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
149  eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
150  eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2),
151  eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2)};
152  FTensor::Tensor1<double, 3> eig_c{eig(0), eig(2), eig(1)};
153  eig(i) = eig_c(i);
154  eigen_vec(i, j) = eigen_vec_c(i, j);
155  } else {
156  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
157  eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2),
158  eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
159  eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2)};
160  FTensor::Tensor1<double, 3> eig_c{eig(1), eig(0), eig(2)};
161  eig(i) = eig_c(i);
162  eigen_vec(i, j) = eigen_vec_c(i, j);
163  }
164 }
165 
166 MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
167  EntityType type,
168  EntData &data) {
170  if (type != MBTET)
172  int nb_integration_pts = data.getN().size1();
177 
178  dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
179  dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
180  dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
181  dataAtPts->streachTensorAtPts.resize(6, nb_integration_pts, false);
182  dataAtPts->diffStreachTensorAtPts.resize(36, nb_integration_pts, false);
183  dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
184  dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
185  dataAtPts->nbUniq.resize(nb_integration_pts, false);
186 
187  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
188  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
189  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
190  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
191  auto t_log_u =
192  getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachTensorAtPts);
193  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
194 
195  auto t_diff_u =
196  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
197  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
198  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
199  auto &nbUniq = dataAtPts->nbUniq;
200 
201  for (int gg = 0; gg != nb_integration_pts; ++gg) {
202 
203  auto t0_diff = get_diff_rotation_form_vector(t_omega);
204  auto t0 = get_rotation_form_vector(t_omega);
205 
206  t_diff_R(i, j, k) = t0_diff(i, j, k);
207  t_R(i, j) = t0(i, j);
208 
211  CHKERR get_eigen_val_and_proj_lapack(t_log_u, eig, eigen_vec);
212 
213  t_eigen_vals(i) = eig(i);
214  t_eigen_vecs(i, j) = eigen_vec(i, j);
215 
216  // rare case when two eigen values are equal
217  nbUniq[gg] = get_uniq_nb(eig);
218  if (nbUniq[gg] == 2)
219  sort_eigen_vals(eig, eigen_vec);
220 
221  auto t_u_tmp = EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, f);
222  auto t_diff_u_tmp =
223  EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs, f, d_f, nbUniq[gg]);
224  for (int ii = 0; ii != 3; ++ii) {
225  for (int jj = ii; jj != 3; ++jj) {
226  for (int kk = 0; kk != 3; ++kk) {
227  for (int ll = 0; ll < kk; ++ll) {
228  t_diff_u_tmp(ii, jj, kk, ll) *= 2;
229  }
230  }
231  }
232  }
233 
234  t_u(i, j) = t_u_tmp(i, j);
235  t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
236  t_h(i, j) = t_R(i, k) * t_u(k, j);
237 
238  ++t_h;
239  ++t_R;
240  ++t_diff_R;
241  ++t_log_u;
242  ++t_u;
243  ++t_diff_u;
244  ++t_eigen_vals;
245  ++t_eigen_vecs;
246  ++t_omega;
247  }
248 
250 }
251 
252 MoFEMErrorCode OpSpatialEquilibrium::integrate(EntData &data) {
254  int nb_dofs = data.getIndices().size();
255  int nb_integration_pts = data.getN().size1();
256  auto v = getVolume();
257  auto t_w = getFTensor0IntegrationWeight();
258  auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
259  auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wDotAtPts);
260  if (dataAtPts->wDotDotAtPts.size1() == 0 &&
261  dataAtPts->wDotDotAtPts.size2() != nb_integration_pts) {
262  dataAtPts->wDotDotAtPts.resize(3, nb_integration_pts);
263  dataAtPts->wDotDotAtPts.clear();
264  }
265  auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wDotDotAtPts);
266 
267  int nb_base_functions = data.getN().size2();
268  auto t_row_base_fun = data.getFTensor0N();
269 
271  auto get_ftensor1 = [](auto &v) {
272  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
273  &v[2]);
274  };
275 
276  for (int gg = 0; gg != nb_integration_pts; ++gg) {
277  double a = v * t_w;
278  auto t_nf = get_ftensor1(nF);
279  int bb = 0;
280  for (; bb != nb_dofs / 3; ++bb) {
281  t_nf(i) += a * t_row_base_fun * t_div_P(i);
282  t_nf(i) -= a * t_row_base_fun * alphaW * t_s_dot_w(i);
283  t_nf(i) -= a * t_row_base_fun * alphaRho * t_s_dot_dot_w(i);
284  ++t_nf;
285  ++t_row_base_fun;
286  }
287  for (; bb != nb_base_functions; ++bb)
288  ++t_row_base_fun;
289  ++t_w;
290  ++t_div_P;
291  ++t_s_dot_w;
292  ++t_s_dot_dot_w;
293  }
294 
296 }
297 
298 MoFEMErrorCode OpSpatialRotation::integrate(EntData &data) {
300  int nb_dofs = data.getIndices().size();
301  int nb_integration_pts = data.getN().size1();
302  auto v = getVolume();
303  auto t_w = getFTensor0IntegrationWeight();
304  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
305  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
306 
307  int nb_base_functions = data.getN().size2();
308  auto t_row_base_fun = data.getFTensor0N();
314  auto get_ftensor1 = [](auto &v) {
315  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
316  &v[2]);
317  };
318 
319  for (int gg = 0; gg != nb_integration_pts; ++gg) {
320  double a = v * t_w;
321  auto t_nf = get_ftensor1(nF);
323  t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
324  FTensor::Tensor1<double, 3> t_leviPRT;
325  t_leviPRT(k) = levi_civita(i, m, k) * t_PRT(i, m);
326  int bb = 0;
327  for (; bb != nb_dofs / 3; ++bb) {
328  t_nf(k) += a * t_row_base_fun * t_leviPRT(k);
329  ++t_nf;
330  ++t_row_base_fun;
331  }
332  for (; bb != nb_base_functions; ++bb)
333  ++t_row_base_fun;
334  ++t_w;
335  ++t_approx_P;
336  ++t_R;
337  }
339 }
340 
341 MoFEMErrorCode OpSpatialPhysical::integrate(EntData &data) {
343  int nb_dofs = data.getIndices().size();
344  int nb_integration_pts = data.getN().size1();
345  auto v = getVolume();
346  auto t_w = getFTensor0IntegrationWeight();
347  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
348  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
349  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
350  auto t_dot_log_u =
351  getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachDotTensorAtPts);
352  auto t_diff_u =
353  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
354 
359  auto get_ftensor2 = [](auto &v) {
361  &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
362  };
363 
364  int nb_base_functions = data.getN().size2();
365  auto t_row_base_fun = data.getFTensor0N();
366  for (int gg = 0; gg != nb_integration_pts; ++gg) {
367  double a = v * t_w;
368  auto t_nf = get_ftensor2(nF);
369 
371  t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
373  t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
374 
375  FTensor::Tensor2_symmetric<double, 3> t_viscous_stress;
376  t_viscous_stress(i, j) = alphaU * t_dot_log_u(i, j);
377 
379  t1(i, j) =
380  a * (t_diff_u(k, l, i, j) * (t_deltaP(k, l) - t_viscous_stress(k, l)));
381 
382  int bb = 0;
383  for (; bb != nb_dofs / 6; ++bb) {
384 
385  t_nf(i, j) += t_row_base_fun * t1(i, j);
386 
387  ++t_nf;
388  ++t_row_base_fun;
389  }
390  for (; bb != nb_base_functions; ++bb)
391  ++t_row_base_fun;
392 
393  ++t_w;
394  ++t_approx_P;
395  ++t_P;
396  ++t_R;
397  ++t_dot_log_u;
398  ++t_diff_u;
399  }
401 }
402 
403 MoFEMErrorCode OpSpatialConsistencyP::integrate(EntData &data) {
405  int nb_dofs = data.getIndices().size();
406  int nb_integration_pts = data.getN().size1();
407  auto v = getVolume();
408  auto t_w = getFTensor0IntegrationWeight();
409  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
410  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
411  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
412 
413  int nb_base_functions = data.getN().size2() / 3;
414  auto t_row_base_fun = data.getFTensor1N<3>();
419 
420  auto get_ftensor1 = [](auto &v) {
421  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
422  &v[2]);
423  };
424 
425  for (int gg = 0; gg != nb_integration_pts; ++gg) {
426  double a = v * t_w;
427  auto t_nf = get_ftensor1(nF);
428 
429  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
431  t_residuum(i, j) = t_R(i, k) * t_u(k, j) - t_kd(i, j);
432  FTensor::Tensor2<double, 3, 3> t_residuum_omega;
433  t_residuum_omega(i, j) =
434  (levi_civita(i, m, k) * t_omega_dot(k)) * t_R(m, j);
435 
436  int bb = 0;
437  for (; bb != nb_dofs / 3; ++bb) {
438  t_nf(i) += a * t_row_base_fun(j) * t_residuum(i, j);
439  t_nf(i) += a * t_row_base_fun(j) * t_residuum_omega(i, j);
440 
441  ++t_nf;
442  ++t_row_base_fun;
443  }
444 
445  for (; bb != nb_base_functions; ++bb)
446  ++t_row_base_fun;
447 
448  ++t_w;
449  ++t_u;
450  ++t_R;
451  ++t_omega_dot;
452  }
453 
455 }
456 
457 MoFEMErrorCode OpSpatialConsistencyBubble::integrate(EntData &data) {
459  int nb_dofs = data.getIndices().size();
460  int nb_integration_pts = data.getN().size1();
461  auto v = getVolume();
462  auto t_w = getFTensor0IntegrationWeight();
463  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
464  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
465  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
466 
467  int nb_base_functions = data.getN().size2() / 9;
468  auto t_row_base_fun = data.getFTensor2N<3, 3>();
473 
474  auto get_ftensor0 = [](auto &v) {
476  };
477 
478  for (int gg = 0; gg != nb_integration_pts; ++gg) {
479  double a = v * t_w;
480  auto t_nf = get_ftensor0(nF);
481 
482  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
484  t_residuum(i, j) = t_R(i, k) * t_u(k, j) - t_kd(i, j);
485  FTensor::Tensor2<double, 3, 3> t_residuum_omega;
486  t_residuum_omega(i, j) = (levi_civita(i, m, k) * t_omega_dot(k)) * t_R(m, j);
487 
488  int bb = 0;
489  for (; bb != nb_dofs; ++bb) {
490  t_nf += a * t_row_base_fun(i, j) * t_residuum(i, j);
491  t_nf += a * t_row_base_fun(i, j) * t_residuum_omega(i, j);
492  ++t_nf;
493  ++t_row_base_fun;
494  }
495  for (; bb != nb_base_functions; ++bb) {
496  ++t_row_base_fun;
497  }
498  ++t_w;
499  ++t_R;
500  ++t_u;
501  ++t_omega_dot;
502  }
503 
505 }
506 
507 MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
509  int nb_dofs = data.getIndices().size();
510  int nb_integration_pts = data.getN().size1();
511  auto v = getVolume();
512  auto t_w = getFTensor0IntegrationWeight();
513  auto t_s_w = getFTensor1FromMat<3>(dataAtPts->wAtPts);
514  int nb_base_functions = data.getN().size2() / 3;
515  auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
517  auto get_ftensor1 = [](auto &v) {
518  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
519  &v[2]);
520  };
521 
522  for (int gg = 0; gg != nb_integration_pts; ++gg) {
523  double a = v * t_w;
524  auto t_nf = get_ftensor1(nF);
525  int bb = 0;
526  for (; bb != nb_dofs / 3; ++bb) {
527  double div_row_base = t_row_diff_base_fun(i, i);
528  t_nf(i) += a * div_row_base * t_s_w(i);
529  ++t_nf;
530  ++t_row_diff_base_fun;
531  }
532  for (; bb != nb_base_functions; ++bb) {
533  ++t_row_diff_base_fun;
534  }
535  ++t_w;
536  ++t_s_w;
537  }
538 
540 }
541 
542 MoFEMErrorCode OpDispBc::integrate(EntData &data) {
544  // get entity of face
545  EntityHandle fe_ent = getFEEntityHandle();
546  // interate over all boundary data
547  for (auto &bc : (*bcDispPtr)) {
548  // check if finite element entity is part of boundary condition
549  if (bc.faces.find(fe_ent) != bc.faces.end()) {
550  int nb_dofs = data.getIndices().size();
551  int nb_integration_pts = data.getN().size1();
552  auto t_normal = getFTensor1Normal();
553  auto t_w = getFTensor0IntegrationWeight();
554  int nb_base_functions = data.getN().size2() / 3;
555  auto t_row_base_fun = data.getFTensor1N<3>();
556 
559  auto get_ftensor1 = [](auto &v) {
560  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
561  &v[2]);
562  };
563 
564  // get bc data
565  FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
566  t_bc_disp(i) *= getFEMethod()->ts_t;
567 
568  for (int gg = 0; gg != nb_integration_pts; ++gg) {
570  t_bc_res(i) = t_bc_disp(i);
571  auto t_nf = get_ftensor1(nF);
572  int bb = 0;
573  for (; bb != nb_dofs / 3; ++bb) {
574  t_nf(i) -=
575  t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
576  ++t_nf;
577  ++t_row_base_fun;
578  }
579  for (; bb != nb_base_functions; ++bb)
580  ++t_row_base_fun;
581 
582  ++t_w;
583  }
584  }
585  }
587 }
588 
589 MoFEMErrorCode OpRotationBc::integrate(EntData &data) {
591  // get entity of face
592  EntityHandle fe_ent = getFEEntityHandle();
593  // interate over all boundary data
594  for (auto &bc : (*bcRotPtr)) {
595  // check if finite element entity is part of boundary condition
596  if (bc.faces.find(fe_ent) != bc.faces.end()) {
597  int nb_dofs = data.getIndices().size();
598  int nb_integration_pts = data.getN().size1();
599  auto t_normal = getFTensor1Normal();
600  auto t_w = getFTensor0IntegrationWeight();
601 
602  int nb_base_functions = data.getN().size2() / 3;
603  auto t_row_base_fun = data.getFTensor1N<3>();
607  auto get_ftensor1 = [](auto &v) {
608  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
609  &v[2]);
610  };
611 
612  // get bc data
613  FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
614 
615  double theta = bc.theta;
616  theta *= getFEMethod()->ts_t;
617 
619  const double a = sqrt(t_normal(i) * t_normal(i));
620  t_omega(i) = t_normal(i) * (theta / a);
621  auto t_R = get_rotation_form_vector(t_omega);
622  auto t_coords = getFTensor1CoordsAtGaussPts();
623 
624  for (int gg = 0; gg != nb_integration_pts; ++gg) {
626  t_delta(i) = t_center(i) - t_coords(i);
628  t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
629 
630  auto t_nf = get_ftensor1(nF);
631  int bb = 0;
632  for (; bb != nb_dofs / 3; ++bb) {
633  t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
634  ++t_nf;
635  ++t_row_base_fun;
636  }
637  for (; bb != nb_base_functions; ++bb)
638  ++t_row_base_fun;
639 
640  ++t_w;
641  ++t_coords;
642  }
643  }
644  }
646 }
647 
648 MoFEMErrorCode FeTractionBc::preProcess() {
650 
651  struct FaceRule {
652  int operator()(int p_row, int p_col, int p_data) const {
653  return 2 * (p_data + 1);
654  }
655  };
656 
657  if (ts_ctx == CTX_TSSETIFUNCTION) {
658 
659  // Loop boundary elements with traction boundary conditions
661  fe.getOpPtrVector().push_back(
662  new OpTractionBc(std::string("P") /* + "_RT"*/, *this));
663  fe.getRuleHook = FaceRule();
664  fe.ts_t = ts_t;
665  CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
666  fe);
667  }
668 
670 }
671 
672 MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type, EntData &data) {
674 
677 
678  auto t_normal = getFTensor1Normal();
679  const double nrm2 = sqrt(t_normal(i) * t_normal(i));
680  FTensor::Tensor1<double, 3> t_unit_normal;
681  t_unit_normal(i) = t_normal(i) / nrm2;
682  int nb_dofs = data.getFieldData().size();
683  int nb_integration_pts = data.getN().size1();
684  int nb_base_functions = data.getN().size2() / 3;
685  double ts_t = getFEMethod()->ts_t;
686 
687  auto integrate_matrix = [&]() {
689 
690  auto t_w = getFTensor0IntegrationWeight();
691  matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
692  matPP.clear();
693 
694  auto t_row_base_fun = data.getFTensor1N<3>();
695  for (int gg = 0; gg != nb_integration_pts; ++gg) {
696 
697  int rr = 0;
698  for (; rr != nb_dofs / 3; ++rr) {
699  const double a = t_row_base_fun(i) * t_unit_normal(i);
700  auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
701  for (int cc = 0; cc != nb_dofs / 3; ++cc) {
702  const double b = t_col_base_fun(i) * t_unit_normal(i);
703  matPP(rr, cc) += t_w * a * b;
704  ++t_col_base_fun;
705  }
706  ++t_row_base_fun;
707  }
708 
709  for (; rr != nb_base_functions; ++rr)
710  ++t_row_base_fun;
711 
712  ++t_w;
713  }
714 
716  };
717 
718  auto integrate_rhs = [&](auto &bc) {
720 
721  auto t_w = getFTensor0IntegrationWeight();
722  vecPv.resize(3, nb_dofs / 3, false);
723  vecPv.clear();
724 
725  auto t_row_base_fun = data.getFTensor1N<3>();
726  double ts_t = getFEMethod()->ts_t;
727 
728  for (int gg = 0; gg != nb_integration_pts; ++gg) {
729  int rr = 0;
730  for (; rr != nb_dofs / 3; ++rr) {
731  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
732  for (int dd = 0; dd != 3; ++dd)
733  if (bc.flags[dd])
734  vecPv(dd, rr) += t * bc.vals[dd];
735  ++t_row_base_fun;
736  }
737  for (; rr != nb_base_functions; ++rr)
738  ++t_row_base_fun;
739  ++t_w;
740  }
742  };
743 
744  auto integrate_rhs_cook = [&](auto &bc) {
746 
747  vecPv.resize(3, nb_dofs / 3, false);
748  vecPv.clear();
749 
750  auto t_w = getFTensor0IntegrationWeight();
751  auto t_row_base_fun = data.getFTensor1N<3>();
752  auto t_coords = getFTensor1CoordsAtGaussPts();
753 
754  for (int gg = 0; gg != nb_integration_pts; ++gg) {
755 
756  auto calc_tau = [](double y) {
757  y -= 44;
758  y /= (60 - 44);
759  return -y * (y - 1) / 0.25;
760  };
761 
762  const double tau = calc_tau(t_coords(1));
763 
764  int rr = 0;
765  for (; rr != nb_dofs / 3; ++rr) {
766  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
767  for (int dd = 0; dd != 3; ++dd)
768  if (bc.flags[dd])
769  vecPv(dd, rr) += t * tau * bc.vals[dd];
770  ++t_row_base_fun;
771  }
772 
773  for (; rr != nb_base_functions; ++rr)
774  ++t_row_base_fun;
775  ++t_w;
776  ++t_coords;
777  }
779  };
780 
781  // get entity of face
782  EntityHandle fe_ent = getFEEntityHandle();
783  for (auto &bc : *(bcFe.bcData)) {
784  if (bc.faces.find(fe_ent) != bc.faces.end()) {
785 
786  int nb_dofs = data.getFieldData().size();
787  if (nb_dofs) {
788 
789  CHKERR integrate_matrix();
790  if (bc.blockName.compare(20, 4, "COOK") == 0)
791  CHKERR integrate_rhs_cook(bc);
792  else
793  CHKERR integrate_rhs(bc);
794 
795  const auto info =
796  lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
797  nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
798  if (info > 0)
799  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
800  "The leading minor of order %i is "
801  "not positive; definite;\nthe "
802  "solution could not be computed",
803  info);
804 
805  for (int dd = 0; dd != 3; ++dd)
806  if (bc.flags[dd])
807  for (int rr = 0; rr != nb_dofs / 3; ++rr)
808  data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
809  }
810  }
811  }
812 
814 }
815 
816 MoFEMErrorCode OpSpatialEquilibrium_dw_dP::integrate(EntData &row_data,
817  EntData &col_data) {
819  int nb_integration_pts = row_data.getN().size1();
820  int row_nb_dofs = row_data.getIndices().size();
821  int col_nb_dofs = col_data.getIndices().size();
822  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
824  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
825  };
827  auto v = getVolume();
828  auto t_w = getFTensor0IntegrationWeight();
829  int row_nb_base_functions = row_data.getN().size2();
830  auto t_row_base_fun = row_data.getFTensor0N();
831  for (int gg = 0; gg != nb_integration_pts; ++gg) {
832  double a = v * t_w;
833  int rr = 0;
834  for (; rr != row_nb_dofs / 3; ++rr) {
835  auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
836  auto t_m = get_ftensor1(K, 3 * rr, 0);
837  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
838  double div_col_base = t_col_diff_base_fun(i, i);
839  t_m(i) += a * t_row_base_fun * div_col_base;
840  ++t_m;
841  ++t_col_diff_base_fun;
842  }
843  ++t_row_base_fun;
844  }
845  for (; rr != row_nb_base_functions; ++rr)
846  ++t_row_base_fun;
847  ++t_w;
848  }
850 }
851 
852 MoFEMErrorCode OpSpatialEquilibrium_dw_dw::integrate(EntData &row_data,
853  EntData &col_data) {
855 
856  if (alphaW < std::numeric_limits<double>::epsilon() &&
857  alphaRho < std::numeric_limits<double>::epsilon())
859 
860  const int nb_integration_pts = row_data.getN().size1();
861  const int row_nb_dofs = row_data.getIndices().size();
862  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
864  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
865 
866  );
867  };
869 
870  auto v = getVolume();
871  auto t_w = getFTensor0IntegrationWeight();
872 
873  int row_nb_base_functions = row_data.getN().size2();
874  auto t_row_base_fun = row_data.getFTensor0N();
875 
876  double ts_scale = alphaW * getTSa();
877  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
878  ts_scale += alphaRho * getTSaa();
879 
880  for (int gg = 0; gg != nb_integration_pts; ++gg) {
881  double a = v * t_w * ts_scale;
882 
883  int rr = 0;
884  for (; rr != row_nb_dofs / 3; ++rr) {
885 
886  auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
887  auto t_m = get_ftensor1(K, 3 * rr, 0);
888  for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
889  const double b = a * t_row_base_fun * t_col_base_fun;
890  t_m(i) -= b;
891  ++t_m;
892  ++t_col_base_fun;
893  }
894 
895  ++t_row_base_fun;
896  }
897 
898  for (; rr != row_nb_base_functions; ++rr)
899  ++t_row_base_fun;
900 
901  ++t_w;
902  }
903 
905 }
906 
907 MoFEMErrorCode OpSpatialPhysical_du_dP::integrate(EntData &row_data,
908  EntData &col_data) {
910  int nb_integration_pts = row_data.getN().size1();
911  int row_nb_dofs = row_data.getIndices().size();
912  int col_nb_dofs = col_data.getIndices().size();
913  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
915 
916  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
917 
918  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
919 
920  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
921 
922  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
923 
924  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
925 
926  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
927  };
928 
935 
936  auto v = getVolume();
937  auto t_w = getFTensor0IntegrationWeight();
938  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
939 
940  auto t_diff_u =
941  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
942 
943  int row_nb_base_functions = row_data.getN().size2();
944  auto t_row_base_fun = row_data.getFTensor0N();
945  for (int gg = 0; gg != nb_integration_pts; ++gg) {
946  double a = v * t_w;
947 
948  int rr = 0;
949  for (; rr != row_nb_dofs / 6; ++rr) {
950 
951  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
952  auto t_m = get_ftensor3(K, 6 * rr, 0);
953 
954  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
955 
957  t_RTP_dP(i, j, k) = t_R(k, i) * t_col_base_fun(j) * t_row_base_fun;
958 
959  for (int kk = 0; kk != 3; ++kk)
960  for (int ii = 0; ii != 3; ++ii)
961  for (int jj = ii; jj != 3; ++jj) {
962  for (int mm = 0; mm != 3; ++mm)
963  for (int nn = 0; nn != 3; ++nn)
964  t_m(ii, jj, kk) +=
965  a * t_diff_u(mm, nn, ii, jj) * t_RTP_dP(mm, nn, kk);
966  }
967 
968  ++t_col_base_fun;
969  ++t_m;
970  }
971 
972  ++t_row_base_fun;
973  }
974  for (; rr != row_nb_base_functions; ++rr)
975  ++t_row_base_fun;
976  ++t_w;
977  ++t_R;
978  ++t_diff_u;
979  }
981 }
982 
983 MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
984  EntData &col_data) {
986  int nb_integration_pts = row_data.getN().size1();
987  int row_nb_dofs = row_data.getIndices().size();
988  int col_nb_dofs = col_data.getIndices().size();
989  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
991  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
992  &m(r + 5, c));
993  };
994 
1000 
1001  auto v = getVolume();
1002  auto t_w = getFTensor0IntegrationWeight();
1003  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1004  auto t_diff_u =
1005  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
1006 
1007  int row_nb_base_functions = row_data.getN().size2();
1008  auto t_row_base_fun = row_data.getFTensor0N();
1009  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1010  double a = v * t_w;
1011 
1012  int rr = 0;
1013  for (; rr != row_nb_dofs / 6; ++rr) {
1014 
1015  auto t_m = get_ftensor2(K, 6 * rr, 0);
1016 
1017  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1018  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1019 
1020  FTensor::Tensor2<double, 3, 3> t_RTP_dBubble;
1021  t_RTP_dBubble(i, j) =
1022  a * t_row_base_fun * t_R(k, i) * t_col_base_fun(k, j);
1023  t_m(i, j) += t_diff_u(m, n, i, j) * t_RTP_dBubble(m, n);
1024 
1025  ++t_m;
1026  ++t_col_base_fun;
1027  }
1028 
1029  ++t_row_base_fun;
1030  }
1031  for (; rr != row_nb_base_functions; ++rr)
1032  ++t_row_base_fun;
1033  ++t_w;
1034  ++t_R;
1035  ++t_diff_u;
1036  }
1038 }
1039 
1040 MoFEMErrorCode OpSpatialPhysical_du_du::integrate(EntData &row_data,
1041  EntData &col_data) {
1043 
1044  int nb_integration_pts = row_data.getN().size1();
1045  int row_nb_dofs = row_data.getIndices().size();
1046  int col_nb_dofs = col_data.getIndices().size();
1047  auto get_ftensor4 = [](MatrixDouble &m, const int r, const int c) {
1049 
1050  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1051  &m(r + 0, c + 4), &m(r + 0, c + 5),
1052 
1053  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1054  &m(r + 1, c + 4), &m(r + 1, c + 5),
1055 
1056  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1057  &m(r + 2, c + 4), &m(r + 2, c + 5),
1058 
1059  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
1060  &m(r + 3, c + 4), &m(r + 3, c + 5),
1061 
1062  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
1063  &m(r + 4, c + 4), &m(r + 4, c + 5),
1064 
1065  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
1066  &m(r + 5, c + 4), &m(r + 5, c + 5)
1067 
1068  );
1069  };
1076 
1080 
1081  auto v = getVolume();
1082  auto t_w = getFTensor0IntegrationWeight();
1083  auto t_P_dh0 = getFTensor3FromMat(dataAtPts->P_dh0);
1084  auto t_P_dh1 = getFTensor3FromMat(dataAtPts->P_dh1);
1085  auto t_P_dh2 = getFTensor3FromMat(dataAtPts->P_dh2);
1086  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1087  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1088  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
1089  auto t_dot_log_u =
1090  getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachDotTensorAtPts);
1091  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1092  auto t_diff_u =
1093  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
1094  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1095  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1096  auto &nbUniq = dataAtPts->nbUniq;
1097 
1098  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
1100  t_one4(i, j, k, l) = t_kd(j, l) * t_kd(i, k);
1101 
1102  FTensor::Tensor4<double, 3, 3, 3, 3> t_one4_symmetric;
1103  t_one4_symmetric(i, j, k, l) = 0;
1104  for (auto ii : {0, 1, 2})
1105  for (auto jj : {0, 1, 2})
1106  for (auto kk : {0, 1, 2})
1107  for (auto ll : {0, 1, 2}) {
1108 
1109  if (ll != kk)
1110  t_one4_symmetric(ii, jj, kk, ll) =
1111  t_one4(ii, jj, kk, ll) + t_one4(ii, jj, ll, kk);
1112  else
1113  t_one4_symmetric(ii, jj, kk, ll) = t_one4(ii, jj, kk, ll);
1114  }
1115 
1116  int row_nb_base_functions = row_data.getN().size2();
1117  auto t_row_base_fun = row_data.getFTensor0N();
1118 
1120  t_one(i, j) = 0;
1121  for (auto ii : {0, 1, 2})
1122  t_one(ii, ii) = 1;
1123 
1124  const double ts_a = getTSa();
1125  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1126  double a = v * t_w;
1127 
1129  t_P_dh(i, j, N0, k) = t_P_dh0(i, j, k);
1130  t_P_dh(i, j, N1, k) = t_P_dh1(i, j, k);
1131  t_P_dh(i, j, N2, k) = t_P_dh2(i, j, k);
1132 
1134  t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
1136  t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
1137 
1139  t_dot_u(i, j) = t_u(i, k) * t_dot_log_u(k, j);
1140 
1142  t_stress_diff(i, j, k, l) = 0;
1143  for (int ii = 0; ii != 3; ++ii)
1144  for (int jj = 0; jj != 3; ++jj)
1145  for (int kk = 0; kk != 3; ++kk)
1146  for (int ll = 0; ll != 3; ++ll)
1147  for (int mm = 0; mm != 3; ++mm)
1148  for (int nn = 0; nn != 3; ++nn)
1149  t_stress_diff(ii, jj, mm, nn) -=
1150  t_P_dh(ii, jj, kk, ll) * t_diff_u(kk, ll, mm, nn);
1151 
1152  FTensor::Tensor2_symmetric<double, 3> t_viscous_stress;
1153  t_viscous_stress(i, j) = alphaU * t_dot_log_u(i, j);
1154  FTensor::Tensor4<double, 3, 3, 3, 3> t_viscous_stress_diff;
1155  t_viscous_stress_diff(i, j, k, l) = t_one4_symmetric(i, j, k, l);
1156  t_viscous_stress_diff(i, j, k, l) *= (ts_a * alphaU);
1157 
1159  t_deltaP2(i, j) = t_deltaP(i, j) - t_viscous_stress(i, j);
1160 
1161  auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
1162  t_eigen_vals, t_eigen_vecs, f, d_f, dd_f, t_deltaP2, nbUniq[gg]);
1163  for (int ii = 0; ii != 3; ++ii) {
1164  for (int jj = 0; jj < ii; ++jj) {
1165  for (int kk = 0; kk != 3; ++kk) {
1166  for (int ll = 0; ll != kk; ++ll) {
1167  t_diff2_uP2(ii, jj, kk, ll) *= 2;
1168  }
1169  }
1170  }
1171  }
1172  for (int ii = 0; ii != 3; ++ii) {
1173  for (int jj = ii; jj != 3; ++jj) {
1174  for (int kk = 0; kk != 3; ++kk) {
1175  for (int ll = 0; ll < kk; ++ll) {
1176  t_diff2_uP2(ii, jj, kk, ll) *= 2;
1177  }
1178  }
1179  }
1180  }
1181 
1182  int rr = 0;
1183  for (; rr != row_nb_dofs / 6; ++rr) {
1184 
1185  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1186  auto t_m = get_ftensor4(K, 6 * rr, 0);
1187 
1188  for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1189 
1190  const double b = a * t_row_base_fun * t_col_base_fun;
1191 
1192  for (int ii = 0; ii != 3; ++ii)
1193  for (int jj = ii; jj != 3; ++jj)
1194  for (int kk = 0; kk != 3; ++kk)
1195  for (int ll = kk; ll != 3; ++ll)
1196  for (int mm = 0; mm != 3; ++mm)
1197  for (int nn = 0; nn != 3; ++nn)
1198  t_m(ii, jj, kk, ll) +=
1199  b * t_diff_u(mm, nn, ii, jj) *
1200  (t_stress_diff(mm, nn, kk, ll) -
1201  t_viscous_stress_diff(mm, nn, kk, ll));
1202 
1203  t_m(i, j, k, l) += b * t_diff2_uP2(i, j, k, l);
1204 
1205  ++t_m;
1206  ++t_col_base_fun;
1207  }
1208 
1209  ++t_row_base_fun;
1210  }
1211  for (; rr != row_nb_base_functions; ++rr)
1212  ++t_row_base_fun;
1213 
1214  ++t_w;
1215  ++t_P_dh0;
1216  ++t_P_dh1;
1217  ++t_P_dh2;
1218  ++t_R;
1219  ++t_approx_P;
1220  ++t_P;
1221  ++t_dot_log_u;
1222  ++t_u;
1223  ++t_diff_u;
1224  ++t_eigen_vals;
1225  ++t_eigen_vecs;
1226  }
1228 }
1229 
1230 MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
1231  EntData &col_data) {
1233  int row_nb_dofs = row_data.getIndices().size();
1234  int col_nb_dofs = col_data.getIndices().size();
1235  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1237 
1238  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1239 
1240  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1241 
1242  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1243 
1244  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1245 
1246  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1247 
1248  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1249 
1250  );
1251  };
1256 
1257  auto v = getVolume();
1258  auto t_w = getFTensor0IntegrationWeight();
1259  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1260  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1261  auto t_diff_u =
1262  getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
1263 
1264  int row_nb_base_functions = row_data.getN().size2();
1265  auto t_row_base_fun = row_data.getFTensor0N();
1266 
1267  int nb_integration_pts = row_data.getN().size1();
1268 
1269  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1270  double a = v * t_w;
1271 
1272  FTensor::Tensor3<double, 3, 3, 3> t_RT_domega_P;
1273  t_RT_domega_P(i, j, m) = t_diff_R(k, i, m) * t_P(k, j);
1275  t(i, j, k) = 0;
1276  for (int ii = 0; ii != 3; ++ii)
1277  for (int jj = ii; jj != 3; ++jj)
1278  for (int mm = 0; mm != 3; ++mm)
1279  for (int nn = 0; nn != 3; ++nn)
1280  for (int kk = 0; kk != 3; ++kk)
1281  t(ii, jj, kk) +=
1282  t_diff_u(mm, nn, ii, jj) * t_RT_domega_P(mm, nn, kk);
1283 
1284  int rr = 0;
1285  for (; rr != row_nb_dofs / 6; ++rr) {
1286  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1287  auto t_m = get_ftensor3(K, 6 * rr, 0);
1288  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1289 
1290  double v = a * t_row_base_fun * t_col_base_fun;
1291  t_m(i, j, k) += v * t(i, j, k);
1292 
1293  ++t_m;
1294  ++t_col_base_fun;
1295  }
1296 
1297  ++t_row_base_fun;
1298  }
1299 
1300  for (; rr != row_nb_base_functions; ++rr)
1301  ++t_row_base_fun;
1302 
1303  ++t_w;
1304  ++t_P;
1305  ++t_diff_R;
1306  ++t_diff_u;
1307  }
1308 
1310 }
1311 
1312 MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
1313  EntData &col_data) {
1315  int nb_integration_pts = row_data.getN().size1();
1316  int row_nb_dofs = row_data.getIndices().size();
1317  int col_nb_dofs = col_data.getIndices().size();
1318  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1320 
1321  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1322 
1323  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1324 
1325  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1326 
1327  );
1328  };
1335 
1336  auto v = getVolume();
1337  auto t_w = getFTensor0IntegrationWeight();
1338  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1339 
1340  int row_nb_base_functions = row_data.getN().size2();
1341  auto t_row_base_fun = row_data.getFTensor0N();
1342 
1343  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1344  double a = v * t_w;
1345 
1346  int rr = 0;
1347  for (; rr != row_nb_dofs / 3; ++rr) {
1348  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1349  auto t_m = get_ftensor2(K, 3 * rr, 0);
1351  t_levi_R(i, j, k) =
1352  (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1353  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1354  FTensor::Tensor1<double, 3> t_pushed_base;
1355  t_m(k, i) += t_levi_R(i, j, k) * t_col_base_fun(j);
1356  ++t_m;
1357  ++t_col_base_fun;
1358  }
1359 
1360  ++t_row_base_fun;
1361  }
1362 
1363  for (; rr != row_nb_base_functions; ++rr)
1364  ++t_row_base_fun;
1365  ++t_w;
1366  ++t_R;
1367  }
1369 }
1370 
1371 MoFEMErrorCode OpSpatialRotation_domega_dBubble::integrate(EntData &row_data,
1372  EntData &col_data) {
1374  int nb_integration_pts = row_data.getN().size1();
1375  int row_nb_dofs = row_data.getIndices().size();
1376  int col_nb_dofs = col_data.getIndices().size();
1377  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1379  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
1380  };
1386 
1387  auto v = getVolume();
1388  auto t_w = getFTensor0IntegrationWeight();
1389  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1390 
1391  int row_nb_base_functions = row_data.getN().size2();
1392  auto t_row_base_fun = row_data.getFTensor0N();
1393 
1394  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1395  double a = v * t_w;
1396 
1397  int rr = 0;
1398  for (; rr != row_nb_dofs / 3; ++rr) {
1399 
1400  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1401  auto t_m = get_ftensor1(K, 3 * rr, 0);
1402 
1404  t_levi_r(i, j, k) =
1405  (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1406 
1407  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1408  FTensor::Tensor2<double, 3, 3> t_pushed_base;
1409  t_m(k) += t_levi_r(i, j, k) * t_col_base_fun(i, j);
1410  ++t_m;
1411  ++t_col_base_fun;
1412  }
1413 
1414  ++t_row_base_fun;
1415  }
1416 
1417  for (; rr != row_nb_base_functions; ++rr)
1418  ++t_row_base_fun;
1419  ++t_w;
1420  ++t_R;
1421  }
1423 }
1424 
1425 MoFEMErrorCode OpSpatialRotation_domega_domega::integrate(EntData &row_data,
1426  EntData &col_data) {
1428  int nb_integration_pts = row_data.getN().size1();
1429  int row_nb_dofs = row_data.getIndices().size();
1430  int col_nb_dofs = col_data.getIndices().size();
1431  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1433 
1434  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1435 
1436  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1437 
1438  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1439 
1440  );
1441  };
1448 
1449  auto v = getVolume();
1450  auto t_w = getFTensor0IntegrationWeight();
1451  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1452  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1453  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1454 
1455  int row_nb_base_functions = row_data.getN().size2();
1456  auto t_row_base_fun = row_data.getFTensor0N();
1457 
1458  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1459  double a = v * t_w;
1460 
1462  t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
1463 
1464  FTensor::Tensor2<double, 3, 3> t_leviPRT_domega;
1465  t_leviPRT_domega(n, l) =
1466  levi_civita(i, j, n) * (t_approx_P(i, k) * t_diff_R(j, k, l));
1467 
1468  int rr = 0;
1469  for (; rr != row_nb_dofs / 3; ++rr) {
1470 
1471  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1472  auto t_m = get_ftensor2(K, 3 * rr, 0);
1473  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1474 
1475  const double b = a * t_row_base_fun * t_col_base_fun;
1476 
1477  t_m(k, m) += b * t_leviPRT_domega(k, m);
1478 
1479  ++t_m;
1480  ++t_col_base_fun;
1481  }
1482 
1483  ++t_row_base_fun;
1484  }
1485 
1486  for (; rr != row_nb_base_functions; ++rr)
1487  ++t_row_base_fun;
1488  ++t_w;
1489  ++t_approx_P;
1490  ++t_R;
1491  ++t_diff_R;
1492  }
1494 }
1495 
1496 MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
1497  EntData &col_data) {
1499 
1500  int nb_integration_pts = row_data.getN().size1();
1501  int row_nb_dofs = row_data.getIndices().size();
1502  int col_nb_dofs = col_data.getIndices().size();
1503 
1504  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1506 
1507  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1508 
1509  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1510 
1511  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1512 
1513  );
1514  };
1515 
1522 
1523  auto v = getVolume();
1524  auto t_w = getFTensor0IntegrationWeight();
1525  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1526  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1527  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1528  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
1529  int row_nb_base_functions = row_data.getN().size2() / 3;
1530  auto t_row_base_fun = row_data.getFTensor1N<3>();
1531 
1532  const double ts_a = getTSa();
1533 
1534  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1535  double a = v * t_w;
1536 
1537  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
1538 
1539  int rr = 0;
1540  for (; rr != row_nb_dofs / 3; ++rr) {
1541 
1543  t_PRT(i, k) = t_row_base_fun(j) * (t_diff_R(i, l, k) * t_u(l, j));
1544 
1545  FTensor::Tensor2<double, 3, 3> t_levi1, t_levi2;
1546  t_levi1(i, k) = (levi_civita(i, m, n) * t_omega_dot(n)) *
1547  (t_row_base_fun(j) * t_diff_R(m, j, k));
1548  t_levi2(i, k) = levi_civita(i, m, k) * (t_row_base_fun(j) * t_R(m, j));
1549 
1550  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1551  auto t_m = get_ftensor2(K, 3 * rr, 0);
1552  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1553  t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
1554  t_m(i, j) += (a * t_col_base_fun) * t_levi1(i, j);
1555  t_m(i, j) += ((a * ts_a) * t_col_base_fun) * t_levi2(i, j);
1556  ++t_m;
1557  ++t_col_base_fun;
1558  }
1559 
1560  ++t_row_base_fun;
1561  }
1562 
1563  for (; rr != row_nb_base_functions; ++rr)
1564  ++t_row_base_fun;
1565  ++t_w;
1566  ++t_R;
1567  ++t_diff_R;
1568  ++t_u;
1569  ++t_omega_dot;
1570  }
1572 }
1573 
1575 OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
1576  EntData &col_data) {
1578 
1579  int nb_integration_pts = row_data.getN().size1();
1580  int row_nb_dofs = row_data.getIndices().size();
1581  int col_nb_dofs = col_data.getIndices().size();
1582 
1583  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1585  &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
1586  };
1587 
1594 
1595  auto v = getVolume();
1596  auto t_w = getFTensor0IntegrationWeight();
1597  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1598  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
1599  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
1600  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
1601  int row_nb_base_functions = row_data.getN().size2() / 9;
1602  auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
1603 
1604  const double ts_a = getTSa();
1605 
1606  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1607  double a = v * t_w;
1608 
1609  int rr = 0;
1610  for (; rr != row_nb_dofs; ++rr) {
1611 
1613  t_PRT(k) = t_row_base_fun(i, j) * (t_diff_R(i, l, k) * t_u(l, j));
1614 
1615  FTensor::Tensor2<double, 3, 3> t_levi_omega;
1616  t_levi_omega(i, m) = levi_civita(i, m, n) * t_omega_dot(n);
1617  FTensor::Tensor3<double, 3, 3, 3> t_levi_omegaDiffR;
1618  t_levi_omegaDiffR(i, j, k) = t_levi_omega(i, m) * t_diff_R(m, j, k);
1619  FTensor::Tensor3<double, 3, 3, 3> t_levi_omegaR;
1620  t_levi_omegaR(i, j, k) = levi_civita(i, m, k) * t_R(m, j);
1621 
1622  FTensor::Tensor1<double, 3> t_levi1, t_levi2;
1623  t_levi1(k) = t_row_base_fun(i, j) * t_levi_omegaDiffR(i, j, k);
1624  t_levi2(k) = t_row_base_fun(i, j) * t_levi_omegaR(i, j, k);
1625 
1626  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1627  auto t_m = get_ftensor2(K, rr, 0);
1628  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1629  t_m(j) += (a * t_col_base_fun) * t_PRT(j);
1630  t_m(j) += (a * t_col_base_fun) * t_levi1(j);
1631  t_m(j) += ((a * ts_a) * t_col_base_fun) * t_levi2(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 
1642  ++t_w;
1643  ++t_R;
1644  ++t_diff_R;
1645  ++t_u;
1646  ++t_omega_dot;
1647  }
1649 }
1650 
1651 MoFEMErrorCode OpSpatialSchurBegin::assemble(int row_side, EntityType row_type,
1652  EntData &data) {
1654  if (row_type != MBTET)
1656  auto &bmc = dataAtPts->blockMatContainor;
1657  for (auto &bit : bmc)
1658  bit.unSetAtElement();
1659  // bmc.clear();
1661 }
1662 
1663 MoFEMErrorCode OpSpatialPreconditionMass::integrate(EntData &row_data) {
1665  const int nb_integration_pts = row_data.getN().size1();
1666  const int row_nb_dofs = row_data.getIndices().size();
1667  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1669  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
1670 
1671  );
1672  };
1674 
1675  auto v = getVolume();
1676  auto t_w = getFTensor0IntegrationWeight();
1677  auto &m = *mPtr;
1678 
1679  m.resize(row_nb_dofs, row_nb_dofs, false);
1680  m.clear();
1681 
1682  int row_nb_base_functions = row_data.getN().size2();
1683  auto t_row_base_fun = row_data.getFTensor0N();
1684 
1685  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1686  double a = v * t_w;
1687 
1688  int rr = 0;
1689  for (; rr != row_nb_dofs / 3; ++rr) {
1690 
1691  auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
1692  auto t_m = get_ftensor1(m, 3 * rr, 0);
1693  for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1694  const double b = a * t_row_base_fun * t_col_base_fun;
1695  t_m(i) += b;
1696  ++t_m;
1697  ++t_col_base_fun;
1698  }
1699 
1700  ++t_row_base_fun;
1701  }
1702 
1703  for (; rr != row_nb_base_functions; ++rr)
1704  ++t_row_base_fun;
1705 
1706  ++t_w;
1707  }
1708 
1710 }
1711 
1712 MoFEMErrorCode OpSpatialSchurEnd::assemble(int row_side, EntityType row_type,
1713  EntData &data) {
1715  if (row_type != MBTET)
1717 
1718  if (auto ep_fe_vol_ptr =
1719  dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
1720  getFEMethod())) {
1721 
1722  SmartPetscObj<AO> aoSuu;
1723  SmartPetscObj<AO> aoSBB;
1724  SmartPetscObj<AO> aoSOO;
1725  SmartPetscObj<AO> aoSww;
1726 
1727  SmartPetscObj<Mat> Suu;
1728  SmartPetscObj<Mat> SBB;
1729  SmartPetscObj<Mat> SOO;
1730  SmartPetscObj<Mat> Sww;
1731 
1732  auto set_data = [&](auto fe) {
1733  aoSuu = fe->aoSuu;
1734  aoSBB = fe->aoSBubble;
1735  aoSOO = fe->aoSOmega;
1736  aoSww = fe->aoSw;
1737 
1738  Suu = fe->Suu;
1739  SBB = fe->SBubble;
1740  SOO = fe->SOmega;
1741  Sww = fe->Sw;
1742  };
1743  set_data(ep_fe_vol_ptr);
1744 
1745  if (Suu) {
1746 
1747  auto find_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1748  auto &row_name, auto &col_name, auto row_side,
1749  auto col_side, auto row_type, auto col_type) {
1750  return add_bmc.get<0>().find(boost::make_tuple(
1751  row_name, col_name, row_type, col_type, row_side, col_side));
1752  };
1753 
1754  auto set_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1755  auto &row_name, auto &col_name, auto row_side,
1756  auto col_side, auto row_type, auto col_type,
1757  const auto &m, const auto &row_ind,
1758  const auto &col_ind) {
1759  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
1760  row_type, col_type);
1761  if (it != add_bmc.get<0>().end()) {
1762  it->setMat(m);
1763  it->setInd(row_ind, col_ind);
1764  it->setSetAtElement();
1765  return it;
1766  } else {
1767  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
1768  row_name, col_name, row_type, col_type, row_side, col_side, m,
1769  row_ind, col_ind));
1770  return p_it.first;
1771  }
1772  };
1773 
1774  auto add_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1775  auto &row_name, auto &col_name, auto row_side,
1776  auto col_side, auto row_type, auto col_type,
1777  const auto &m, const auto &row_ind,
1778  const auto &col_ind) {
1779  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
1780  row_type, col_type);
1781  if (it != add_bmc.get<0>().end()) {
1782  it->addMat(m);
1783  it->setInd(row_ind, col_ind);
1784  it->setSetAtElement();
1785  return it;
1786  } else {
1787  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
1788  row_name, col_name, row_type, col_type, row_side, col_side, m,
1789  row_ind, col_ind));
1790  return p_it.first;
1791  }
1792  };
1793 
1794  auto assemble_block = [&](auto &bit, Mat S) {
1796  const VectorInt &rind = bit.rowInd;
1797  const VectorInt &cind = bit.colInd;
1798  const MatrixDouble &m = bit.M;
1799  CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
1800  &*cind.begin(), &*m.data().begin(), ADD_VALUES);
1801 
1803  };
1804 
1805  auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
1807  const int nb = m.size1();
1808 
1809  inv.resize(nb, nb, false);
1810  inv.clear();
1811  for (int cc = 0; cc != nb; ++cc)
1812  inv(cc, cc) = -1;
1813 
1814  iPIV.resize(nb, false);
1815  lapackWork.resize(nb * nb, false);
1816  const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
1817  &*iPIV.begin(), &*inv.data().begin(), nb,
1818  &*lapackWork.begin(), nb * nb);
1819  // if (info != 0)
1820  // SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1821  // "Can not invert matrix info = %d", info);
1822 
1824  };
1825 
1826  auto invert_symm_schur = [&](DataAtIntegrationPts::BlockMatContainor &bmc,
1827  std::string field, auto &inv) {
1829 
1830  auto bit =
1831  bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
1832  if (bit != bmc.get<1>().end()) {
1833 
1834  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
1835  CHKERR invert_symm_mat(m, inv);
1836 
1837  } else
1838  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1839  "%s matrix not found", field.c_str());
1840 
1842  };
1843 
1844  auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
1846 
1847  const int nb = m.size1();
1848 
1849  MatrixDouble trans_m = trans(m);
1850  MatrixDouble trans_inv;
1851  trans_inv.resize(nb, nb, false);
1852  trans_inv.clear();
1853  for (int c = 0; c != nb; ++c)
1854  trans_inv(c, c) = -1;
1855 
1856  iPIV.resize(nb, false);
1857  const auto info =
1858  lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb, &*iPIV.begin(),
1859  &*trans_inv.data().begin(), nb);
1860  if (info != 0)
1861  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1862  "Can not invert matrix info = %d", info);
1863 
1864  inv.resize(nb, nb, false);
1865  noalias(inv) = trans(trans_inv);
1866 
1868  };
1869 
1870  auto invert_nonsymm_schur =
1871  [&](DataAtIntegrationPts::BlockMatContainor &bmc, std::string field,
1872  auto &inv, const bool debug = false) {
1874 
1875  auto bit = bmc.get<1>().find(
1876  boost::make_tuple(field, field, MBTET, MBTET));
1877  if (bit != bmc.get<1>().end()) {
1878 
1879  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
1880  CHKERR invert_nonsymm_mat(m, inv);
1881 
1882  if (debug) {
1883  std::cerr << prod(m, inv) << endl;
1884  std::cerr << endl;
1885  }
1886 
1887  } else
1888  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1889  "%s matrix not found", field.c_str());
1890 
1892  };
1893 
1894  auto create_block_schur =
1897  std::string field, AO ao, MatrixDouble &inv) {
1899 
1900  for (auto &bit : add_bmc) {
1901  bit.unSetAtElement();
1902  bit.clearMat();
1903  }
1904 
1905  for (auto &bit : bmc) {
1906  if (bit.setAtElement && bit.rowField != field &&
1907  bit.colField != field) {
1908  VectorInt rind = bit.rowInd;
1909  VectorInt cind = bit.colInd;
1910  const MatrixDouble &m = bit.M;
1911  if (ao) {
1912  CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
1913  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1914  }
1915  auto it = set_block(add_bmc, bit.rowField, bit.colField,
1916  bit.rowSide, bit.colSide, bit.rowType,
1917  bit.colType, m, rind, cind);
1918  }
1919  }
1920 
1921  for (auto &bit_col : bmc) {
1922  if (bit_col.setAtElement && bit_col.rowField == field &&
1923  bit_col.colField != field) {
1924  const MatrixDouble &cm = bit_col.M;
1925  VectorInt cind = bit_col.colInd;
1926  invMat.resize(inv.size1(), cm.size2(), false);
1927  noalias(invMat) = prod(inv, cm);
1928  if (ao)
1929  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
1930  for (auto &bit_row : bmc) {
1931  if (bit_row.setAtElement && bit_row.rowField != field &&
1932  bit_row.colField == field) {
1933  const MatrixDouble &rm = bit_row.M;
1934  VectorInt rind = bit_row.rowInd;
1935  K.resize(rind.size(), cind.size(), false);
1936  noalias(K) = prod(rm, invMat);
1937  if (ao)
1938  CHKERR AOApplicationToPetsc(ao, rind.size(),
1939  &*rind.begin());
1940  auto it = add_block(add_bmc, bit_row.rowField,
1941  bit_col.colField, bit_row.rowSide,
1942  bit_col.colSide, bit_row.rowType,
1943  bit_col.colType, K, rind, cind);
1944  }
1945  }
1946  }
1947  }
1948 
1950  };
1951 
1952  auto assemble_schur =
1953  [&](DataAtIntegrationPts::BlockMatContainor &add_bmc, Mat S,
1954  bool debug = false) {
1956  for (auto &bit : add_bmc) {
1957  if (bit.setAtElement)
1958  CHKERR assemble_block(bit, S);
1959  }
1960  if (debug) {
1961  for (auto &bit : add_bmc) {
1962  if (bit.setAtElement) {
1963  std::cerr << "assemble: " << bit.rowField << " "
1964  << bit.colField << endl;
1965  std::cerr << bit.M << endl;
1966  }
1967  }
1968  std::cerr << std::endl;
1969  }
1971  };
1972 
1973  auto precondition_schur =
1976  const std::string field, const MatrixDouble &diag_mat,
1977  const double eps) {
1979 
1980  for (auto &bit : add_bmc) {
1981  bit.unSetAtElement();
1982  bit.clearMat();
1983  }
1984 
1985  for (auto &bit : bmc) {
1986  if (bit.setAtElement) {
1987  if (bit.rowField != field || bit.colField != field)
1988  auto it =
1989  set_block(add_bmc, bit.rowField, bit.colField,
1990  bit.rowSide, bit.colSide, bit.rowType,
1991  bit.colType, bit.M, bit.rowInd, bit.colInd);
1992  }
1993  }
1994 
1995  auto bit = bmc.get<1>().find(
1996  boost::make_tuple(field, field, MBTET, MBTET));
1997  if (bit->setAtElement && bit != bmc.get<1>().end()) {
1998  auto it =
1999  set_block(add_bmc, bit->rowField, bit->colField, bit->rowSide,
2000  bit->colSide, bit->rowType, bit->colType, bit->M,
2001  bit->rowInd, bit->colInd);
2002  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2003  m += eps * diag_mat;
2004  } else {
2005  auto row_it = bmc.get<3>().lower_bound(field);
2006  for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
2007  if (row_it->setAtElement) {
2008  auto it = set_block(add_bmc, field, field, 0, 0, MBTET, MBTET,
2009  diag_mat, row_it->rowInd, row_it->rowInd);
2010  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2011  m *= eps;
2012  break;
2013  }
2014  }
2015  if (row_it == bmc.get<3>().end())
2016  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2017  "row field not found %s", field.c_str());
2018  }
2019 
2021  };
2022 
2023  CHKERR invert_symm_schur(dataAtPts->blockMatContainor, "u",
2024  invBlockMat["uu"]);
2025  CHKERR create_block_schur(dataAtPts->blockMatContainor, blockMat["uu"],
2026  "u", aoSuu, invBlockMat["uu"]);
2027  CHKERR assemble_schur(blockMat["uu"], Suu);
2028 
2029  if (SBB) {
2030  CHKERR invert_symm_schur(blockMat["uu"], "BUBBLE", invBlockMat["BB"]);
2031  CHKERR create_block_schur(blockMat["uu"], blockMat["BB"], "BUBBLE",
2032  aoSBB, invBlockMat["BB"]);
2033  CHKERR precondition_schur(blockMat["BB"], blockMat["precBBOO"], "omega",
2034  *dataAtPts->ooMatPtr, eps);
2035  CHKERR assemble_schur(blockMat["precBBOO"], SBB);
2036 
2037  if (SOO) {
2038  CHKERR invert_nonsymm_schur(blockMat["precBBOO"], "omega",
2039  invBlockMat["OO"]);
2040  CHKERR create_block_schur(blockMat["precBBOO"], blockMat["OO"],
2041  "omega", aoSOO, invBlockMat["OO"]);
2042  if (dataAtPts->wwMatPtr) {
2043  CHKERR precondition_schur(blockMat["OO"], blockMat["precOOww"], "w",
2044  *dataAtPts->wwMatPtr, -eps);
2045  } else {
2046  blockMat["precOOww"] = blockMat["OO"];
2047  }
2048  CHKERR assemble_schur(blockMat["precOOww"], SOO);
2049 
2050  if (Sww) {
2051  CHKERR invert_symm_schur(blockMat["precOOww"], "w",
2052  invBlockMat["ww"]);
2053  CHKERR create_block_schur(blockMat["precOOww"], blockMat["ww"], "w",
2054  aoSww, invBlockMat["ww"]);
2055  CHKERR assemble_schur(blockMat["ww"], Sww);
2056  }
2057  }
2058  }
2059  }
2060  }
2061 
2063 }
2064 
2065 MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
2066  EntData &data) {
2068  if (type != MBVERTEX)
2070 
2071  auto create_tag = [this](const std::string tag_name, const int size) {
2072  double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2073  Tag th;
2074  CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
2075  th, MB_TAG_CREAT | MB_TAG_SPARSE,
2076  def_VAL);
2077  return th;
2078  };
2079 
2080  Tag th_w = create_tag("SpatialDisplacement", 3);
2081  Tag th_omega = create_tag("Omega", 3);
2082  Tag th_approxP = create_tag("Piola1Stress", 9);
2083  Tag th_sigma = create_tag("CauchyStress", 9);
2084  Tag th_log_u = create_tag("LogSpatialStreach", 9);
2085  Tag th_u = create_tag("SpatialStreach", 9);
2086  Tag th_h = create_tag("h", 9);
2087  Tag th_X = create_tag("X", 3);
2088  Tag th_detF = create_tag("detF", 1);
2089  Tag th_angular_momentum = create_tag("AngularMomentum", 3);
2090 
2091  Tag th_u_eig_vec = create_tag("SpatialStreachEigenVec", 9);
2092  Tag th_u_eig_vals = create_tag("SpatialStreachEigenVals", 3);
2093 
2094  int nb_gauss_pts = data.getN().size1();
2095  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
2096  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2097  "Nb. of integration points is not equal to number points on "
2098  "post-processing mesh");
2099  }
2100 
2101  auto t_w = getFTensor1FromMat<3>(dataAtPts->wAtPts);
2102  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
2103  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2104  auto t_log_u =
2105  getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachTensorAtPts);
2106  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
2107  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
2108  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2109  auto t_coords = getFTensor1CoordsAtGaussPts();
2110 
2115 
2116  // vectors
2117  VectorDouble3 v(3);
2118  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
2119  auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
2121  t_v(i) = t_d(i);
2122  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2123  &*v.data().begin());
2125  };
2126 
2127  MatrixDouble3by3 m(3, 3);
2129  &m(0, 0), &m(0, 1), &m(0, 2),
2130 
2131  &m(1, 0), &m(1, 1), &m(1, 2),
2132 
2133  &m(2, 0), &m(2, 1), &m(2, 2));
2134 
2135  auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
2137  t_m(i, j) = t_d(i, j);
2138  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2139  &*m.data().begin());
2141  };
2142 
2143  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2144 
2145  // vetors
2146  CHKERR save_vec_tag(th_w, t_w, gg);
2147  CHKERR save_vec_tag(th_X, t_coords, gg);
2148  CHKERR save_vec_tag(th_omega, t_omega, gg);
2149 
2150  // tensors
2151  CHKERR save_mat_tag(th_h, t_h, gg);
2152 
2153  FTensor::Tensor2<double, 3, 3> t_log_u_tmp;
2154  for (int ii = 0; ii != 3; ++ii)
2155  for (int jj = 0; jj != 3; ++jj)
2156  t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
2157 
2158  CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
2159 
2161  for (int ii = 0; ii != 3; ++ii)
2162  for (int jj = 0; jj != 3; ++jj)
2163  t_u_tmp(ii, jj) = t_u(ii, jj);
2164 
2165  CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2166  CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2167 
2168  const double jac = dEterminant(t_h);
2170  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
2171  CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2172  CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
2173 
2175  t_PhT(i, k) = t_approx_P(i, j) * t_R(k, j);
2176  FTensor::Tensor1<double, 3> t_leviPRT;
2177  t_leviPRT(k) = levi_civita(i, l, k) * t_PhT(i, l);
2178 
2179  CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
2180  &t_leviPRT(0));
2181 
2182  auto get_eiegn_vector_symmetric = [&](auto &t_u) {
2184 
2185  for (int ii = 0; ii != 3; ++ii)
2186  for (int jj = 0; jj != 3; ++jj)
2187  t_m(ii, jj) = t_u(ii, jj);
2188 
2189  VectorDouble3 eigen_values(3);
2190 
2191  // LAPACK - eigenvalues and vectors. Applied twice for initial creates
2192  // memory space
2193  int n = 3, lda = 3, info, lwork = -1;
2194  double wkopt;
2195  info = lapack_dsyev('V', 'U', n, &(m.data()[0]), lda,
2196  &(eigen_values.data()[0]), &wkopt, lwork);
2197  if (info != 0)
2198  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2199  "is something wrong with lapack_dsyev info = %d", info);
2200  lwork = (int)wkopt;
2201  double work[lwork];
2202  info = lapack_dsyev('V', 'U', n, &(m.data()[0]), lda,
2203  &(eigen_values.data()[0]), work, lwork);
2204  if (info != 0)
2205  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2206  "is something wrong with lapack_dsyev info = %d", info);
2207 
2208  CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
2209  &*m.data().begin());
2210  CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
2211  &*eigen_values.data().begin());
2212 
2214  };
2215 
2216  CHKERR get_eiegn_vector_symmetric(t_u);
2217 
2218  ++t_w;
2219  ++t_h;
2220  ++t_log_u;
2221  ++t_u;
2222  ++t_omega;
2223  ++t_R;
2224  ++t_approx_P;
2225  ++t_coords;
2226  }
2227 
2229 }
2230 
2232 OpCalculateStrainEnergy::doWork(int side, EntityType type,
2235  if (type == MBTET) {
2236  int nb_integration_pts = data.getN().size1();
2237  auto v = getVolume();
2238  auto t_w = getFTensor0IntegrationWeight();
2239  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2240  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2241 
2244 
2245  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2246  const double a = t_w * v;
2247  (*energy) += a * t_P(i, J) * t_h(i, J);
2248  ++t_w;
2249  ++t_P;
2250  ++t_h;
2251  }
2252  }
2254 }
2255 
2256 MoFEMErrorCode OpL2Transform::doWork(int side, EntityType type,
2259  if (type == MBTET) {
2260  auto v = getVolume();
2261  for (auto &base : data.getN(AINSWORTH_LEGENDRE_BASE).data())
2262  base /= v;
2263  }
2265 }
2266 
2267 } // namespace EshelbianPlasticity
static Index< 'm', 3 > m
static Index< 'n', 3 > n
static Index< 'J', 3 > J
static Number< 2 > N2
static Number< 1 > N1
static Index< 'K', 3 > K
static Number< 0 > N0
MoFEMErrorCode get_eigen_val_and_proj_lapack(const Tensor2_symmetric< T, 3 > &X, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec)
bool is_eq(const double &a, const double &b)
Eshelbian plasticity interface.
Kronecker Delta class.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:126
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
@ MOFEM_INVALID_DATA
Definition: definitions.h:128
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
static const bool debug
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
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
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:273
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::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
double tol
const double T
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.
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
Get the Diff Diff Mat object.
static FTensor::Tensor3< typename TensorTypeExtractor< T >::Type, 3, 3, 3 > get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
double f(const double v)
void sort_eigen_vals(T1 &eig, T2 &eigen_vec)
double dd_f(const double v)
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
double d_f(const double v)
JSON compatible output.
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
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:102
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:875
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
const double r
rate factor
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
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
const VectorDouble & getFieldData() const
get dofs values
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
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.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.