v0.14.0
EshelbianADOL-C.cpp
Go to the documentation of this file.
1 /**
2  * \file EshelbianPlasticity.cpp
3  * \brief Implementation of automatic differentiation
4  */
5 
6 #include <MoFEM.hpp>
7 using namespace MoFEM;
8 
9 #include <BasicFiniteElements.hpp>
10 
11 #include <EshelbianPlasticity.hpp>
12 #include <boost/math/constants/constants.hpp>
13 
14 constexpr double third = boost::math::constants::third<double>();
15 
16 namespace EshelbianPlasticity {
17 
18 struct OpHMHH : public OpJacobian {
19  OpHMHH(const int tag, const bool eval_rhs, const bool eval_lhs,
20  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
21  boost::shared_ptr<PhysicalEquations> &physics_ptr)
22  : OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr) {}
23 
24  MoFEMErrorCode evaluateRhs(EntData &data);
25  MoFEMErrorCode evaluateLhs(EntData &data);
26 };
27 
28 MoFEMErrorCode OpHMHH::evaluateRhs(EntData &ent_data) {
29  MoFEMFunctionBegin;
30  FTensor::Index<'i', 3> i;
31  FTensor::Index<'j', 3> j;
32  FTensor::Index<'k', 3> k;
33  FTensor::Index<'l', 3> l;
34  FTensor::Index<'m', 3> m;
35  FTensor::Index<'n', 3> n;
36 
37  const auto nb_integration_pts = getGaussPts().size2();
38  auto iu = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
39  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
40 
41  auto create_data_vec = [nb_integration_pts](auto &v) {
42  v.resize(nb_integration_pts, false);
43  };
44  auto create_data_mat = [nb_integration_pts](auto &m) {
45  m.resize(9, nb_integration_pts, false);
46  };
47 
48  create_data_mat(dataAtPts->PAtPts);
49 
50  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
51  auto r_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
52  for (unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
53 
54  FTensor::Tensor2<double, 3, 3> t_h1;
55 
56  switch (EshelbianCore::gradApperoximator) {
57  case LARGE_ROT:
58  case MODERATE_ROT:
59  t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
60  physicsPtr->get_h()(i, j) = iu(i, k) * t_h1(k, j);
61  break;
62  case SMALL_ROT:
63  physicsPtr->get_h()(i, j) = iu(i, j);
64  break;
65  };
66 
67  int r = ::function(tAg, physicsPtr->dependentVariablesPiola.size(),
68  physicsPtr->activeVariables.size(),
69  &physicsPtr->activeVariables[0],
70  &physicsPtr->dependentVariablesPiola[0]);
71  if (r < 0) { // function is locally analytic
72  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
73  "ADOL-C function evaluation with error r = %d", r);
74  }
75 
76  switch (EshelbianCore::gradApperoximator) {
77  case LARGE_ROT:
78  case MODERATE_ROT: {
79  FTensor::Tensor4<double, 3, 3, 3, 3> t_h1_du;
80  t_h1_du(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * t_h1(k, j));
81  r_P(k, l) = physicsPtr->get_P()(i, j) * t_h1_du(i, j, k, l);
82  }; break;
83  case SMALL_ROT:
84  r_P(i, j) = physicsPtr->get_P()(i, j);
85  break;
86  };
87 
88  ++iu;
89  ++t_grad_h1;
90  ++r_P;
91  }
92  MoFEMFunctionReturn(0);
93 }
94 
95 MoFEMErrorCode OpHMHH::evaluateLhs(EntData &ent_data) {
96  MoFEMFunctionBegin;
97  FTensor::Index<'i', 3> i;
98  FTensor::Index<'j', 3> j;
99  FTensor::Index<'k', 3> k;
100  FTensor::Index<'l', 3> l;
101  FTensor::Index<'m', 3> m;
102  FTensor::Index<'n', 3> n;
103  FTensor::Index<'o', 3> o;
104  FTensor::Index<'p', 3> p;
105 
106  const int number_of_active_variables = physicsPtr->activeVariables.size();
107  const int number_of_dependent_variables =
108  physicsPtr->dependentVariablesPiola.size();
109  std::vector<double *> jac_ptr(number_of_dependent_variables);
110  for (unsigned int n = 0; n != number_of_dependent_variables; ++n) {
111  jac_ptr[n] =
112  &(physicsPtr
113  ->dependentVariablesPiolaDirevatives[n *
114  number_of_active_variables]);
115  }
116 
117  const auto nb_integration_pts = getGaussPts().size2();
118 
119  auto create_data_mat = [nb_integration_pts](auto &m) {
120  m.resize(9, nb_integration_pts, false);
121  };
122 
123  dataAtPts->P_du.resize(81, nb_integration_pts, false);
124 
125  auto iu = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
126  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
127  auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
128 
129  FTensor::Number<0> N0;
130  FTensor::Number<1> N1;
131  FTensor::Number<2> N2;
132 
133  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
134 
135  for (unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
136 
137  FTensor::Tensor2<double, 3, 3> t_h1;
138 
139  switch (EshelbianCore::gradApperoximator) {
140  case LARGE_ROT:
141  case MODERATE_ROT:
142  t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
143  physicsPtr->get_h()(i, j) = iu(i, k) * t_h1(k, j);
144  break;
145  case SMALL_ROT:
146  physicsPtr->get_h()(i, j) = iu(i, j);
147  break;
148  };
149 
150  // play recorder for jacobians
151  int r = ::jacobian(tAg, number_of_dependent_variables,
152  number_of_active_variables,
153  &physicsPtr->activeVariables[0], &jac_ptr[0]);
154  if (r < 0) {
155  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
156  "ADOL-C function evaluation with error");
157  }
158 
159  FTensor::Tensor4<double, 3, 3, 3, 3> t_P_dh_tmp;
160  t_P_dh_tmp(i, j, N0, k) = physicsPtr->get_P_dh0()(i, j, k);
161  t_P_dh_tmp(i, j, N1, k) = physicsPtr->get_P_dh1()(i, j, k);
162  t_P_dh_tmp(i, j, N2, k) = physicsPtr->get_P_dh2()(i, j, k);
163 
164  switch (EshelbianCore::gradApperoximator) {
165  case LARGE_ROT:
166  case MODERATE_ROT: {
167  FTensor::Tensor4<double, 3, 3, 3, 3> t_h1_du;
168  t_h1_du(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * t_h1(k, j));
169  r_P_du(k, l, m, n) =
170  (t_P_dh_tmp(i, j, o, p) * t_h1_du(o, p, m, n)) * t_h1_du(i, j, k, l);
171  } break;
172  case SMALL_ROT:
173  r_P_du(i, j, m, n) = t_P_dh_tmp(i, j, m, n);
174  break;
175  };
176 
177  ++iu;
178  ++t_grad_h1;
179  ++r_P_du;
180  }
181  MoFEMFunctionReturn(0);
182 }
183 
184 struct HMHStVenantKirchhoff : public PhysicalEquations {
185 
186  static constexpr int numberOfActiveVariables = 9;
187  static constexpr int numberOfDependentVariables = 9;
188 
189  HMHStVenantKirchhoff(const double lambda, const double mu)
190  : PhysicalEquations(numberOfActiveVariables, numberOfDependentVariables),
191  lambda(lambda), mu(mu) {}
192 
193  MoFEMErrorCode getOptions() {
194  MoFEMFunctionBegin;
195 
196  double E = 1;
197  double nu = 0;
198 
199  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "stvenant_", "", "none");
200 
201  CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
202  PETSC_NULL);
203  CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
204  PETSC_NULL);
205 
206  ierr = PetscOptionsEnd();
207  CHKERRG(ierr);
208 
209  MOFEM_LOG("EP", Sev::inform) << "St Venant Kirchhoff model parameters: "
210  << "E = " << E << ", nu = " << nu;
211 
212  lambda = LAMBDA(E, nu);
213  mu = MU(E, nu);
214 
215  MoFEMFunctionReturn(0);
216  }
217 
218  virtual OpJacobian *
219  returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
220  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
221  boost::shared_ptr<PhysicalEquations> &physics_ptr) {
222  return (new OpHMHH(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
223  }
224 
225  double lambda;
226  double mu;
227  double sigmaY;
228 
229  ATensor2 th;
230  ATensor2 tH;
231  ATensor2 tP;
232  ATensor2 tSigma;
233  adouble phi;
234  adouble s2;
235  adouble f;
236  adouble energy;
237 
238  adouble detH;
239  adouble detF;
240  adouble trE;
241  ATensor2 tInvH;
242  ATensor2 tF;
243  ATensor2 tInvF;
244  ATensor2 tC;
245  ATensor2 tE;
246  ATensor2 tS;
247 
248  ATensor2 tPulledP;
249 
250  MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr) {
251  MoFEMFunctionBegin;
252  FTensor::Index<'i', 3> i;
253  FTensor::Index<'j', 3> j;
254  FTensor::Index<'I', 3> I;
255  FTensor::Index<'J', 3> J;
256  FTensor::Number<0> N0;
257  FTensor::Number<1> N1;
258  FTensor::Number<2> N2;
259 
260  CHKERR getOptions();
261 
262  auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
263 
264  auto ih = get_h();
265  // auto iH = get_H();
266  if (t_h_ptr)
267  ih(i, j) = (*t_h_ptr)(i, j);
268  else {
269  ih(i, j) = t_kd(i, j);
270  }
271 
272  auto r_P = get_P();
273 
274  enableMinMaxUsingAbs();
275  trace_on(tape);
276 
277  // Set active variables to ADOL-C
278  th(i, j) <<= get_h()(i, j);
279 
280  tF(i, I) = th(i, I);
281  CHKERR determinantTensor3by3(tF, detF);
282  CHKERR invertTensor3by3(tF, detF, tInvF);
283 
284  // Deformation and strain
285  tC(I, J) = tF(i, I) * tF(i, J);
286  tE(I, J) = tC(I, J);
287  tE(N0, N0) -= 1;
288  tE(N1, N1) -= 1;
289  tE(N2, N2) -= 1;
290  tE(I, J) *= 0.5;
291 
292  // Energy
293  trE = tE(I, I);
294  energy = 0.5 * lambda * trE * trE + mu * (tE(I, J) * tE(I, J));
295 
296  // Stress Piola II
297  tS(I, J) = tE(I, J);
298  tS(I, J) *= 2 * mu;
299  tS(N0, N0) += lambda * trE;
300  tS(N1, N1) += lambda * trE;
301  tS(N2, N2) += lambda * trE;
302  // Stress Piola I
303  tP(i, J) = tF(i, I) * tS(I, J);
304 
305  // Set dependent variables to ADOL-C
306  tP(i, j) >>= r_P(i, j);
307 
308  trace_off();
309 
310  MoFEMFunctionReturn(0);
311  }
312 };
313 
314 struct HMHPMooneyRivlinWriggersEq63 : public PhysicalEquations {
315 
316  static constexpr int numberOfActiveVariables = 9;
317  static constexpr int numberOfDependentVariables = 9;
318 
319  HMHPMooneyRivlinWriggersEq63(const double alpha, const double beta,
320  const double lambda)
321  : PhysicalEquations(numberOfActiveVariables, numberOfDependentVariables),
322  alpha(alpha), beta(beta), lambda(lambda), epsilon(0) {}
323 
324  virtual OpJacobian *
325  returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
326  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
327  boost::shared_ptr<PhysicalEquations> &physics_ptr) {
328  return (new OpHMHH(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
329  }
330 
331  MoFEMErrorCode getOptions() {
332  MoFEMFunctionBegin;
333  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "mooneyrivlin_", "", "none");
334 
335  alpha = 1;
336  CHKERR PetscOptionsScalar("-alpha", "Alpha", "", alpha, &alpha, PETSC_NULL);
337  beta = 1;
338  CHKERR PetscOptionsScalar("-beta", "Beta", "", beta, &beta, PETSC_NULL);
339 
340  lambda = 1;
341  CHKERR PetscOptionsScalar("-lambda", "Lambda", "", lambda, &lambda,
342  PETSC_NULL);
343 
344  epsilon = 0;
345  CHKERR PetscOptionsScalar("-epsilon", "Epsilon", "", epsilon, &epsilon,
346  PETSC_NULL);
347 
348  ierr = PetscOptionsEnd();
349  CHKERRG(ierr);
350  MoFEMFunctionReturn(0);
351  }
352 
353  double alpha;
354  double beta;
355  double lambda;
356  double epsilon;
357 
358  ATensor2 th;
359  ATensor2 tH;
360  ATensor2 tF;
361 
362  adouble detH;
363  adouble detF;
364  ATensor2 tInvH;
365  ATensor2 tInvF;
366 
367  ATensor2 tP;
368  ATensor2 tSigma;
369 
370  ATensor2 tCof;
371  // ATensor3 tCofT1;
372  // ATensor3 tCofT2;
373  ATensor2 tBF;
374  ATensor2 tBCof;
375  adouble tBj;
376 
377  adouble energy;
378  adouble phi;
379  adouble A;
380  adouble B;
381 
382  ATensor2 tPulledP;
383 
384  MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr) {
385  MoFEMFunctionBegin;
386 
387  FTensor::Index<'i', 3> i;
388  FTensor::Index<'j', 3> j;
389  FTensor::Index<'k', 3> k;
390  FTensor::Index<'I', 3> I;
391  FTensor::Index<'J', 3> J;
392  FTensor::Index<'K', 3> K;
393  FTensor::Number<0> N0;
394  FTensor::Number<1> N1;
395  FTensor::Number<2> N2;
396 
397  CHKERR getOptions();
398 
399  auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
400 
401  auto ih = get_h();
402  // auto iH = get_H();
403  if (t_h_ptr)
404  ih(i, j) = (*t_h_ptr)(i, j);
405  else {
406  ih(i, j) = t_kd(i, j);
407  }
408 
409  auto r_P = get_P();
410 
411  enableMinMaxUsingAbs();
412  trace_on(tape);
413 
414  // Set active variables to ADOL-C
415  th(i, j) <<= get_h()(i, j);
416 
417  tF(i, I) = th(i, I);
418  CHKERR determinantTensor3by3(tF, detF);
419  CHKERR invertTensor3by3(tF, detF, tInvF);
420 
421  tCof(i, I) = detF * tInvF(I, i);
422 
423  A = tF(k, K) * tF(k, K);
424  B = tCof(k, K) * tCof(k, K);
425 
426  tBF(i, I) = 4 * alpha * (A * tF(i, I));
427  tBCof(i, I) = 4 * beta * (B * tCof(i, I));
428  tBj = (-12 * alpha - 24 * beta) / detF +
429  0.5 * (lambda / epsilon) *
430  (pow(detF, epsilon - 1) - pow(detF, -epsilon - 1));
431 
432  tP(i, I) = tBF(i, I);
433  tP(i, I) += (levi_civita(i, j, k) * tBCof(j, J)) *
434  (levi_civita(I, J, K) * tF(k, K));
435  tP(i, I) += tCof(i, I) * tBj;
436 
437  // Set dependent variables to ADOL-C
438  tP(i, j) >>= r_P(i, j);
439 
440  trace_off();
441 
442  MoFEMFunctionReturn(0);
443  }
444 };
445 
446 MoFEMErrorCode EshelbianCore::addMaterial_HMHHStVenantKirchhoff(
447  const int tape, const double lambda, const double mu,
448  const double sigma_y) {
449  MoFEMFunctionBegin;
450  physicalEquations = boost::make_shared<HMHStVenantKirchhoff>(lambda, mu);
451  CHKERR physicalEquations->recordTape(tape, nullptr);
452  MoFEMFunctionReturn(0);
453 }
454 
455 MoFEMErrorCode EshelbianCore::addMaterial_HMHMooneyRivlin(
456  const int tape, const double alpha, const double beta, const double lambda,
457  const double sigma_y) {
458  MoFEMFunctionBegin;
459  physicalEquations =
460  boost::make_shared<HMHPMooneyRivlinWriggersEq63>(alpha, beta, lambda);
461  CHKERR physicalEquations->recordTape(tape, nullptr);
462  MoFEMFunctionReturn(0);
463 }
464 
465 struct HMHHencky : public PhysicalEquations {
466 
467  HMHHencky(MoFEM::Interface &m_field, const double E, const double nu)
468  : PhysicalEquations(0, 0), mField(m_field), E(E), nu(nu) {}
469 
470  MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h) { return 0; }
471 
472  struct OpHenckyJacobian : public OpJacobian {
473  OpHenckyJacobian(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
474  boost::shared_ptr<HMHHencky> hencky_ptr)
475  : OpJacobian(), dataAtGaussPts(data_ptr), henckyPtr(hencky_ptr) {
476  CHK_THROW_MESSAGE(henckyPtr->getOptions(dataAtGaussPts),
477  "getOptions failed");
478  CHK_THROW_MESSAGE(henckyPtr->extractBlockData(Sev::inform),
479  "Can not get data from block");
480  }
481 
482  MoFEMErrorCode doWork(int side, EntityType type,
483  EntitiesFieldData::EntData &data) {
484  MoFEMFunctionBegin;
485 
486  for (auto &b : henckyPtr->blockData) {
487 
488  if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
489  dataAtGaussPts->mu = b.bulkModulusK;
490  dataAtGaussPts->lambda = b.bulkModulusK;
491  CHKERR henckyPtr->getMatDPtr(dataAtGaussPts->getMatDPtr(),
492  dataAtGaussPts->getMatAxiatorDPtr(),
493  dataAtGaussPts->getMatDeviatorDPtr(),
494  b.bulkModulusK, b.shearModulusG);
495  MoFEMFunctionReturnHot(0);
496  }
497  }
498 
499  const auto E = henckyPtr->E;
500  const auto nu = henckyPtr->nu;
501 
502  double bulk_modulus_K = E / (3 * (1 - 2 * nu));
503  double shear_modulus_G = E / (2 * (1 + nu));
504 
505  dataAtGaussPts->mu = shear_modulus_G;
506  dataAtGaussPts->lambda = bulk_modulus_K;
507 
508  CHKERR henckyPtr->getMatDPtr(dataAtGaussPts->getMatDPtr(),
509  dataAtGaussPts->getMatAxiatorDPtr(),
510  dataAtGaussPts->getMatDeviatorDPtr(),
511  bulk_modulus_K, shear_modulus_G);
512 
513  MoFEMFunctionReturn(0);
514  }
515 
516  MoFEMErrorCode evaluateRhs(EntData &data) { return 0; }
517  MoFEMErrorCode evaluateLhs(EntData &data) { return 0; }
518 
519  private:
520  boost::shared_ptr<DataAtIntegrationPts> dataAtGaussPts;
521  boost::shared_ptr<HMHHencky> henckyPtr;
522  };
523 
524  virtual OpJacobian *
525  returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
526  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
527  boost::shared_ptr<PhysicalEquations> &physics_ptr) {
528  return (new OpHenckyJacobian(
529  data_ptr, boost::dynamic_pointer_cast<HMHHencky>(physics_ptr)));
530  }
531 
532  MoFEMErrorCode getOptions(boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
533  MoFEMFunctionBegin;
534  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
535 
536  CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
537  PETSC_NULL);
538  CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
539  PETSC_NULL);
540 
541  ierr = PetscOptionsEnd();
542  CHKERRG(ierr);
543 
544  MoFEMFunctionReturn(0);
545  }
546 
547  MoFEMErrorCode extractBlockData(Sev sev) {
548  return extractBlockData(
549 
550  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
551 
552  (boost::format("%s(.*)") % "MAT_ELASTIC").str()
553 
554  )),
555 
556  sev);
557  }
558 
559  MoFEMErrorCode
560  extractBlockData(std::vector<const CubitMeshSets *> meshset_vec_ptr,
561  Sev sev) {
562  MoFEMFunctionBegin;
563 
564  for (auto m : meshset_vec_ptr) {
565  MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
566  std::vector<double> block_data;
567  CHKERR m->getAttributes(block_data);
568  if (block_data.size() != 2) {
569  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
570  "Expected that block has two attribute");
571  }
572  auto get_block_ents = [&]() {
573  Range ents;
574  CHKERR mField.get_moab().get_entities_by_handle(m->meshset, ents, true);
575  return ents;
576  };
577 
578  double young_modulus = block_data[0];
579  double poisson_ratio = block_data[1];
580  double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
581  double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
582 
583  MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
584  << "E = " << young_modulus << " nu = " << poisson_ratio;
585 
586  blockData.push_back({bulk_modulus_K, shear_modulus_G, get_block_ents()});
587  }
588  MOFEM_LOG_CHANNEL("WORLD");
589  MoFEMFunctionReturn(0);
590  }
591 
592  MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
593  boost::shared_ptr<MatrixDouble> mat_axiator_D_ptr,
594  boost::shared_ptr<MatrixDouble> mat_deviator_D_ptr,
595  double bulk_modulus_K, double shear_modulus_G) {
596  MoFEMFunctionBegin;
597  //! [Calculate elasticity tensor]
598  auto set_material_stiffness = [&]() {
599  FTensor::Index<'i', SPACE_DIM> i;
600  FTensor::Index<'j', SPACE_DIM> j;
601  FTensor::Index<'k', SPACE_DIM> k;
602  FTensor::Index<'l', SPACE_DIM> l;
603  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
604  double A = (SPACE_DIM == 2)
605  ? 2 * shear_modulus_G /
606  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
607  : 1;
608  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
609  auto t_axiator_D =
610  getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_axiator_D_ptr);
611  auto t_deviator_D =
612  getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_deviator_D_ptr);
613  t_axiator_D(i, j, k, l) = A *
614  (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
615  t_kd(i, j) * t_kd(k, l);
616  t_deviator_D(i, j, k, l) =
617  2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
618  t_D(i, j, k, l) = t_axiator_D(i, j, k, l) + t_deviator_D(i, j, k, l);
619  };
620  //! [Calculate elasticity tensor]
621  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
622  mat_D_ptr->resize(size_symm * size_symm, 1);
623  mat_axiator_D_ptr->resize(size_symm * size_symm, 1);
624  mat_deviator_D_ptr->resize(size_symm * size_symm, 1);
625  set_material_stiffness();
626  MoFEMFunctionReturn(0);
627  }
628 
629 private:
630  MoFEM::Interface &mField;
631 
632  struct BlockData {
633  double bulkModulusK;
634  double shearModulusG;
635  Range blockEnts;
636  };
637  std::vector<BlockData> blockData;
638 
639  double E;
640  double nu;
641 };
642 
643 MoFEMErrorCode EshelbianCore::addMaterial_Hencky(double E, double nu) {
644  MoFEMFunctionBegin;
645  physicalEquations = boost::make_shared<HMHHencky>(mField, E, nu);
646  MoFEMFunctionReturn(0);
647 }
648 
649 }; // namespace EshelbianPlasticity
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
EshelbianPlasticity::HMHStVenantKirchhoff::sigmaY
double sigmaY
Definition: EshelbianADOL-C.cpp:227
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::energy
adouble energy
Definition: EshelbianADOL-C.cpp:377
EshelbianPlasticity::HMHHencky::OpHenckyJacobian::evaluateRhs
MoFEMErrorCode evaluateRhs(EntData &data)
Definition: EshelbianADOL-C.cpp:516
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
EshelbianPlasticity::HMHHencky::getOptions
MoFEMErrorCode getOptions(boost::shared_ptr< DataAtIntegrationPts > data_ptr)
Definition: EshelbianADOL-C.cpp:532
EshelbianPlasticity::HMHHencky::recordTape
MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h)
Definition: EshelbianADOL-C.cpp:470
EshelbianPlasticity::HMHStVenantKirchhoff::detF
adouble detF
Definition: EshelbianADOL-C.cpp:239
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:172
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tSigma
ATensor2 tSigma
Definition: EshelbianADOL-C.cpp:368
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::alpha
double alpha
Definition: EshelbianADOL-C.cpp:353
EshelbianPlasticity::HMHStVenantKirchhoff::tInvH
ATensor2 tInvH
Definition: EshelbianADOL-C.cpp:241
EshelbianPlasticity::HMHHencky::OpHenckyJacobian
Definition: EshelbianADOL-C.cpp:472
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
E
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tBCof
ATensor2 tBCof
Definition: EshelbianADOL-C.cpp:374
EshelbianPlasticity::HMHStVenantKirchhoff::phi
adouble phi
Definition: EshelbianADOL-C.cpp:233
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
EshelbianPlasticity::HMHStVenantKirchhoff::returnOpJacobian
virtual OpJacobian * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
Definition: EshelbianADOL-C.cpp:219
EshelbianPlasticity::HMHHencky::OpHenckyJacobian::henckyPtr
boost::shared_ptr< HMHHencky > henckyPtr
Definition: EshelbianADOL-C.cpp:521
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tP
ATensor2 tP
Definition: EshelbianADOL-C.cpp:367
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tInvH
ATensor2 tInvH
Definition: EshelbianADOL-C.cpp:364
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:96
EshelbianPlasticity::HMHHencky::blockData
std::vector< BlockData > blockData
Definition: EshelbianADOL-C.cpp:637
BasicFiniteElements.hpp
EshelbianPlasticity::PhysicalEquations
Definition: EshelbianPlasticity.hpp:181
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
EshelbianPlasticity::HMHHencky::BlockData
Definition: EshelbianADOL-C.cpp:632
EshelbianPlasticity::HMHStVenantKirchhoff::trE
adouble trE
Definition: EshelbianADOL-C.cpp:240
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::recordTape
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
Definition: EshelbianADOL-C.cpp:384
EshelbianPlasticity::HMHStVenantKirchhoff::tInvF
ATensor2 tInvF
Definition: EshelbianADOL-C.cpp:243
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::beta
double beta
Definition: EshelbianADOL-C.cpp:354
EshelbianPlasticity::HMHHencky::nu
double nu
Definition: EshelbianADOL-C.cpp:640
EshelbianPlasticity::HMHHencky::returnOpJacobian
virtual OpJacobian * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
Definition: EshelbianADOL-C.cpp:525
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tInvF
ATensor2 tInvF
Definition: EshelbianADOL-C.cpp:365
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
EshelbianPlasticity::HMHStVenantKirchhoff::detH
adouble detH
Definition: EshelbianADOL-C.cpp:238
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
EshelbianPlasticity::HMHStVenantKirchhoff::tP
ATensor2 tP
Definition: EshelbianADOL-C.cpp:231
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::invertTensor3by3
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:1559
EshelbianPlasticity::HMHHencky::extractBlockData
MoFEMErrorCode extractBlockData(std::vector< const CubitMeshSets * > meshset_vec_ptr, Sev sev)
Definition: EshelbianADOL-C.cpp:560
EshelbianPlasticity::HMHHencky::OpHenckyJacobian::evaluateLhs
MoFEMErrorCode evaluateLhs(EntData &data)
Definition: EshelbianADOL-C.cpp:517
EshelbianPlasticity::HMHHencky::OpHenckyJacobian::OpHenckyJacobian
OpHenckyJacobian(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< HMHHencky > hencky_ptr)
Definition: EshelbianADOL-C.cpp:473
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::HMHPMooneyRivlinWriggersEq63
HMHPMooneyRivlinWriggersEq63(const double alpha, const double beta, const double lambda)
Definition: EshelbianADOL-C.cpp:319
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::th
ATensor2 th
Definition: EshelbianADOL-C.cpp:358
EshelbianPlasticity::HMHStVenantKirchhoff::tPulledP
ATensor2 tPulledP
Definition: EshelbianADOL-C.cpp:248
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
EshelbianPlasticity::HMHHencky::BlockData::bulkModulusK
double bulkModulusK
Definition: EshelbianADOL-C.cpp:633
EshelbianPlasticity::HMHStVenantKirchhoff::tF
ATensor2 tF
Definition: EshelbianADOL-C.cpp:242
EshelbianPlasticity::HMHStVenantKirchhoff::s2
adouble s2
Definition: EshelbianADOL-C.cpp:234
convert.type
type
Definition: convert.py:64
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tH
ATensor2 tH
Definition: EshelbianADOL-C.cpp:359
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tBF
ATensor2 tBF
Definition: EshelbianADOL-C.cpp:373
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::getOptions
MoFEMErrorCode getOptions()
Definition: EshelbianADOL-C.cpp:331
EshelbianPlasticity::HMHHencky
Definition: EshelbianADOL-C.cpp:465
EshelbianPlasticity::HMHHencky::getMatDPtr
MoFEMErrorCode getMatDPtr(boost::shared_ptr< MatrixDouble > mat_D_ptr, boost::shared_ptr< MatrixDouble > mat_axiator_D_ptr, boost::shared_ptr< MatrixDouble > mat_deviator_D_ptr, double bulk_modulus_K, double shear_modulus_G)
Definition: EshelbianADOL-C.cpp:592
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::returnOpJacobian
virtual OpJacobian * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
Definition: EshelbianADOL-C.cpp:325
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:173
EshelbianPlasticity::HMHHencky::E
double E
Definition: EshelbianADOL-C.cpp:639
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::A
adouble A
Definition: EshelbianADOL-C.cpp:379
EshelbianPlasticity::HMHHencky::OpHenckyJacobian::dataAtGaussPts
boost::shared_ptr< DataAtIntegrationPts > dataAtGaussPts
Definition: EshelbianADOL-C.cpp:520
EshelbianPlasticity::OpHMHH
Definition: EshelbianADOL-C.cpp:18
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
EshelbianPlasticity::HMHStVenantKirchhoff::tS
ATensor2 tS
Definition: EshelbianADOL-C.cpp:246
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::detH
adouble detH
Definition: EshelbianADOL-C.cpp:362
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
EshelbianPlasticity::HMHStVenantKirchhoff
Definition: EshelbianADOL-C.cpp:184
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
EshelbianPlasticity::HMHHencky::mField
MoFEM::Interface & mField
Definition: EshelbianADOL-C.cpp:630
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
EshelbianPlasticity::HMHStVenantKirchhoff::lambda
double lambda
Definition: EshelbianADOL-C.cpp:225
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
EshelbianPlasticity::HMHStVenantKirchhoff::energy
adouble energy
Definition: EshelbianADOL-C.cpp:236
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
EshelbianPlasticity::HMHStVenantKirchhoff::tE
ATensor2 tE
Definition: EshelbianADOL-C.cpp:245
EshelbianPlasticity::HMHHencky::extractBlockData
MoFEMErrorCode extractBlockData(Sev sev)
Definition: EshelbianADOL-C.cpp:547
EshelbianPlasticity::HMHHencky::HMHHencky
HMHHencky(MoFEM::Interface &m_field, const double E, const double nu)
Definition: EshelbianADOL-C.cpp:467
FTensor::Index< 'i', 3 >
convert.n
n
Definition: convert.py:82
EshelbianPlasticity::HMHStVenantKirchhoff::HMHStVenantKirchhoff
HMHStVenantKirchhoff(const double lambda, const double mu)
Definition: EshelbianADOL-C.cpp:189
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1511
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
adouble
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::lambda
double lambda
Definition: EshelbianADOL-C.cpp:355
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:20
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
EshelbianPlasticity::HMHStVenantKirchhoff::getOptions
MoFEMErrorCode getOptions()
Definition: EshelbianADOL-C.cpp:193
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:20
EshelbianPlasticity::HMHStVenantKirchhoff::th
ATensor2 th
Definition: EshelbianADOL-C.cpp:229
EshelbianPlasticity::HMHHencky::BlockData::shearModulusG
double shearModulusG
Definition: EshelbianADOL-C.cpp:634
EshelbianPlasticity::OpJacobian
Definition: EshelbianPlasticity.hpp:325
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::detF
adouble detF
Definition: EshelbianADOL-C.cpp:363
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tF
ATensor2 tF
Definition: EshelbianADOL-C.cpp:360
EshelbianPlasticity::HMHStVenantKirchhoff::f
adouble f
Definition: EshelbianADOL-C.cpp:235
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::B
adouble B
Definition: EshelbianADOL-C.cpp:380
EshelbianPlasticity::HMHStVenantKirchhoff::tH
ATensor2 tH
Definition: EshelbianADOL-C.cpp:230
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tCof
ATensor2 tCof
Definition: EshelbianADOL-C.cpp:370
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tBj
adouble tBj
Definition: EshelbianADOL-C.cpp:375
third
constexpr double third
Definition: EshelbianADOL-C.cpp:14
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:575
EshelbianPlasticity::HMHStVenantKirchhoff::recordTape
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
Definition: EshelbianADOL-C.cpp:250
EshelbianPlasticity::HMHHencky::BlockData::blockEnts
Range blockEnts
Definition: EshelbianADOL-C.cpp:635
EshelbianPlasticity::HMHStVenantKirchhoff::tSigma
ATensor2 tSigma
Definition: EshelbianADOL-C.cpp:232
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:20
EshelbianPlasticity::HMHStVenantKirchhoff::mu
double mu
Definition: EshelbianADOL-C.cpp:226
EshelbianPlasticity::HMHStVenantKirchhoff::tC
ATensor2 tC
Definition: EshelbianADOL-C.cpp:244
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::phi
adouble phi
Definition: EshelbianADOL-C.cpp:378
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63
Definition: EshelbianADOL-C.cpp:314
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::tPulledP
ATensor2 tPulledP
Definition: EshelbianADOL-C.cpp:382
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MU
#define MU(E, NU)
Definition: fem_tools.h:23
EshelbianPlasticity::HMHHencky::OpHenckyJacobian::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: EshelbianADOL-C.cpp:482
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::epsilon
double epsilon
Definition: EshelbianADOL-C.cpp:356
EshelbianPlasticity::OpHMHH::OpHMHH
OpHMHH(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
Definition: EshelbianADOL-C.cpp:19
Generated by Doxygen 1.8.17 and hosted at University of Glasgow