v0.9.0
EshelbianOperators.cpp
Go to the documentation of this file.
1 /**
2  * \file EshelbianOperators.cpp
3  * \example EshelbianOperators.cpp
4  *
5  * \brief Implementation of operators
6  */
7 
8 #include <MoFEM.hpp>
9 using namespace MoFEM;
10 
11 #include <BasicFiniteElements.hpp>
12 
13 #include <EshelbianPlasticity.hpp>
14 #include <boost/math/constants/constants.hpp>
15 
16 #include <lapack_wrap.h>
17 
18 namespace EshelbianPlasticity {
19 
20 MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
21  EntityType type,
22  EntData &data) {
24  if (type != MBTET)
26  int nb_integration_pts = data.getN().size1();
31 
32  dataAtPts->hAtPts->resize(9, nb_integration_pts, false);
33  dataAtPts->hDeltaAtPts->resize(9, nb_integration_pts, false);
34  dataAtPts->xDetAtPts->resize(nb_integration_pts, false);
35  dataAtPts->xInvGradAtPts->resize(9, nb_integration_pts, false);
36  dataAtPts->rotMatAtPts->resize(9, nb_integration_pts, false);
37  dataAtPts->diffRotMatAtPts->resize(27, nb_integration_pts, false);
38 
39  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
40  auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
41  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
42  auto t_x_det = getFTensor0FromVec(*(dataAtPts->xDetAtPts));
43  auto t_x_inv_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xInvGradAtPts));
44  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
45  auto t_omega = getFTensor1FromMat<3>(*(dataAtPts->rotAxisAtPts));
46  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
47  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
48 
49  for (int gg = 0; gg != nb_integration_pts; ++gg) {
50 
51  CHKERR determinantTensor3by3(t_x_grad, t_x_det);
52  CHKERR invertTensor3by3(t_x_grad, t_x_det, t_x_inv_grad);
53 
54  auto t0_diff = getDiffRotationFormVector(t_omega);
55  auto t0 = getRotationFormVector(t_omega);
56 
57  t_diff_R(i, j, k) = t0_diff(i, j, k);
58  t_R(i, j) = t0(i, j);
59  t_delta_h(i, j) = t_R(i, k) * t_u(k, j) + t_R(i, j);
60  t_h(i, j) = t_delta_h(i, k) * t_x_grad(k, j);
61 
62  ++t_h;
63  ++t_delta_h;
64  ++t_x_grad;
65  ++t_x_det;
66  ++t_x_inv_grad;
67  ++t_R;
68  ++t_diff_R;
69  ++t_u;
70  ++t_omega;
71  }
72 
74 }
75 
76 MoFEMErrorCode OpSpatialEquilibrium::integrate(EntData &data) {
78  int nb_dofs = data.getIndices().size();
79  int nb_integration_pts = data.getN().size1();
80  auto v = getVolume();
81  auto t_w = getFTensor0IntegrationWeight();
82  auto t_div_P = getFTensor1FromMat<3>(*(dataAtPts->divPAtPts));
83  int nb_base_functions = data.getN().size2();
84  auto t_row_base_fun = data.getFTensor0N();
85 
87  auto get_ftensor1 = [](auto &v) {
88  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
89  &v[2]);
90  };
91 
92  for (int gg = 0; gg != nb_integration_pts; ++gg) {
93  double a = v * t_w;
94  auto t_nf = get_ftensor1(nF);
95  int bb = 0;
96  for (; bb != nb_dofs / 3; ++bb) {
97  t_nf(i) += a * t_row_base_fun * t_div_P(i);
98  ++t_nf;
99  ++t_row_base_fun;
100  }
101  for (; bb != nb_base_functions; ++bb)
102  ++t_row_base_fun;
103  ++t_w;
104  ++t_div_P;
105  }
106 
108 }
109 
110 MoFEMErrorCode OpSpatialRotation::integrate(EntData &data) {
112  int nb_dofs = data.getIndices().size();
113  int nb_integration_pts = data.getN().size1();
114  auto v = getVolume();
115  auto t_w = getFTensor0IntegrationWeight();
116  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
117  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
118  auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
119  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
120 
121  int nb_base_functions = data.getN().size2();
122  auto t_row_base_fun = data.getFTensor0N();
128  auto get_ftensor1 = [](auto &v) {
129  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
130  &v[2]);
131  };
132 
133  for (int gg = 0; gg != nb_integration_pts; ++gg) {
134  double a = v * t_w;
135  auto t_nf = get_ftensor1(nF);
136  // FTensor::Tensor2<double, 3, 3> t_PhT;
137  // t_PhT(i, m) = (t_x_grad(j, l) * t_approx_P(i, l)) * t_delta_h(m, j);
139  t_PRT(i, m) = (t_x_grad(j, l) * t_approx_P(i, l)) * t_R(m, j);
140  FTensor::Tensor1<double, 3> t_leviPRT;
141  t_leviPRT(k) = levi_civita(i, m, k) * t_PRT(i, m);
142  int bb = 0;
143  for (; bb != nb_dofs / 3; ++bb) {
144  t_nf(k) += a * t_row_base_fun * t_leviPRT(k);
145  ++t_nf;
146  ++t_row_base_fun;
147  }
148  for (; bb != nb_base_functions; ++bb)
149  ++t_row_base_fun;
150  ++t_w;
151  ++t_approx_P;
152  ++t_x_grad;
153  ++t_delta_h;
154  ++t_R;
155  }
157 }
158 
159 MoFEMErrorCode OpSpatialPhysical::integrate(EntData &data) {
161  int nb_dofs = data.getIndices().size();
162  int nb_integration_pts = data.getN().size1();
163  auto v = getVolume();
164  auto t_w = getFTensor0IntegrationWeight();
165  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
166  auto t_P = getFTensor2FromMat<3, 3>(*(dataAtPts->PAtPts));
167  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
168  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
169  auto t_dot_u =
170  getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachDotTensorAtPts));
175  auto get_ftensor2 = [](auto &v) {
177  &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
178  };
179 
180  int nb_base_functions = data.getN().size2();
181  auto t_row_base_fun = data.getFTensor0N();
182  for (int gg = 0; gg != nb_integration_pts; ++gg) {
183  double a = v * t_w;
184  auto t_nf = get_ftensor2(nF);
185 
187  t_xP(i, j) = t_x_grad(j, k) * t_P(i, k);
189  t_RTxP(i, j) = t_R(k, i) * t_x_grad(j, l) * t_approx_P(k, l);
190 
191  int bb = 0;
192  for (; bb != nb_dofs / 6; ++bb) {
193 
194  for (int ii : {0, 1, 2})
195  for (int jj : {0, 1, 2}) {
196  t_nf(ii, jj) += a * t_row_base_fun * (t_RTxP(ii, jj) - t_xP(ii, jj));
197  t_nf(ii, jj) -= a * t_row_base_fun * alpha_u * t_dot_u(ii, jj);
198  }
199 
200  ++t_nf;
201  ++t_row_base_fun;
202  }
203  for (; bb != nb_base_functions; ++bb)
204  ++t_row_base_fun;
205 
206  ++t_w;
207  ++t_approx_P;
208  ++t_P;
209  ++t_R;
210  ++t_x_grad;
211  ++t_dot_u;
212  }
214 }
215 
216 MoFEMErrorCode OpSpatialConsistencyP::integrate(EntData &data) {
218  int nb_dofs = data.getIndices().size();
219  int nb_integration_pts = data.getN().size1();
220  auto v = getVolume();
221  auto t_w = getFTensor0IntegrationWeight();
222  auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
223  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
224 
225  int nb_base_functions = data.getN().size2() / 3;
226  auto t_row_base_fun = data.getFTensor1N<3>();
230  auto get_ftensor1 = [](auto &v) {
231  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
232  &v[2]);
233  };
234 
235  for (int gg = 0; gg != nb_integration_pts; ++gg) {
236  double a = v * t_w;
237  auto t_nf = get_ftensor1(nF);
238  int bb = 0;
239  for (; bb != nb_dofs / 3; ++bb) {
240  FTensor::Tensor1<double, 3> t_pushed_base;
241  t_pushed_base(i) = t_x_grad(i, k) * t_row_base_fun(k);
242  t_nf(i) += a * t_pushed_base(j) * t_delta_h(i, j);
243  t_nf(i) -= a * t_pushed_base(i);
244  ++t_nf;
245  ++t_row_base_fun;
246  }
247  for (; bb != nb_base_functions; ++bb) {
248  ++t_row_base_fun;
249  }
250  ++t_w;
251  ++t_delta_h;
252  ++t_x_grad;
253  }
254 
256 }
257 
258 MoFEMErrorCode OpSpatialConsistencyBubble::integrate(EntData &data) {
260  int nb_dofs = data.getIndices().size();
261  int nb_integration_pts = data.getN().size1();
262  auto v = getVolume();
263  auto t_w = getFTensor0IntegrationWeight();
264  auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
265  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
266 
267  int nb_base_functions = data.getN().size2() / 9;
268  auto t_row_base_fun = data.getFTensor2N<3, 3>();
272  auto get_ftensor0 = [](auto &v) {
274  };
275 
276  for (int gg = 0; gg != nb_integration_pts; ++gg) {
277  double a = v * t_w;
278  auto t_nf = get_ftensor0(nF);
279  int bb = 0;
280  for (; bb != nb_dofs; ++bb) {
281  FTensor::Tensor2<double, 3, 3> t_pushed_base;
282  t_pushed_base(i, j) = t_x_grad(j, k) * t_row_base_fun(i, k);
283  t_nf += a * t_pushed_base(i, j) * t_delta_h(i, j);
284  t_nf -= a * t_pushed_base(i, i);
285  ++t_nf;
286  ++t_row_base_fun;
287  }
288  for (; bb != nb_base_functions; ++bb) {
289  ++t_row_base_fun;
290  }
291  ++t_w;
292  ++t_delta_h;
293  ++t_x_grad;
294  }
295 
297 }
298 
299 MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
301  int nb_dofs = data.getIndices().size();
302  int nb_integration_pts = data.getN().size1();
303  auto v = getVolume();
304  auto t_w = getFTensor0IntegrationWeight();
305  auto t_s_w = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
306  int nb_base_functions = data.getN().size2() / 3;
307  auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
309  auto get_ftensor1 = [](auto &v) {
310  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
311  &v[2]);
312  };
313 
314  for (int gg = 0; gg != nb_integration_pts; ++gg) {
315  double a = v * t_w;
316  auto t_nf = get_ftensor1(nF);
317  int bb = 0;
318  for (; bb != nb_dofs / 3; ++bb) {
319  double div_row_base = t_row_diff_base_fun(i, i);
320  t_nf(i) += a * div_row_base * t_s_w(i);
321  ++t_nf;
322  ++t_row_diff_base_fun;
323  }
324  for (; bb != nb_base_functions; ++bb) {
325  ++t_row_diff_base_fun;
326  }
327  ++t_w;
328  ++t_s_w;
329  }
330 
332 }
333 
334 MoFEMErrorCode OpDispBc::integrate(EntData &data) {
336  // get entity of face
337  EntityHandle fe_ent = getFEEntityHandle();
338  // interate over all boundary data
339  for (auto &bc : (*bcDispPtr)) {
340  // check if finite element entity is part of boundary condition
341  if (bc.faces.find(fe_ent) != bc.faces.end()) {
342  int nb_dofs = data.getIndices().size();
343  int nb_integration_pts = data.getN().size1();
344  auto t_normal = getFTensor1Normal();
345  auto t_w = getFTensor0IntegrationWeight();
346  int nb_base_functions = data.getN().size2() / 3;
347  auto t_row_base_fun = data.getFTensor1N<3>();
348 
349  auto t_coords = getFTensor1CoordsAtGaussPts();
350  auto t_x = getFTensor1FromMat<3>(*(dataAtPts->xAtPts));
351 
354  auto get_ftensor1 = [](auto &v) {
355  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
356  &v[2]);
357  };
358 
359  // get bc data
360  FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
361  t_bc_disp(i) *= getFEMethod()->ts_t;
362 
363  for (int gg = 0; gg != nb_integration_pts; ++gg) {
365  t_bc_res(i) = t_bc_disp(i) - t_x(i) + t_coords(i);
366  auto t_nf = get_ftensor1(nF);
367  int bb = 0;
368  for (; bb != nb_dofs / 3; ++bb) {
369  t_nf(i) -=
370  t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
371  ++t_nf;
372  ++t_row_base_fun;
373  }
374  for (; bb != nb_base_functions; ++bb)
375  ++t_row_base_fun;
376 
377  ++t_w;
378  ++t_coords;
379  ++t_x;
380  }
381  }
382  }
384 }
385 
386 MoFEMErrorCode OpDispBc_dx::integrate(EntData &row_data, EntData &col_data) {
388 
389  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
391  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
392 
393  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
394 
395  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
396  };
397 
398  // get entity of face
399  EntityHandle fe_ent = getFEEntityHandle();
400  // interate over all boundary data
401  for (auto &bc : (*bcDispPtr)) {
402  // check if finite element entity is part of boundary condition
403  if (bc.faces.find(fe_ent) != bc.faces.end()) {
404 
405  int nb_integration_pts = row_data.getN().size1();
406  int row_nb_dofs = row_data.getIndices().size();
407  int col_nb_dofs = col_data.getIndices().size();
408 
409  auto t_normal = getFTensor1Normal();
410  auto t_w = getFTensor0IntegrationWeight();
411  int nb_base_functions = row_data.getN().size2() / 3;
412  auto t_row_base_fun = row_data.getFTensor1N<3>();
413 
416 
417  for (int gg = 0; gg != nb_integration_pts; ++gg) {
418 
419  int rr = 0;
420  for (; rr != row_nb_dofs / 3; ++rr) {
421 
422  auto t_m = get_ftensor2(K, 3 * rr, 0);
423  const double val = t_w * (t_row_base_fun(i) * t_normal(i)) * 0.5;
424  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
425 
426  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
427  const double set_val = val * t_col_base_fun;
428  t_m(0, 0) += set_val;
429  t_m(1, 1) += set_val;
430  t_m(2, 2) += set_val;
431  ++t_m;
432  ++t_col_base_fun;
433  }
434 
435  ++t_row_base_fun;
436  }
437  for (; rr != nb_base_functions; ++rr)
438  ++t_row_base_fun;
439 
440  ++t_w;
441  }
442  }
443  }
444 
446 }
447 
448 MoFEMErrorCode OpRotationBc::integrate(EntData &data) {
450  // get entity of face
451  EntityHandle fe_ent = getFEEntityHandle();
452  // interate over all boundary data
453  for (auto &bc : (*bcRotPtr)) {
454  // check if finite element entity is part of boundary condition
455  if (bc.faces.find(fe_ent) != bc.faces.end()) {
456  int nb_dofs = data.getIndices().size();
457  int nb_integration_pts = data.getN().size1();
458  auto t_normal = getFTensor1Normal();
459  auto t_w = getFTensor0IntegrationWeight();
460 
461  int nb_base_functions = data.getN().size2() / 3;
462  auto t_row_base_fun = data.getFTensor1N<3>();
466  auto get_ftensor1 = [](auto &v) {
467  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
468  &v[2]);
469  };
470 
471  // get bc data
472  FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
473 
474  double theta = bc.theta;
475  theta *= getFEMethod()->ts_t;
476 
478  const double a = sqrt(t_normal(i) * t_normal(i));
479  t_omega(i) = t_normal(i) * (theta / a);
480  auto t_R = getRotationFormVector(t_omega);
481  auto t_coords = getFTensor1CoordsAtGaussPts();
482  auto t_x = getFTensor1FromMat<3>(*(dataAtPts->xAtPts));
483 
484  for (int gg = 0; gg != nb_integration_pts; ++gg) {
486  t_delta(i) = t_center(i) - t_coords(i);
488  t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
489  t_disp(i) -= t_x(i) - t_coords(i);
490 
491  auto t_nf = get_ftensor1(nF);
492  int bb = 0;
493  for (; bb != nb_dofs / 3; ++bb) {
494  t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
495  ++t_nf;
496  ++t_row_base_fun;
497  }
498  for (; bb != nb_base_functions; ++bb)
499  ++t_row_base_fun;
500 
501  ++t_w;
502  ++t_coords;
503  ++t_x;
504  }
505  }
506  }
508 }
509 
510 MoFEMErrorCode OpRotationBc_dx::integrate(EntData &row_data,
511  EntData &col_data) {
513 
514  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
516  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
517 
518  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
519 
520  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
521  };
522 
523  // get entity of face
524  EntityHandle fe_ent = getFEEntityHandle();
525  // interate over all boundary data
526  for (auto &bc : (*bcRotPtr)) {
527  // check if finite element entity is part of boundary condition
528  if (bc.faces.find(fe_ent) != bc.faces.end()) {
529 
530  int nb_integration_pts = row_data.getN().size1();
531  int row_nb_dofs = row_data.getIndices().size();
532  int col_nb_dofs = col_data.getIndices().size();
533 
534  int row_nb_base_functions = row_data.getN().size2() / 3;
535  auto t_row_base_fun = row_data.getFTensor1N<3>();
536 
538 
539  auto t_normal = getFTensor1Normal();
540  auto t_w = getFTensor0IntegrationWeight();
541 
542  // get bc data
543  for (int gg = 0; gg != nb_integration_pts; ++gg) {
544 
545  int rr = 0;
546  for (; rr != row_nb_dofs / 3; ++rr) {
547 
548  const double val = t_w * (t_row_base_fun(i) * t_normal(i)) * 0.5;
549 
550  auto t_m = get_ftensor2(K, 3 * rr, 0);
551  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
552 
553  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
554  const double set_val = val * t_col_base_fun;
555  t_m(0, 0) += set_val;
556  t_m(1, 1) += set_val;
557  t_m(2, 2) += set_val;
558  ++t_m;
559  ++t_col_base_fun;
560  }
561 
562  ++t_row_base_fun;
563  }
564  for (; rr != row_nb_base_functions; ++rr)
565  ++t_row_base_fun;
566 
567  ++t_w;
568  }
569  }
570  }
571 
573 }
574 
575 MoFEMErrorCode FeTractionBc::preProcess() {
577 
578  struct FaceRule {
579  int operator()(int p_row, int p_col, int p_data) const {
580  return 2 * (p_data + 1);
581  }
582  };
583 
584  if (ts_ctx == CTX_TSSETIFUNCTION) {
585 
586  // Loop boundary elements with traction boundary conditions
588  fe.getOpPtrVector().push_back(
589  new OpTractionBc(std::string("P")/* + "_RT"*/, *this));
590  fe.getRuleHook = FaceRule();
591  fe.ts_t = ts_t;
592  CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
593  fe);
594  }
595 
597 }
598 
599 MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type, EntData &data) {
601 
604 
605  auto t_normal = getFTensor1Normal();
606  const double nrm2 = sqrt(t_normal(i) * t_normal(i));
607  FTensor::Tensor1<double, 3> t_unit_normal;
608  t_unit_normal(i) = t_normal(i) / nrm2;
609  int nb_dofs = data.getFieldData().size();
610  int nb_integration_pts = data.getN().size1();
611  int nb_base_functions = data.getN().size2() / 3;
612  double ts_t = getFEMethod()->ts_t;
613 
614  auto integrate_matrix = [&]() {
616 
617  auto t_w = getFTensor0IntegrationWeight();
618  matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
619  matPP.clear();
620 
621  auto t_row_base_fun = data.getFTensor1N<3>();
622  for (int gg = 0; gg != nb_integration_pts; ++gg) {
623 
624  int rr = 0;
625  for (; rr != nb_dofs / 3; ++rr) {
626  const double a = t_row_base_fun(i) * t_unit_normal(i);
627  auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
628  for (int cc = 0; cc != nb_dofs / 3; ++cc) {
629  const double b = t_col_base_fun(i) * t_unit_normal(i);
630  matPP(rr, cc) += t_w * a * b;
631  ++t_col_base_fun;
632  }
633  ++t_row_base_fun;
634  }
635 
636  for (; rr != nb_base_functions; ++rr)
637  ++t_row_base_fun;
638 
639  ++t_w;
640  }
641 
643  };
644 
645  auto integrate_rhs = [&](auto &bc) {
647 
648  auto t_w = getFTensor0IntegrationWeight();
649  vecPv.resize(3, nb_dofs / 3, false);
650  vecPv.clear();
651 
652  auto t_row_base_fun = data.getFTensor1N<3>();
653  double ts_t = getFEMethod()->ts_t;
654 
655  for (int gg = 0; gg != nb_integration_pts; ++gg) {
656  int rr = 0;
657  for (; rr != nb_dofs / 3; ++rr) {
658  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
659  for (int dd = 0; dd != 3; ++dd)
660  if (bc.flags[dd])
661  vecPv(dd, rr) += t * bc.vals[dd];
662  ++t_row_base_fun;
663  }
664  for (; rr != nb_base_functions; ++rr)
665  ++t_row_base_fun;
666  ++t_w;
667  }
669  };
670 
671  auto integrate_rhs_cook = [&](auto &bc) {
673 
674  vecPv.resize(3, nb_dofs / 3, false);
675  vecPv.clear();
676 
677  auto t_w = getFTensor0IntegrationWeight();
678  auto t_row_base_fun = data.getFTensor1N<3>();
679  auto t_coords = getFTensor1CoordsAtGaussPts();
680 
681  for (int gg = 0; gg != nb_integration_pts; ++gg) {
682 
683  auto calc_tau = [](double y) {
684  y -= 44;
685  y /= (60 - 44);
686  return -y * (y - 1) / 0.25;
687  };
688 
689  const double tau = calc_tau(t_coords(1));
690 
691  int rr = 0;
692  for (; rr != nb_dofs / 3; ++rr) {
693  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
694  for (int dd = 0; dd != 3; ++dd)
695  if (bc.flags[dd])
696  vecPv(dd, rr) += t * tau * bc.vals[dd];
697  ++t_row_base_fun;
698  }
699 
700  for (; rr != nb_base_functions; ++rr)
701  ++t_row_base_fun;
702  ++t_w;
703  ++t_coords;
704  }
706  };
707 
708  // get entity of face
709  EntityHandle fe_ent = getFEEntityHandle();
710  for (auto &bc : *(bcFe.bcData)) {
711  if (bc.faces.find(fe_ent) != bc.faces.end()) {
712 
713  int nb_dofs = data.getFieldData().size();
714  if (nb_dofs) {
715 
716  CHKERR integrate_matrix();
717  if (bc.blockName.compare(20, 4, "COOK") == 0)
718  CHKERR integrate_rhs_cook(bc);
719  else
720  CHKERR integrate_rhs(bc);
721 
722  const auto info =
723  lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
724  nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
725  if (info > 0)
726  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
727  "The leading minor of order %i is "
728  "not positive; definite;\nthe "
729  "solution could not be computed",
730  info);
731 
732  for (int dd = 0; dd != 3; ++dd)
733  if (bc.flags[dd])
734  for (int rr = 0; rr != nb_dofs / 3; ++rr)
735  data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
736 
737  // int nb_integration_pts = data.getN().size1();
738  // auto t_row_base_fun = data.getFTensor1N<3>();
739  // auto t_coords = getFTensor1CoordsAtGaussPts();
740  // for (int gg = 0; gg != nb_integration_pts; ++gg) {
741 
742  // FTensor::Tensor1<double, 3> t_v;
743  // t_v(i) = 0;
744 
745  // int rr = 0;
746  // for (; rr != nb_dofs / 3; ++rr) {
747  // for (int dd = 0; dd != 3; ++dd) {
748  // const double v = data.getFieldDofs()[3 * rr +
749  // dd]->getFieldData(); t_v(dd) += t_row_base_fun(i) *
750  // t_unit_normal(i) * v;
751  // }
752  // ++t_row_base_fun;
753  // }
754  // for (; rr != nb_base_functions; ++rr)
755  // ++t_row_base_fun;
756  // ++t_coords;
757 
758  // cerr << t_v << endl;
759  // }
760  }
761  }
762  }
763 
765 }
766 
767 MoFEMErrorCode OpSpatialEquilibrium_dw_dP::integrate(EntData &row_data,
768  EntData &col_data) {
770  int nb_integration_pts = row_data.getN().size1();
771  int row_nb_dofs = row_data.getIndices().size();
772  int col_nb_dofs = col_data.getIndices().size();
773  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
775  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
776  };
778  auto v = getVolume();
779  auto t_w = getFTensor0IntegrationWeight();
780  int row_nb_base_functions = row_data.getN().size2();
781  auto t_row_base_fun = row_data.getFTensor0N();
782  for (int gg = 0; gg != nb_integration_pts; ++gg) {
783  double a = v * t_w;
784  int rr = 0;
785  for (; rr != row_nb_dofs / 3; ++rr) {
786  auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
787  auto t_m = get_ftensor1(K, 3 * rr, 0);
788  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
789  double div_col_base = t_col_diff_base_fun(i, i);
790  t_m(i) += a * t_row_base_fun * div_col_base;
791  ++t_m;
792  ++t_col_diff_base_fun;
793  }
794  ++t_row_base_fun;
795  }
796  for (; rr != row_nb_base_functions; ++rr)
797  ++t_row_base_fun;
798  ++t_w;
799  }
801 }
802 
803 MoFEMErrorCode OpSpatialPhysical_du_dP::integrate(EntData &row_data,
804  EntData &col_data) {
806  int nb_integration_pts = row_data.getN().size1();
807  int row_nb_dofs = row_data.getIndices().size();
808  int col_nb_dofs = col_data.getIndices().size();
809  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
811 
812  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
813 
814  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
815 
816  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
817 
818  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
819 
820  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
821 
822  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
823  };
824 
828  auto v = getVolume();
829  auto t_w = getFTensor0IntegrationWeight();
830  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
831  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
832 
833  int row_nb_base_functions = row_data.getN().size2();
834  auto t_row_base_fun = row_data.getFTensor0N();
835  for (int gg = 0; gg != nb_integration_pts; ++gg) {
836  double a = v * t_w;
837 
838  int rr = 0;
839  for (; rr != row_nb_dofs / 6; ++rr) {
840 
841  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
842  auto t_m = get_ftensor3(K, 6 * rr, 0);
843 
844  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
845 
846  FTensor::Tensor1<double, 3> t_pushed_base;
847  t_pushed_base(i) = t_x_grad(i, k) * t_col_base_fun(k);
849  t_RTP_dP(i, j, k) = t_R(k, i) * t_pushed_base(j) * t_row_base_fun;
850 
851  for (int kk = 0; kk != 3; ++kk)
852  for (int ii = 0; ii != 3; ++ii)
853  for (int jj = 0; jj != 3; ++jj) {
854  t_m(ii, jj, kk) += a * t_RTP_dP(ii, jj, kk);
855  }
856 
857  ++t_col_base_fun;
858  ++t_m;
859  }
860 
861  ++t_row_base_fun;
862  }
863  for (; rr != row_nb_base_functions; ++rr)
864  ++t_row_base_fun;
865  ++t_w;
866  ++t_R;
867  ++t_x_grad;
868  }
870 }
871 
872 MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
873  EntData &col_data) {
875  int nb_integration_pts = row_data.getN().size1();
876  int row_nb_dofs = row_data.getIndices().size();
877  int col_nb_dofs = col_data.getIndices().size();
878  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
880  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
881  &m(r + 5, c));
882  };
883 
887  auto v = getVolume();
888  auto t_w = getFTensor0IntegrationWeight();
889  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
890  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
891 
892  int row_nb_base_functions = row_data.getN().size2();
893  auto t_row_base_fun = row_data.getFTensor0N();
894  for (int gg = 0; gg != nb_integration_pts; ++gg) {
895  double a = v * t_w;
896 
897  int rr = 0;
898  for (; rr != row_nb_dofs / 6; ++rr) {
899 
900  auto t_m = get_ftensor2(K, 6 * rr, 0);
901 
902  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
903  for (int cc = 0; cc != col_nb_dofs; ++cc) {
904 
905  FTensor::Tensor2<double, 3, 3> t_pushed_base;
906  t_pushed_base(i, j) = t_x_grad(j, k) * t_col_base_fun(i, k);
907  FTensor::Tensor2<double, 3, 3> t_RTP_dBubble;
908  t_RTP_dBubble(i, j) = t_R(k, i) * t_pushed_base(k, j) * t_row_base_fun;
909 
910  for (int ii = 0; ii != 3; ++ii)
911  for (int jj = 0; jj != 3; ++jj) {
912  t_m(ii, jj) += a * t_RTP_dBubble(ii, jj);
913  }
914 
915  ++t_m;
916  ++t_col_base_fun;
917  }
918 
919  ++t_row_base_fun;
920  }
921  for (; rr != row_nb_base_functions; ++rr)
922  ++t_row_base_fun;
923  ++t_w;
924  ++t_R;
925  ++t_x_grad;
926  }
928 }
929 
930 MoFEMErrorCode OpSpatialPhysical_du_du::integrate(EntData &row_data,
931  EntData &col_data) {
933 
934  int nb_integration_pts = row_data.getN().size1();
935  int row_nb_dofs = row_data.getIndices().size();
936  int col_nb_dofs = col_data.getIndices().size();
937  auto get_ftensor4 = [](MatrixDouble &m, const int r, const int c) {
939 
940  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
941  &m(r + 0, c + 4), &m(r + 0, c + 5),
942 
943  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
944  &m(r + 1, c + 4), &m(r + 1, c + 5),
945 
946  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
947  &m(r + 2, c + 4), &m(r + 2, c + 5),
948 
949  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
950  &m(r + 3, c + 4), &m(r + 3, c + 5),
951 
952  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
953  &m(r + 4, c + 4), &m(r + 4, c + 5),
954 
955  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
956  &m(r + 5, c + 4), &m(r + 5, c + 5)
957 
958  );
959  };
966 
970 
971  auto v = getVolume();
972  auto t_w = getFTensor0IntegrationWeight();
973  auto t_P_dh0 = getFTensor3FromMat(*(dataAtPts->P_dh0));
974  auto t_P_dh1 = getFTensor3FromMat(*(dataAtPts->P_dh1));
975  auto t_P_dh2 = getFTensor3FromMat(*(dataAtPts->P_dh2));
976  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
977  int row_nb_base_functions = row_data.getN().size2();
978  auto t_row_base_fun = row_data.getFTensor0N();
979 
980  const double ts_a = getFEMethod()->ts_a;
981  for (int gg = 0; gg != nb_integration_pts; ++gg) {
982  double a = v * t_w;
983 
985  t_P_dh(i, j, N0, k) = t_x_grad(j, m) * t_P_dh0(i, m, k);
986  t_P_dh(i, j, N1, k) = t_x_grad(j, m) * t_P_dh1(i, m, k);
987  t_P_dh(i, j, N2, k) = t_x_grad(j, m) * t_P_dh2(i, m, k);
988 
990  t_pushed_P_du(i, j, k, l) = t_P_dh(i, j, k, n) * t_x_grad(l, n);
991 
992  int rr = 0;
993  for (; rr != row_nb_dofs / 6; ++rr) {
994 
995  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
996  auto t_m = get_ftensor4(K, 6 * rr, 0);
997 
998  for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
999 
1000  const double b = a * t_row_base_fun * t_col_base_fun;
1001 
1002  for (int ii : {0, 1, 2})
1003  for (int jj : {0, 1, 2})
1004  for (int kk : {0, 1, 2})
1005  for (int ll : {0, 1, 2})
1006  t_m(ii, jj, kk, ll) -= b * t_pushed_P_du(ii, jj, kk, ll);
1007 
1008  for (int ii : {0, 1, 2})
1009  for (int jj : {0, 1, 2})
1010  t_m(ii, jj, ii, jj) -= b * alpha_u * ts_a;
1011 
1012  ++t_m;
1013  ++t_col_base_fun;
1014  }
1015 
1016  ++t_row_base_fun;
1017  }
1018  for (; rr != row_nb_base_functions; ++rr)
1019  ++t_row_base_fun;
1020 
1021  ++t_w;
1022  ++t_P_dh0;
1023  ++t_P_dh1;
1024  ++t_P_dh2;
1025  ++t_x_grad;
1026  }
1028 }
1029 
1030 MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
1031  EntData &col_data) {
1033  int row_nb_dofs = row_data.getIndices().size();
1034  int col_nb_dofs = col_data.getIndices().size();
1035  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1037 
1038  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1039 
1040  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1041 
1042  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1043 
1044  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1045 
1046  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1047 
1048  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1049 
1050  );
1051  };
1058 
1059  FTensor::Number<0> N0;
1060  FTensor::Number<1> N1;
1061  FTensor::Number<2> N2;
1062 
1063  auto v = getVolume();
1064  auto t_w = getFTensor0IntegrationWeight();
1065  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
1066  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
1067  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
1068 
1069  int row_nb_base_functions = row_data.getN().size2();
1070  auto t_row_base_fun = row_data.getFTensor0N();
1071 
1072  int nb_integration_pts = row_data.getN().size1();
1073 
1074  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1075  double a = v * t_w;
1076 
1078  t_P(i, j) = t_x_grad(j, k) * t_approx_P(i, k);
1079  FTensor::Tensor3<double, 3, 3, 3> t_RT_domega_P;
1080  t_RT_domega_P(i, j, m) = t_diff_R(k, i, m) * t_P(k, j);
1081 
1082  int rr = 0;
1083  for (; rr != row_nb_dofs / 6; ++rr) {
1084  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1085  auto t_m = get_ftensor3(K, 6 * rr, 0);
1086  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1087 
1088  double v = a * t_row_base_fun * t_col_base_fun;
1089 
1090  for (int ii = 0; ii != 3; ++ii)
1091  for (int jj = 0; jj != 3; ++jj)
1092  for (int kk = 0; kk != 3; ++kk)
1093  t_m(ii, jj, kk) += v * t_RT_domega_P(ii, jj, kk);
1094 
1095  ++t_m;
1096  ++t_col_base_fun;
1097  }
1098 
1099  ++t_row_base_fun;
1100  }
1101 
1102  for (; rr != row_nb_base_functions; ++rr)
1103  ++t_row_base_fun;
1104 
1105  ++t_w;
1106  ++t_approx_P;
1107  ++t_diff_R;
1108  ++t_x_grad;
1109  }
1110 
1112 }
1113 
1114 MoFEMErrorCode OpSpatialPhysical_du_dx::integrate(EntData &row_data,
1115  EntData &col_data) {
1117 
1118  int row_nb_dofs = row_data.getIndices().size();
1119  int col_nb_dofs = col_data.getIndices().size();
1120 
1121  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1123 
1124  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1125 
1126  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1127 
1128  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
1129 
1130  &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
1131 
1132  &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
1133 
1134  &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
1135 
1136  );
1137  };
1138 
1139  auto v = getVolume();
1140  auto t_w = getFTensor0IntegrationWeight();
1141  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
1142  auto t_P = getFTensor2FromMat<3, 3>(*(dataAtPts->PAtPts));
1143  auto t_P_dh0 = getFTensor3FromMat(*(dataAtPts->P_dh0));
1144  auto t_P_dh1 = getFTensor3FromMat(*(dataAtPts->P_dh1));
1145  auto t_P_dh2 = getFTensor3FromMat(*(dataAtPts->P_dh2));
1146  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1147  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
1148  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
1149 
1156 
1157  FTensor::Number<0> N0;
1158  FTensor::Number<1> N1;
1159  FTensor::Number<2> N2;
1160 
1161  int nb_integration_pts = row_data.getN().size1();
1162  int row_nb_base_functions = row_data.getN().size2();
1163  auto t_row_base_fun = row_data.getFTensor0N();
1164 
1165  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1166  double a = v * t_w;
1167 
1169  t_delta_P(i, j) = t_R(l, i) * t_approx_P(l, j) - t_P(i, j);
1170 
1172  t_P_dh(i, j, N0, l) = t_x_grad(j, m) * t_P_dh0(i, m, l);
1173  t_P_dh(i, j, N1, l) = t_x_grad(j, m) * t_P_dh1(i, m, l);
1174  t_P_dh(i, j, N2, l) = t_x_grad(j, m) * t_P_dh2(i, m, l);
1175 
1177  for (int ii = 0; ii != 3; ++ii)
1178  for (int jj = 0; jj != 3; ++jj)
1179  t_delta_h(ii, jj) = t_u(ii, jj);
1180 
1181  for (int ii = 0; ii != 3; ++ii)
1182  t_delta_h(ii, ii) += 1;
1183 
1184  FTensor::Tensor4<double, 3, 3, 3, 3> t_P_intermidiate_dh;
1185  t_P_intermidiate_dh(i, j, k, l) = t_P_dh(i, j, m, l) * t_delta_h(m, k);
1186 
1187  int rr = 0;
1188  for (; rr != row_nb_dofs / 6; ++rr) {
1189 
1190  auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
1191  auto t_m = get_ftensor3(K, 6 * rr, 0);
1192  const double b = a * t_row_base_fun;
1193 
1194  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1195 
1196  FTensor::Tensor1<double, 3> t_dgrad_x_deltaP;
1197  t_dgrad_x_deltaP(i) = b * t_col_diff_base_fun(k) * t_delta_P(i, k);
1198 
1199  FTensor::Tensor3<double, 3, 3, 3> t_col_P_delta_dh;
1200  t_col_P_delta_dh(i, j, k) =
1201  b * (t_P_intermidiate_dh(i, j, k, l) * t_col_diff_base_fun(l));
1202 
1203  for (int ii = 0; ii != 3; ++ii)
1204  for (int jj = 0; jj != 3; ++jj) {
1205  t_m(ii, jj, jj) += t_dgrad_x_deltaP(ii);
1206  for (int kk = 0; kk != 3; ++kk)
1207  t_m(ii, jj, kk) -= t_col_P_delta_dh(ii, jj, kk);
1208  }
1209 
1210  ++t_m;
1211  ++t_col_diff_base_fun;
1212  }
1213 
1214  ++t_row_base_fun;
1215  }
1216 
1217  for (; rr != row_nb_base_functions; ++rr)
1218  ++t_row_base_fun;
1219 
1220  ++t_approx_P;
1221  ++t_P;
1222  ++t_P_dh0;
1223  ++t_P_dh1;
1224  ++t_P_dh2;
1225  ++t_u;
1226  ++t_x_grad;
1227  ++t_R;
1228  ++t_w;
1229  }
1230 
1232 }
1233 
1234 MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
1235  EntData &col_data) {
1237  int nb_integration_pts = row_data.getN().size1();
1238  int row_nb_dofs = row_data.getIndices().size();
1239  int col_nb_dofs = col_data.getIndices().size();
1240  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1242 
1243  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1244 
1245  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1246 
1247  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1248 
1249  );
1250  };
1257 
1258  auto v = getVolume();
1259  auto t_w = getFTensor0IntegrationWeight();
1260  auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
1261  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1262  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
1263 
1264  int row_nb_base_functions = row_data.getN().size2();
1265  auto t_row_base_fun = row_data.getFTensor0N();
1266 
1267  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1268  double a = v * t_w;
1269 
1270  int rr = 0;
1271  for (; rr != row_nb_dofs / 3; ++rr) {
1272  auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
1273  auto t_m = get_ftensor2(K, 3 * rr, 0);
1275  // t_levi_h(i, j, k) =
1276  // (a * t_row_base_fun) * (levi_civita(i, m, k) * t_delta_h(m, j));
1277  t_levi_h(i, j, k) =
1278  (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1279  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1280  FTensor::Tensor1<double, 3> t_pushed_base;
1281  t_pushed_base(i) = t_x_grad(i, k) * t_col_base_fun(k);
1282  t_m(k, i) += t_levi_h(i, j, k) * t_pushed_base(j);
1283  ++t_m;
1284  ++t_col_base_fun;
1285  }
1286 
1287  ++t_row_base_fun;
1288  }
1289 
1290  for (; rr != row_nb_base_functions; ++rr)
1291  ++t_row_base_fun;
1292  ++t_w;
1293  ++t_delta_h;
1294  ++t_R;
1295  ++t_x_grad;
1296  }
1298 }
1299 
1300 MoFEMErrorCode OpSpatialRotation_domega_dBubble::integrate(EntData &row_data,
1301  EntData &col_data) {
1303  int nb_integration_pts = row_data.getN().size1();
1304  int row_nb_dofs = row_data.getIndices().size();
1305  int col_nb_dofs = col_data.getIndices().size();
1306  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
1308  &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
1309  };
1315 
1316  auto v = getVolume();
1317  auto t_w = getFTensor0IntegrationWeight();
1318  auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
1319  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1320  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
1321 
1322  int row_nb_base_functions = row_data.getN().size2();
1323  auto t_row_base_fun = row_data.getFTensor0N();
1324 
1325  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1326  double a = v * t_w;
1327 
1328  int rr = 0;
1329  for (; rr != row_nb_dofs / 3; ++rr) {
1330 
1331  auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
1332  auto t_m = get_ftensor1(K, 3 * rr, 0);
1333 
1335  // t_levi_h(i, j, k) =
1336  // (a * t_row_base_fun) * (levi_civita(i, m, k) * t_delta_h(m, j));
1337  t_levi_h(i, j, k) =
1338  (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
1339 
1340  for (int cc = 0; cc != col_nb_dofs; ++cc) {
1341  FTensor::Tensor2<double, 3, 3> t_pushed_base;
1342  t_pushed_base(i, j) = t_x_grad(j, k) * t_col_base_fun(i, k);
1343  t_m(k) += t_levi_h(i, j, k) * t_pushed_base(i, j);
1344  ++t_m;
1345  ++t_col_base_fun;
1346  }
1347 
1348  ++t_row_base_fun;
1349  }
1350 
1351  for (; rr != row_nb_base_functions; ++rr)
1352  ++t_row_base_fun;
1353  ++t_w;
1354  ++t_delta_h;
1355  ++t_R;
1356  ++t_x_grad;
1357  }
1359 }
1360 
1361 MoFEMErrorCode OpSpatialRotation_domega_domega::integrate(EntData &row_data,
1362  EntData &col_data) {
1364  int nb_integration_pts = row_data.getN().size1();
1365  int row_nb_dofs = row_data.getIndices().size();
1366  int col_nb_dofs = col_data.getIndices().size();
1367  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1369 
1370  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1371 
1372  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1373 
1374  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1375 
1376  );
1377  };
1384 
1385  auto v = getVolume();
1386  auto t_w = getFTensor0IntegrationWeight();
1387  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
1388  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1389  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
1390  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
1391  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
1392 
1393  int row_nb_base_functions = row_data.getN().size2();
1394  auto t_row_base_fun = row_data.getFTensor0N();
1395 
1396  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1397  double a = v * t_w;
1398 
1399  FTensor::Tensor2<double, 3, 3> t_pushed_P;
1400  t_pushed_P(i, j) = t_x_grad(j, k) * t_approx_P(i, k);
1401 
1402  // FTensor::Tensor2<double, 3, 3> t_leviPhT_domega;
1403  // t_leviPhT_domega(n, l) =
1404  // levi_civita(i, j, n) *
1405  // (t_pushed_P(i, k) *
1406 
1407  // (t_diff_R(j, m, l) * t_u(m, k) + t_diff_R(j, k, l))
1408 
1409  // );
1410 
1411  FTensor::Tensor2<double, 3, 3> t_leviPRT_domega;
1412  t_leviPRT_domega(n, l) =
1413  levi_civita(i, j, n) * (t_pushed_P(i, k) * t_diff_R(j, k, l));
1414 
1415  int rr = 0;
1416  for (; rr != row_nb_dofs / 3; ++rr) {
1417 
1418  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1419  auto t_m = get_ftensor2(K, 3 * rr, 0);
1420  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1421  t_m(k, m) +=
1422  (a * t_row_base_fun * t_col_base_fun) * t_leviPRT_domega(k, m);
1423  ++t_m;
1424  ++t_col_base_fun;
1425  }
1426 
1427  ++t_row_base_fun;
1428  }
1429 
1430  for (; rr != row_nb_base_functions; ++rr)
1431  ++t_row_base_fun;
1432  ++t_w;
1433  ++t_approx_P;
1434  ++t_R;
1435  ++t_u;
1436  ++t_diff_R;
1437  ++t_x_grad;
1438  }
1440 }
1441 
1442 MoFEMErrorCode OpSpatialRotation_domega_du::integrate(EntData &row_data,
1443  EntData &col_data) {
1445  // int nb_integration_pts = row_data.getN().size1();
1446  // int row_nb_dofs = row_data.getIndices().size();
1447  // int col_nb_dofs = col_data.getIndices().size();
1448 
1449  // auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
1450  // return FTensor::Christof<FTensor::PackPtr<double *, 6>, 3, 3>(
1451 
1452  // &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
1453  // &m(r + 0, c + 4), &m(r + 0, c + 5),
1454 
1455  // &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
1456  // &m(r + 1, c + 4), &m(r + 1, c + 5),
1457 
1458  // &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
1459  // &m(r + 2, c + 4), &m(r + 2, c + 5)
1460 
1461  // );
1462  // };
1463 
1464  // FTensor::Index<'i', 3> i;
1465  // FTensor::Index<'j', 3> j;
1466  // FTensor::Index<'k', 3> k;
1467  // FTensor::Index<'l', 3> l;
1468  // FTensor::Index<'m', 3> m;
1469  // FTensor::Index<'n', 3> n;
1470 
1471  // auto v = getVolume();
1472  // auto t_w = getFTensor0IntegrationWeight();
1473  // auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
1474  // auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1475  // auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
1476  // auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
1477  // auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
1478 
1479  // int row_nb_base_functions = row_data.getN().size2();
1480  // auto t_row_base_fun = row_data.getFTensor0N();
1481 
1482  // for (int gg = 0; gg != nb_integration_pts; ++gg) {
1483  // double a = v * t_w;
1484 
1485  // FTensor::Tensor2<double, 3, 3> t_pushed_P;
1486  // t_pushed_P(i, j) = t_x_grad(j, k) * t_approx_P(i, k);
1487  // FTensor::Tensor4<double, 3, 3, 3, 3> t_pushed_PhT_du;
1488  // t_pushed_PhT_du(i, j, k, l) = t_pushed_P(i, l) * t_R(j, k);
1489 
1490  // FTensor::Tensor3<double, 3, 3, 3> t_levi;
1491  // t_levi(i, j, k) = levi_civita(i, j, k);
1492  // FTensor::Tensor3<double, 3, 3, 3> t_leviPhT_du;
1493  // t_leviPhT_du(n, k, l) = 0;
1494  // for (int ii : {0, 1, 2})
1495  // for (int jj : {0, 1, 2})
1496  // for (int nn : {0, 1, 2})
1497  // for (int kk : {0, 1, 2})
1498  // for (int ll : {0, 1, 2})
1499  // t_leviPhT_du(nn, kk, ll) += t_levi(ii, jj, nn) *
1500  // t_pushed_PhT_du(ii, jj, kk, ll);
1501 
1502 
1503  // int rr = 0;
1504  // for (; rr != row_nb_dofs / 3; ++rr) {
1505 
1506  // auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1507  // auto t_m = get_ftensor3(K, 3 * rr, 0);
1508  // for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1509 
1510  // const double b = a * t_row_base_fun * t_col_base_fun;
1511 
1512  // for (int ii : {0, 1, 2})
1513  // for (int jj : {0, 1, 2})
1514  // for (int kk : {0, 1, 2})
1515  // t_m(ii, jj, kk) += b * t_leviPhT_du(ii, jj, kk);
1516 
1517  // ++t_m;
1518  // ++t_col_base_fun;
1519  // }
1520 
1521  // ++t_row_base_fun;
1522  // }
1523 
1524  // for (; rr != row_nb_base_functions; ++rr)
1525  // ++t_row_base_fun;
1526  // ++t_w;
1527  // ++t_approx_P;
1528  // ++t_R;
1529  // ++t_u;
1530  // ++t_x_grad;
1531  // }
1533 }
1534 
1535 MoFEMErrorCode OpSpatialRotation_domega_dx::integrate(EntData &row_data,
1536  EntData &col_data) {
1538  int nb_integration_pts = row_data.getN().size1();
1539  int row_nb_dofs = row_data.getIndices().size();
1540  int col_nb_dofs = col_data.getIndices().size();
1541  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1543 
1544  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1545 
1546  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1547 
1548  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1549 
1550  );
1551  };
1552 
1558 
1559  auto v = getVolume();
1560  auto t_w = getFTensor0IntegrationWeight();
1561  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
1562  auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
1563  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
1564 
1565  int row_nb_base_functions = row_data.getN().size2();
1566  auto t_row_base_fun = row_data.getFTensor0N();
1567 
1568  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1569  double a = v * t_w;
1570 
1571  int rr = 0;
1572  for (; rr != row_nb_dofs / 3; ++rr) {
1573 
1574  // FTensor::Tensor3<double, 3, 3, 3> t_levihT;
1575  // t_levihT(i, l, k) = levi_civita(i, j, k) * t_delta_h(j, l);
1577  t_leviRT(i, l, k) = levi_civita(i, j, k) * t_R(j, l);
1578 
1579  auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
1580  auto t_m = get_ftensor2(K, 3 * rr, 0);
1581 
1582  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1583  t_m(k, l) +=
1584  a * t_row_base_fun *
1585  ((t_col_diff_base_fun(m) * t_approx_P(i, m)) * t_leviRT(i, l, k));
1586  ++t_col_diff_base_fun;
1587  ++t_m;
1588  }
1589 
1590  ++t_row_base_fun;
1591  }
1592  for (; rr != row_nb_base_functions; ++rr)
1593  ++t_row_base_fun;
1594 
1595  ++t_w;
1596  ++t_approx_P;
1597  ++t_delta_h;
1598  ++t_R;
1599  }
1601 }
1602 
1603 MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
1604  EntData &col_data) {
1606 
1607  int nb_integration_pts = row_data.getN().size1();
1608  int row_nb_dofs = row_data.getIndices().size();
1609  int col_nb_dofs = col_data.getIndices().size();
1610 
1611  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1613 
1614  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1615 
1616  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1617 
1618  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1619 
1620  );
1621  };
1622 
1627 
1628  auto v = getVolume();
1629  auto t_w = getFTensor0IntegrationWeight();
1630  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
1631  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
1632  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
1633  int row_nb_base_functions = row_data.getN().size2() / 3;
1634  auto t_row_base_fun = row_data.getFTensor1N<3>();
1635 
1636  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1637  double a = v * t_w;
1638 
1639  int rr = 0;
1640  for (; rr != row_nb_dofs / 3; ++rr) {
1641 
1642  FTensor::Tensor1<double, 3> t_pushed_base;
1643  t_pushed_base(i) = t_x_grad(i, k) * t_row_base_fun(k);
1644 
1646  t_PRT(i, k) = t_pushed_base(j) *
1647  (t_diff_R(i, l, k) * t_u(l, j) + t_diff_R(i, j, k));
1648 
1649  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1650  auto t_m = get_ftensor2(K, 3 * rr, 0);
1651  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1652  t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
1653  ++t_m;
1654  ++t_col_base_fun;
1655  }
1656 
1657  ++t_row_base_fun;
1658  }
1659 
1660  for (; rr != row_nb_base_functions; ++rr)
1661  ++t_row_base_fun;
1662 
1663  ++t_w;
1664  ++t_diff_R;
1665  ++t_u;
1666  ++t_x_grad;
1667  }
1669 }
1670 
1672 OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
1673  EntData &col_data) {
1675 
1676  int nb_integration_pts = row_data.getN().size1();
1677  int row_nb_dofs = row_data.getIndices().size();
1678  int col_nb_dofs = col_data.getIndices().size();
1679 
1680  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1682  &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
1683  };
1684 
1689 
1690  auto v = getVolume();
1691  auto t_w = getFTensor0IntegrationWeight();
1692  auto t_diff_R = getFTensor3FromMat(*(dataAtPts->diffRotMatAtPts));
1693  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
1694  auto t_x_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xGradAtPts));
1695  int row_nb_base_functions = row_data.getN().size2() / 9;
1696  auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
1697 
1698  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1699  double a = v * t_w;
1700 
1701  int rr = 0;
1702  for (; rr != row_nb_dofs; ++rr) {
1703 
1704  FTensor::Tensor2<double, 3, 3> t_pushed_base;
1705  t_pushed_base(i, j) = t_x_grad(j, k) * t_row_base_fun(i, k);
1706 
1708  t_PRT(k) = t_pushed_base(i, j) *
1709  (t_diff_R(i, l, k) * t_u(l, j) + t_diff_R(i, j, k));
1710 
1711  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1712  auto t_m = get_ftensor2(K, rr, 0);
1713  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1714  t_m(j) += (a * t_col_base_fun) * t_PRT(j);
1715  ++t_m;
1716  ++t_col_base_fun;
1717  }
1718 
1719  ++t_row_base_fun;
1720  }
1721 
1722  for (; rr != row_nb_base_functions; ++rr)
1723  ++t_row_base_fun;
1724 
1725  ++t_w;
1726  ++t_diff_R;
1727  ++t_x_grad;
1728  ++t_u;
1729  }
1731 }
1732 
1733 MoFEMErrorCode OpSpatialConsistency_dP_dx::integrate(EntData &row_data,
1734  EntData &col_data) {
1736 
1737  int nb_integration_pts = row_data.getN().size1();
1738  int row_nb_dofs = row_data.getIndices().size();
1739  int col_nb_dofs = col_data.getIndices().size();
1740 
1741  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1743 
1744  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1745 
1746  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1747 
1748  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1749 
1750  );
1751  };
1752 
1757 
1758  auto v = getVolume();
1759  auto t_w = getFTensor0IntegrationWeight();
1760  auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
1761  int row_nb_base_functions = row_data.getN().size2() / 3;
1762  auto t_row_base_fun = row_data.getFTensor1N<3>();
1763 
1764  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1765  double a = v * t_w;
1766 
1767  int rr = 0;
1768  for (; rr != row_nb_dofs / 3; ++rr) {
1769 
1770  auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
1771  auto t_m = get_ftensor2(K, 3 * rr, 0);
1772  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1773 
1774  const double b = a * t_col_diff_base_fun(l) * t_row_base_fun(l);
1775 
1776  t_m(i, j) += b * t_delta_h(i, j);
1777  t_m(0, 0) -= b;
1778  t_m(1, 1) -= b;
1779  t_m(2, 2) -= b;
1780 
1781  ++t_m;
1782  ++t_col_diff_base_fun;
1783  }
1784 
1785  ++t_row_base_fun;
1786  }
1787 
1788  for (; rr != row_nb_base_functions; ++rr)
1789  ++t_row_base_fun;
1790 
1791  ++t_w;
1792  ++t_delta_h;
1793  }
1795 }
1796 
1797 MoFEMErrorCode OpSpatialConsistency_dBubble_dx::integrate(EntData &row_data,
1798  EntData &col_data) {
1800 
1801  int nb_integration_pts = row_data.getN().size1();
1802  int row_nb_dofs = row_data.getIndices().size();
1803  int col_nb_dofs = col_data.getIndices().size();
1804 
1805  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1807  &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
1808  };
1809 
1814 
1815  auto v = getVolume();
1816  auto t_w = getFTensor0IntegrationWeight();
1817  auto t_delta_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hDeltaAtPts));
1818  int row_nb_base_functions = row_data.getN().size2() / 9;
1819  auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
1820 
1821  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1822  double a = v * t_w;
1823 
1824  int rr = 0;
1825  for (; rr != row_nb_dofs; ++rr) {
1826 
1827  auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
1828  auto t_m = get_ftensor2(K, rr, 0);
1829 
1830  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1831 
1832  FTensor::Tensor1<double, 3> t_pushed_base_dx;
1833  t_pushed_base_dx(i) = a * t_col_diff_base_fun(l) * t_row_base_fun(i, l);
1834  t_m(j) += t_pushed_base_dx(i) * t_delta_h(i, j);
1835  t_m(i) -= t_pushed_base_dx(i);
1836 
1837  ++t_m;
1838  ++t_col_diff_base_fun;
1839  }
1840 
1841  ++t_row_base_fun;
1842  }
1843 
1844  for (; rr != row_nb_base_functions; ++rr)
1845  ++t_row_base_fun;
1846 
1847  ++t_w;
1848  ++t_delta_h;
1849  }
1851 }
1852 
1853 MoFEMErrorCode OpSpatialPrj::integrate(EntData &row_data) {
1855 
1856  int nb_integration_pts = row_data.getN().size1();
1857  int row_nb_dofs = row_data.getIndices().size();
1858 
1859  auto get_ftensor1 = [](VectorDouble &v, const int r) {
1861  &v[r + 0], &v[r + 1], &v[r + 2]);
1862  };
1863 
1867 
1868  int row_nb_base_functions = row_data.getN().size2();
1869  auto t_row_base_fun = row_data.getFTensor0N();
1870 
1871  auto v = getVolume();
1872  auto t_w = getFTensor0IntegrationWeight();
1873 
1874  auto t_x_det = getFTensor0FromVec(*(dataAtPts->xDetAtPts));
1875  auto t_spatial_disp = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
1876 
1877  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1878 
1879  auto t_rhs = get_ftensor1(nF, 0);
1880 
1881  double a = v * t_w * t_x_det;
1882 
1883  int rr = 0;
1884  for (; rr != row_nb_dofs / 3; ++rr) {
1885  t_rhs(i) += a * t_row_base_fun * (-t_spatial_disp(i));
1886  ++t_row_base_fun;
1887  ++t_rhs;
1888  }
1889 
1890  for (; rr < row_nb_base_functions; ++rr)
1891  ++t_row_base_fun;
1892 
1893  ++t_spatial_disp;
1894  ++t_w;
1895  ++t_x_det;
1896  }
1897 
1899 }
1900 
1901 MoFEMErrorCode OpSpatialPrj_dx_dx::integrate(EntData &row_data,
1902  EntData &col_data) {
1904 
1905  int nb_integration_pts = row_data.getN().size1();
1906  int row_nb_dofs = row_data.getIndices().size();
1907  int col_nb_dofs = col_data.getIndices().size();
1908 
1909  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1911 
1912  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1913 
1914  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1915 
1916  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1917 
1918  );
1919  };
1920 
1927 
1928  auto t_x_det = getFTensor0FromVec(*(dataAtPts->xDetAtPts));
1929  auto t_x_inv_grad = getFTensor2FromMat<3, 3>(*(dataAtPts->xInvGradAtPts));
1930  auto t_spatial_disp = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
1931 
1932  int row_nb_base_functions = row_data.getN().size2();
1933  auto t_row_base_fun = row_data.getFTensor0N();
1934 
1935  auto v = getVolume();
1936  auto t_w = getFTensor0IntegrationWeight();
1937 
1938  const double ts_a = getFEMethod()->ts_a;
1939 
1940  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1941 
1942  double a = v * t_w * t_x_det;
1943 
1944  int rr = 0;
1945  for (; rr != row_nb_dofs / 3; ++rr) {
1946 
1947  auto t_m = get_ftensor2(K, 3 * rr, 0);
1948  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
1949  auto t_col_diff_base_fun = col_data.getFTensor1DiffN<3>(gg, 0);
1950 
1952  t_val(i) = a * t_row_base_fun * (-t_spatial_disp(i));
1953 
1954  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1955 
1956  t_m(i, j) += t_val(i) * t_x_inv_grad(k, j) * t_col_diff_base_fun(k);
1957 
1958  ++t_m;
1959  ++t_col_base_fun;
1960  ++t_col_diff_base_fun;
1961  }
1962 
1963  ++t_row_base_fun;
1964  }
1965  for (; rr < row_nb_base_functions; ++rr)
1966  ++t_row_base_fun;
1967 
1968  ++t_spatial_disp;
1969  ++t_w;
1970  ++t_x_det;
1971  ++t_x_inv_grad;
1972  }
1973 
1975 }
1976 
1977 MoFEMErrorCode OpSpatialPrj_dx_dw::integrate(EntData &row_data,
1978  EntData &col_data) {
1980 
1981  int nb_integration_pts = row_data.getN().size1();
1982  int row_nb_dofs = row_data.getIndices().size();
1983  int col_nb_dofs = col_data.getIndices().size();
1984 
1985  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
1987 
1988  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
1989 
1990  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
1991 
1992  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
1993 
1994  );
1995  };
1996 
1999 
2000  int row_nb_base_functions = row_data.getN().size2();
2001  auto t_row_base_fun = row_data.getFTensor0N();
2002 
2003  auto v = getVolume();
2004  auto t_w = getFTensor0IntegrationWeight();
2005  auto t_x_det = getFTensor0FromVec(*(dataAtPts->xDetAtPts));
2006 
2007  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2008 
2009  double a = v * t_w * t_x_det;
2010 
2011  int rr = 0;
2012  for (; rr != row_nb_dofs / 3; ++rr) {
2013 
2014  auto t_m = get_ftensor2(K, 3 * rr, 0);
2015  auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
2016 
2017  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2018 
2019  const double val = a * t_row_base_fun * t_col_base_fun;
2020  t_m(0, 0) -= val;
2021  t_m(1, 1) -= val;
2022  t_m(2, 2) -= val;
2023 
2024  ++t_m;
2025  ++t_col_base_fun;
2026  }
2027 
2028  ++t_row_base_fun;
2029  }
2030  for (; rr < row_nb_base_functions; ++rr)
2031  ++t_row_base_fun;
2032 
2033  ++t_w;
2034  ++t_x_det;
2035  }
2036 
2038 };
2039 
2040 MoFEMErrorCode OpSpatialSchurBegin::assemble(int row_side, EntityType row_type,
2041  EntData &data) {
2043  if (row_type != MBTET)
2045  auto &bmc = dataAtPts->blockMatContainor;
2046  for (auto &bit : bmc)
2047  bit.unSetAtElement();
2048  // bmc.clear();
2050 }
2051 
2052 MoFEMErrorCode OpSpatialPreconditionMass::integrate(EntData &row_data) {
2054  const int nb_integration_pts = row_data.getN().size1();
2055  const int row_nb_dofs = row_data.getIndices().size();
2056  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
2058  &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
2059 
2060  );
2061  };
2063 
2064  auto v = getVolume();
2065  auto t_w = getFTensor0IntegrationWeight();
2066  auto &m = *mPtr;
2067 
2068  m.resize(row_nb_dofs, row_nb_dofs, false);
2069  m.clear();
2070 
2071  int row_nb_base_functions = row_data.getN().size2();
2072  auto t_row_base_fun = row_data.getFTensor0N();
2073 
2074  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2075  double a = v * t_w;
2076 
2077  int rr = 0;
2078  for (; rr != row_nb_dofs / 3; ++rr) {
2079 
2080  auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
2081  auto t_m = get_ftensor1(m, 3 * rr, 0);
2082  for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
2083  const double b = a * t_row_base_fun * t_col_base_fun;
2084  t_m(i) += b;
2085  ++t_m;
2086  ++t_col_base_fun;
2087  }
2088 
2089  ++t_row_base_fun;
2090  }
2091 
2092  for (; rr != row_nb_base_functions; ++rr)
2093  ++t_row_base_fun;
2094 
2095  ++t_w;
2096  }
2097 
2099 }
2100 
2101 MoFEMErrorCode OpSpatialSchurEnd::assemble(int row_side, EntityType row_type,
2102  EntData &data) {
2104  if (row_type != MBTET)
2106 
2107  if (auto ep_fe_vol_ptr =
2108  dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
2109  getFEMethod())) {
2110 
2111  AO aoSuu = PETSC_NULL;
2112  AO aoSBB = PETSC_NULL;
2113  AO aoSOO = PETSC_NULL;
2114  AO aoSww = PETSC_NULL;
2115 
2116  Mat Suu = PETSC_NULL;
2117  Mat SBB = PETSC_NULL;
2118  Mat SOO = PETSC_NULL;
2119  Mat Sww = PETSC_NULL;
2120 
2121  auto set_data = [&](auto fe) {
2122  aoSuu = fe->aoSuu;
2123  aoSBB = fe->aoSBubble;
2124  aoSOO = fe->aoSOmega;
2125  aoSww = fe->aoSw;
2126 
2127  Suu = fe->Suu;
2128  SBB = fe->SBubble;
2129  SOO = fe->SOmega;
2130  Sww = fe->Sw;
2131  };
2132  set_data(ep_fe_vol_ptr);
2133 
2134  if (Suu) {
2135 
2136  auto find_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2137  auto &row_name, auto &col_name, auto row_side,
2138  auto col_side, auto row_type, auto col_type) {
2139  return add_bmc.get<0>().find(boost::make_tuple(
2140  row_name, col_name, row_type, col_type, row_side, col_side));
2141  };
2142 
2143  auto set_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2144  auto &row_name, auto &col_name, auto row_side,
2145  auto col_side, auto row_type, auto col_type,
2146  const auto &m, const auto &row_ind,
2147  const auto &col_ind) {
2148  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2149  row_type, col_type);
2150  if (it != add_bmc.get<0>().end()) {
2151  it->setMat(m);
2152  it->setInd(row_ind, col_ind);
2153  it->setSetAtElement();
2154  return it;
2155  } else {
2156  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2157  row_name, col_name, row_type, col_type, row_side, col_side, m,
2158  row_ind, col_ind));
2159  return p_it.first;
2160  }
2161  };
2162 
2163  auto add_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2164  auto &row_name, auto &col_name, auto row_side,
2165  auto col_side, auto row_type, auto col_type,
2166  const auto &m, const auto &row_ind,
2167  const auto &col_ind) {
2168  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2169  row_type, col_type);
2170  if (it != add_bmc.get<0>().end()) {
2171  it->addMat(m);
2172  it->setInd(row_ind, col_ind);
2173  it->setSetAtElement();
2174  return it;
2175  } else {
2176  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2177  row_name, col_name, row_type, col_type, row_side, col_side, m,
2178  row_ind, col_ind));
2179  return p_it.first;
2180  }
2181  };
2182 
2183  auto assemble_block = [&](auto &bit, Mat S) {
2185  const VectorInt &rind = bit.rowInd;
2186  const VectorInt &cind = bit.colInd;
2187  const MatrixDouble &m = bit.M;
2188  CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
2189  &*cind.begin(), &*m.data().begin(), ADD_VALUES);
2190 
2192  };
2193 
2194  auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
2196  const int nb = m.size1();
2197 
2198  inv.resize(nb, nb, false);
2199  inv.clear();
2200  for (int cc = 0; cc != nb; ++cc)
2201  inv(cc, cc) = -1;
2202 
2203  iPIV.resize(nb, false);
2204  lapackWork.resize(nb * nb, false);
2205  const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
2206  &*iPIV.begin(), &*inv.data().begin(), nb,
2207  &*lapackWork.begin(), nb * nb);
2208  if (info != 0)
2209  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2210  "Can not invert matrix info = %d", info);
2211 
2213  };
2214 
2215  auto invert_symm_schur = [&](DataAtIntegrationPts::BlockMatContainor &bmc,
2216  std::string field, auto &inv) {
2218 
2219  auto bit =
2220  bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
2221  if (bit != bmc.get<1>().end()) {
2222 
2223  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2224  CHKERR invert_symm_mat(m, inv);
2225 
2226  } else
2227  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2228  "%s matrix not found", field.c_str());
2229 
2231  };
2232 
2233  auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
2235 
2236  const int nb = m.size1();
2237 
2238  MatrixDouble trans_m = trans(m);
2239  MatrixDouble trans_inv;
2240  trans_inv.resize(nb, nb, false);
2241  trans_inv.clear();
2242  for (int c = 0; c != nb; ++c)
2243  trans_inv(c, c) = -1;
2244 
2245  iPIV.resize(nb, false);
2246  const auto info =
2247  lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb, &*iPIV.begin(),
2248  &*trans_inv.data().begin(), nb);
2249  if (info != 0)
2250  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2251  "Can not invert matrix info = %d", info);
2252 
2253  inv.resize(nb, nb, false);
2254  noalias(inv) = trans(trans_inv);
2255 
2257  };
2258 
2259  auto invert_nonsymm_schur =
2260  [&](DataAtIntegrationPts::BlockMatContainor &bmc, std::string field,
2261  auto &inv, const bool debug = false) {
2263 
2264  auto bit = bmc.get<1>().find(
2265  boost::make_tuple(field, field, MBTET, MBTET));
2266  if (bit != bmc.get<1>().end()) {
2267 
2268  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2269  CHKERR invert_nonsymm_mat(m, inv);
2270 
2271  if (debug) {
2272  std::cerr << prod(m, inv) << endl;
2273  std::cerr << endl;
2274  }
2275 
2276  } else
2277  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2278  "%s matrix not found", field.c_str());
2279 
2281  };
2282 
2283  auto create_block_schur =
2286  std::string field, AO ao, MatrixDouble &inv) {
2288 
2289  for (auto &bit : add_bmc) {
2290  bit.unSetAtElement();
2291  bit.clearMat();
2292  }
2293 
2294  for (auto &bit : bmc) {
2295  if (bit.setAtElement && bit.rowField != field &&
2296  bit.colField != field) {
2297  VectorInt rind = bit.rowInd;
2298  VectorInt cind = bit.colInd;
2299  const MatrixDouble &m = bit.M;
2300  if (ao) {
2301  CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
2302  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2303  }
2304  auto it = set_block(add_bmc, bit.rowField, bit.colField,
2305  bit.rowSide, bit.colSide, bit.rowType,
2306  bit.colType, m, rind, cind);
2307  }
2308  }
2309 
2310  for (auto &bit_col : bmc) {
2311  if (bit_col.setAtElement && bit_col.rowField == field &&
2312  bit_col.colField != field) {
2313  const MatrixDouble &cm = bit_col.M;
2314  VectorInt cind = bit_col.colInd;
2315  invMat.resize(inv.size1(), cm.size2(), false);
2316  noalias(invMat) = prod(inv, cm);
2317  if (ao)
2318  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2319  for (auto &bit_row : bmc) {
2320  if (bit_row.setAtElement && bit_row.rowField != field &&
2321  bit_row.colField == field) {
2322  const MatrixDouble &rm = bit_row.M;
2323  VectorInt rind = bit_row.rowInd;
2324  K.resize(rind.size(), cind.size(), false);
2325  noalias(K) = prod(rm, invMat);
2326  if (ao)
2327  CHKERR AOApplicationToPetsc(ao, rind.size(),
2328  &*rind.begin());
2329  auto it = add_block(add_bmc, bit_row.rowField,
2330  bit_col.colField, bit_row.rowSide,
2331  bit_col.colSide, bit_row.rowType,
2332  bit_col.colType, K, rind, cind);
2333  }
2334  }
2335  }
2336  }
2337 
2339  };
2340 
2341  auto assemble_schur =
2342  [&](DataAtIntegrationPts::BlockMatContainor &add_bmc, Mat S,
2343  bool debug = false) {
2345  for (auto &bit : add_bmc) {
2346  if (bit.setAtElement)
2347  CHKERR assemble_block(bit, S);
2348  }
2349  if (debug) {
2350  for (auto &bit : add_bmc) {
2351  if (bit.setAtElement) {
2352  std::cerr << "assemble: " << bit.rowField << " "
2353  << bit.colField << endl;
2354  std::cerr << bit.M << endl;
2355  }
2356  }
2357  std::cerr << std::endl;
2358  }
2360  };
2361 
2362  auto precondition_schur =
2365  const std::string field, const MatrixDouble &diag_mat,
2366  const double eps) {
2368 
2369  for (auto &bit : add_bmc) {
2370  bit.unSetAtElement();
2371  bit.clearMat();
2372  }
2373 
2374  for (auto &bit : bmc) {
2375  if (bit.setAtElement) {
2376  if (bit.rowField != field || bit.colField != field)
2377  auto it =
2378  set_block(add_bmc, bit.rowField, bit.colField,
2379  bit.rowSide, bit.colSide, bit.rowType,
2380  bit.colType, bit.M, bit.rowInd, bit.colInd);
2381  }
2382  }
2383 
2384  auto bit = bmc.get<1>().find(
2385  boost::make_tuple(field, field, MBTET, MBTET));
2386  if (bit->setAtElement && bit != bmc.get<1>().end()) {
2387  auto it =
2388  set_block(add_bmc, bit->rowField, bit->colField, bit->rowSide,
2389  bit->colSide, bit->rowType, bit->colType, bit->M,
2390  bit->rowInd, bit->colInd);
2391  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2392  m += eps * diag_mat;
2393  } else {
2394  auto row_it = bmc.get<3>().lower_bound(field);
2395  for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
2396  if (row_it->setAtElement) {
2397  auto it = set_block(add_bmc, field, field, 0, 0, MBTET, MBTET,
2398  diag_mat, row_it->rowInd, row_it->rowInd);
2399  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2400  m *= eps;
2401  break;
2402  }
2403  }
2404  if (row_it == bmc.get<3>().end())
2405  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2406  "row field not found %s", field.c_str());
2407  }
2408 
2410  };
2411 
2412  CHKERR invert_symm_schur(dataAtPts->blockMatContainor, "u",
2413  invBlockMat["uu"]);
2414  CHKERR create_block_schur(dataAtPts->blockMatContainor, blockMat["uu"],
2415  "u", aoSuu, invBlockMat["uu"]);
2416  CHKERR assemble_schur(blockMat["uu"], Suu);
2417 
2418  if (SBB) {
2419  CHKERR invert_symm_schur(blockMat["uu"], "BUBBLE", invBlockMat["BB"]);
2420  CHKERR create_block_schur(blockMat["uu"], blockMat["BB"], "BUBBLE",
2421  aoSBB, invBlockMat["BB"]);
2422  CHKERR precondition_schur(blockMat["BB"], blockMat["precBBOO"], "omega",
2423  *dataAtPts->ooMatPtr, eps);
2424  CHKERR assemble_schur(blockMat["precBBOO"], SBB);
2425 
2426  if (SOO) {
2427  CHKERR invert_nonsymm_schur(blockMat["precBBOO"], "omega",
2428  invBlockMat["OO"]);
2429  CHKERR create_block_schur(blockMat["precBBOO"], blockMat["OO"],
2430  "omega", aoSOO, invBlockMat["OO"]);
2431  CHKERR precondition_schur(blockMat["OO"], blockMat["precOOww"], "w",
2432  *dataAtPts->wwMatPtr, eps);
2433  CHKERR assemble_schur(blockMat["precOOww"], SOO);
2434 
2435  if (Sww) {
2436  CHKERR invert_symm_schur(blockMat["precOOww"], "w",
2437  invBlockMat["ww"]);
2438  CHKERR create_block_schur(blockMat["precOOww"], blockMat["ww"], "w",
2439  aoSww, invBlockMat["ww"]);
2440  CHKERR assemble_schur(blockMat["ww"], Sww);
2441  }
2442  }
2443  }
2444  }
2445  }
2446 
2448 } // namespace EshelbianPlasticity
2449 
2450 MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
2451  EntData &data) {
2453  if (type != MBVERTEX)
2455 
2456  auto create_tag = [this](const std::string tag_name, const int size) {
2457  double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2458  Tag th;
2459  CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
2460  th, MB_TAG_CREAT | MB_TAG_SPARSE,
2461  def_VAL);
2462  return th;
2463  };
2464 
2465  Tag th_w = create_tag("SpatialDisplacement", 3);
2466  Tag th_omega = create_tag("Omega", 3);
2467  Tag th_approxP = create_tag("Piola2Stress", 9);
2468  Tag th_sigma = create_tag("CauchyStress", 9);
2469  Tag th_u = create_tag("SpatialStreach", 9);
2470  Tag th_h = create_tag("h", 9);
2471  Tag th_x = create_tag("x", 3);
2472  Tag th_X = create_tag("X", 3);
2473  Tag th_detF = create_tag("detF", 1);
2474  Tag th_angular_momentum = create_tag("AngularMomentum", 3);
2475 
2476  int nb_gauss_pts = data.getN().size1();
2477  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
2478  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2479  "Nb. of integration points is not equal to number points on "
2480  "post-processing mesh");
2481  }
2482 
2483  auto t_w = getFTensor1FromMat<3>(*(dataAtPts->wAtPts));
2484  auto t_x = getFTensor1FromMat<3>(*(dataAtPts->xAtPts));
2485  auto t_omega = getFTensor1FromMat<3>(*(dataAtPts->rotAxisAtPts));
2486  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
2487  auto t_u = getFTensor2SymmetricFromMat<3>(*(dataAtPts->streachTensorAtPts));
2488  auto t_R = getFTensor2FromMat<3, 3>(*(dataAtPts->rotMatAtPts));
2489  auto t_approx_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
2490  auto t_coords = getFTensor1CoordsAtGaussPts();
2491 
2496 
2497  // vectors
2498  VectorDouble3 v(3);
2499  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
2500  auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
2502  t_v(i) = t_d(i);
2503  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2504  &*v.data().begin());
2506  };
2507 
2508  MatrixDouble3by3 m(3, 3);
2510  &m(0, 0), &m(0, 1), &m(0, 2),
2511 
2512  &m(1, 0), &m(1, 1), &m(1, 2),
2513 
2514  &m(2, 0), &m(2, 1), &m(2, 2));
2515 
2516  auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
2518  t_m(i, j) = t_d(i, j);
2519  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
2520  &*m.data().begin());
2522  };
2523 
2524  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2525 
2526  // vetors
2527  CHKERR save_vec_tag(th_w, t_w, gg);
2528  CHKERR save_vec_tag(th_x, t_x, gg);
2529  CHKERR save_vec_tag(th_X, t_coords, gg);
2530  CHKERR save_vec_tag(th_omega, t_omega, gg);
2531 
2532  // tensors
2533  CHKERR save_mat_tag(th_h, t_h, gg);
2535  for (int ii = 0; ii != 3; ++ii)
2536  for (int jj = 0; jj != 3; ++jj)
2537  t_u_tmp(ii, jj) = t_u(ii, jj);
2538  CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2539  CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2540 
2541  const double jac = dEterminant(t_h);
2543  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
2544  CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2545  CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
2546 
2548  t_PhT(i, k) = t_approx_P(i, j) * t_R(k, j);
2549  FTensor::Tensor1<double, 3> t_leviPRT;
2550  t_leviPRT(k) = levi_civita(i, l, k) * t_PhT(i, l);
2551 
2552  CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
2553  &t_leviPRT(0));
2554 
2555  ++t_w;
2556  ++t_x;
2557  ++t_h;
2558  ++t_u;
2559  ++t_omega;
2560  ++t_R;
2561  ++t_approx_P;
2562  ++t_coords;
2563  }
2564 
2566 }
2567 
2569 OpCalculateStrainEnergy::doWork(int side, EntityType type,
2572  int nb_integration_pts = data.getN().size1();
2573  auto v = getVolume();
2574  auto t_w = getFTensor0IntegrationWeight();
2575  auto t_P = getFTensor2FromMat<3, 3>(*(dataAtPts->approxPAtPts));
2576  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hAtPts));
2577 
2580 
2581  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2582  const double a = t_w * v;
2583  (*energy) += a * t_P(i, J) * t_h(i, J);
2584  ++t_w;
2585  ++t_P;
2586  ++t_h;
2587  }
2588 
2590 }
2591 
2592 } // namespace EshelbianPlasticity
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:142
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:415
multi_index_container< BlockMatData, indexed_by< ordered_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType >, member< BlockMatData, int, &BlockMatData::rowSide >, member< BlockMatData, int, &BlockMatData::colSide > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField > > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::rowField > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::colField > > > > BlockMatContainor
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
FTensor::Tensor3< double, 3, 3, 3 > getDiffRotationFormVector(FTensor::Tensor1< T, 3 > &t_omega)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:396
const VectorInt & getIndices() const
Get global indices of dofs on entity.
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
static __CLPK_integer lapack_dposv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:223
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:385
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
static const bool debug
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:95
#define CHKERR
Inline error check.
Definition: definitions.h:596
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:71
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Eshelbian plasticity interface.
Data on single entity (This is passed as argument to DataOperator::doWork)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
static const double eps
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
static __CLPK_integer lapack_dsysv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:232
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor2< double, 3, 3 > getRotationFormVector(FTensor::Tensor1< T, 3 > &t_omega)
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:188
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
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