v0.10.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 std::string &field_name, const int tag, const bool eval_rhs,
20  const bool eval_lhs, boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
21  boost::shared_ptr<PhysicalEquations> &physics_ptr)
22  : OpJacobian(field_name, tag, eval_rhs, eval_lhs, data_ptr, physics_ptr) {
23  }
24 
25  MoFEMErrorCode evaluateRhs(EntData &data);
26  MoFEMErrorCode evaluateLhs(EntData &data);
27 };
28 
29 MoFEMErrorCode OpHMHH::evaluateRhs(EntData &ent_data) {
30  MoFEMFunctionBegin;
31  FTensor::Index<'i', 3> i;
32  FTensor::Index<'j', 3> j;
33  FTensor::Index<'k', 3> k;
34 
35  const auto nb_integration_pts = ent_data.getN().size1();
36  auto iu = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
37  // FIXME: that should work with material streach
38  auto iG = getFTensor2FromMat<3, 3>(dataAtPts->GAtPts);
39 
40  auto create_data_vec = [nb_integration_pts](auto &v) {
41  v.resize(nb_integration_pts, false);
42  };
43  auto create_data_mat = [nb_integration_pts](auto &m) {
44  m.resize(9, nb_integration_pts, false);
45  };
46 
47  create_data_mat(dataAtPts->PAtPts);
48  create_data_mat(dataAtPts->SigmaAtPts);
49  create_data_vec(dataAtPts->phiAtPts);
50  create_data_mat(dataAtPts->flowAtPts);
51 
52  auto r_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
53  auto r_Sigma = getFTensor2FromMat<3, 3>(dataAtPts->SigmaAtPts);
54  auto r_phi = getFTensor0FromVec(dataAtPts->phiAtPts);
55  auto r_flow = getFTensor2FromMat<3, 3>(dataAtPts->flowAtPts);
56 
57  for (unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
58 
59  auto t_h = physicsPtr->get_h();
60  for (auto ii : {0, 1, 2})
61  for (auto jj : {0, 1, 2})
62  t_h(ii, jj) = iu(ii, jj);
63 
64  physicsPtr->get_H()(i, j) = iG(i, j);
65  for (int ii = 0; ii != 3; ++ii)
66  physicsPtr->get_H()(ii, ii) += 1;
67 
68  // CHKERR physicsPtr->recordTape(tAg, &t_h);
69  int r = ::function(tAg, physicsPtr->dependentVariables.size(),
70  physicsPtr->activeVariables.size(),
71  &physicsPtr->activeVariables[0],
72  &physicsPtr->dependentVariables[0]);
73  if (r < 0) { // function is locally analytic
74  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
75  "ADOL-C function evaluation with error r = %d", r);
76  }
77 
78  r_P(i, j) = physicsPtr->get_P()(i, j);
79  r_Sigma(i, j) = physicsPtr->get_Sigma()(i, j);
80  r_phi = physicsPtr->get_Phi();
81  r_flow(i, j) = physicsPtr->get_Flow()(i, j);
82 
83  ++iu;
84  ++iG;
85  ++r_P;
86  ++r_Sigma;
87  ++r_phi;
88  ++r_flow;
89  }
90  MoFEMFunctionReturn(0);
91 }
92 
93 MoFEMErrorCode OpHMHH::evaluateLhs(EntData &ent_data) {
94  MoFEMFunctionBegin;
95  FTensor::Index<'i', 3> i;
96  FTensor::Index<'j', 3> j;
97  FTensor::Index<'k', 3> k;
98  const int number_of_active_variables = physicsPtr->activeVariables.size();
99  const int number_of_dependent_variables =
100  physicsPtr->dependentVariables.size();
101  std::vector<double *> jac_ptr(number_of_dependent_variables);
102  for (unsigned int n = 0; n != number_of_dependent_variables; ++n) {
103  jac_ptr[n] =
104  &(physicsPtr
105  ->dependentVariablesDirevatives[n * number_of_active_variables]);
106  }
107 
108  const auto nb_integration_pts = ent_data.getN().size1();
109 
110  auto create_data_mat = [nb_integration_pts](auto &m) {
111  m.resize(9, nb_integration_pts, false);
112  };
113 
114  auto create_data_ten = [nb_integration_pts](auto &m) {
115  m.resize(27, nb_integration_pts, false);
116  };
117 
118  create_data_ten(dataAtPts->P_dh0);
119  create_data_ten(dataAtPts->P_dh1);
120  create_data_ten(dataAtPts->P_dh2);
121  create_data_ten(dataAtPts->P_dH0);
122  create_data_ten(dataAtPts->P_dH1);
123  create_data_ten(dataAtPts->P_dH2);
124  create_data_ten(dataAtPts->Sigma_dh0);
125  create_data_ten(dataAtPts->Sigma_dh1);
126  create_data_ten(dataAtPts->Sigma_dh2);
127  create_data_ten(dataAtPts->Sigma_dH0);
128  create_data_ten(dataAtPts->Sigma_dH1);
129  create_data_ten(dataAtPts->Sigma_dH2);
130  create_data_mat(dataAtPts->phi_dh);
131  create_data_mat(dataAtPts->phi_dH);
132  create_data_ten(dataAtPts->Flow_dh0);
133  create_data_ten(dataAtPts->Flow_dh1);
134  create_data_ten(dataAtPts->Flow_dh2);
135  create_data_ten(dataAtPts->Flow_dH0);
136  create_data_ten(dataAtPts->Flow_dH1);
137  create_data_ten(dataAtPts->Flow_dH2);
138 
139  auto iu = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
140  auto iG = getFTensor2FromMat<3, 3>(dataAtPts->GAtPts);
141 
142  auto r_P_dh0 = getFTensor3FromMat(dataAtPts->P_dh0);
143  auto r_P_dh1 = getFTensor3FromMat(dataAtPts->P_dh1);
144  auto r_P_dh2 = getFTensor3FromMat(dataAtPts->P_dh2);
145  auto r_P_dH0 = getFTensor3FromMat(dataAtPts->P_dH0);
146  auto r_P_dH1 = getFTensor3FromMat(dataAtPts->P_dH1);
147  auto r_P_dH2 = getFTensor3FromMat(dataAtPts->P_dH2);
148  auto r_Sigma_dh0 = getFTensor3FromMat(dataAtPts->Sigma_dh0);
149  auto r_Sigma_dh1 = getFTensor3FromMat(dataAtPts->Sigma_dh1);
150  auto r_Sigma_dh2 = getFTensor3FromMat(dataAtPts->Sigma_dh2);
151  auto r_Sigma_dH0 = getFTensor3FromMat(dataAtPts->Sigma_dH0);
152  auto r_Sigma_dH1 = getFTensor3FromMat(dataAtPts->Sigma_dH1);
153  auto r_Sigma_dH2 = getFTensor3FromMat(dataAtPts->Sigma_dH2);
154  auto r_phi_dh = getFTensor2FromMat<3, 3>(dataAtPts->phi_dh);
155  auto r_phi_dH = getFTensor2FromMat<3, 3>(dataAtPts->phi_dH);
156  auto r_Flow_dh0 = getFTensor3FromMat(dataAtPts->Flow_dh0);
157  auto r_Flow_dh1 = getFTensor3FromMat(dataAtPts->Flow_dh1);
158  auto r_Flow_dh2 = getFTensor3FromMat(dataAtPts->Flow_dh2);
159  auto r_Flow_dH0 = getFTensor3FromMat(dataAtPts->Flow_dH0);
160  auto r_Flow_dH1 = getFTensor3FromMat(dataAtPts->Flow_dH1);
161  auto r_Flow_dH2 = getFTensor3FromMat(dataAtPts->Flow_dH2);
162 
163  for (unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
164 
165  auto t_h = physicsPtr->get_h();
166  for(auto ii : {0, 1, 2})
167  for(auto jj : {0, 1, 2})
168  t_h(ii, jj) = iu(ii, jj);
169 
170  physicsPtr->get_H()(i, j) = iG(i, j);
171  for (int ii = 0; ii != 3; ++ii)
172  physicsPtr->get_H()(ii, ii) += 1;
173 
174  // play recorder for jacobians
175  // CHKERR physicsPtr->recordTape(tAg, &t_h);
176  int r = ::jacobian(tAg, number_of_dependent_variables,
177  number_of_active_variables,
178  &physicsPtr->activeVariables[0], &jac_ptr[0]);
179  if (r < 0) {
180  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
181  "ADOL-C function evaluation with error");
182  }
183 
184  r_P_dh0(i, j, k) = physicsPtr->get_P_dh0()(i, j, k);
185  r_P_dh1(i, j, k) = physicsPtr->get_P_dh1()(i, j, k);
186  r_P_dh2(i, j, k) = physicsPtr->get_P_dh2()(i, j, k);
187  r_P_dH0(i, j, k) = physicsPtr->get_P_dH0()(i, j, k);
188  r_P_dH1(i, j, k) = physicsPtr->get_P_dH1()(i, j, k);
189  r_P_dH2(i, j, k) = physicsPtr->get_P_dH2()(i, j, k);
190  r_Sigma_dh0(i, j, k) = physicsPtr->get_Sigma_dh0()(i, j, k);
191  r_Sigma_dh1(i, j, k) = physicsPtr->get_Sigma_dh1()(i, j, k);
192  r_Sigma_dh2(i, j, k) = physicsPtr->get_Sigma_dh2()(i, j, k);
193  r_Sigma_dH0(i, j, k) = physicsPtr->get_Sigma_dH0()(i, j, k);
194  r_Sigma_dH1(i, j, k) = physicsPtr->get_Sigma_dH1()(i, j, k);
195  r_Sigma_dH2(i, j, k) = physicsPtr->get_Sigma_dH2()(i, j, k);
196  r_phi_dh(i, j) = physicsPtr->get_Phi_dh()(i, j);
197  r_phi_dH(i, j) = physicsPtr->get_Phi_dH()(i, j);
198  r_Flow_dh0(i, j, k) = physicsPtr->get_Flow_dh0()(i, j, k);
199  r_Flow_dh1(i, j, k) = physicsPtr->get_Flow_dh1()(i, j, k);
200  r_Flow_dh2(i, j, k) = physicsPtr->get_Flow_dh2()(i, j, k);
201  r_Flow_dH0(i, j, k) = physicsPtr->get_Flow_dH0()(i, j, k);
202  r_Flow_dH1(i, j, k) = physicsPtr->get_Flow_dH1()(i, j, k);
203  r_Flow_dH2(i, j, k) = physicsPtr->get_Flow_dH2()(i, j, k);
204 
205  ++iu;
206  ++iG;
207  ++r_P_dh0;
208  ++r_P_dh1;
209  ++r_P_dh2;
210  ++r_P_dH0;
211  ++r_P_dH1;
212  ++r_P_dH2;
213  ++r_Sigma_dh0;
214  ++r_Sigma_dh1;
215  ++r_Sigma_dh2;
216  ++r_Sigma_dH0;
217  ++r_Sigma_dH1;
218  ++r_Sigma_dH2;
219  ++r_phi_dh;
220  ++r_phi_dH;
221  ++r_Flow_dh0;
222  ++r_Flow_dh1;
223  ++r_Flow_dh2;
224  ++r_Flow_dH0;
225  ++r_Flow_dH1;
226  ++r_Flow_dH2;
227  }
228  MoFEMFunctionReturn(0);
229 }
230 
231 struct HMHStVenantKirchhoff : public PhysicalEquations {
232 
233  static constexpr int numberOfActiveVariables = 9 + 9;
234  static constexpr int numberOfDependentVariables = 9 + 9 + 1 + 9;
235 
236  HMHStVenantKirchhoff(const double lambda, const double mu,
237  const double sigma_y)
238  : PhysicalEquations(numberOfActiveVariables, numberOfDependentVariables),
239  lambda(lambda), mu(mu), sigmaY(sigma_y) {}
240 
241  MoFEMErrorCode getOptions() {
242  MoFEMFunctionBegin;
243 
244  double E = 1;
245  double nu = 0;
246 
247  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "stvenant_", "", "none");
248 
249  CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
250  PETSC_NULL);
251  CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
252  PETSC_NULL);
253  CHKERR PetscOptionsScalar("-sigmaY", "plastic strain", "", sigmaY, &sigmaY,
254  PETSC_NULL);
255 
256  ierr = PetscOptionsEnd();
257  CHKERRG(ierr);
258 
259  lambda = LAMBDA(E, nu);
260  mu = MU(E, nu);
261 
262  MoFEMFunctionReturn(0);
263  }
264 
265  virtual OpJacobian *
266  returnOpJacobian(const std::string &field_name, const int tag,
267  const bool eval_rhs, const bool eval_lhs,
268  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
269  boost::shared_ptr<PhysicalEquations> &physics_ptr) {
270  return (
271  new OpHMHH(field_name, tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
272  }
273 
274  double lambda;
275  double mu;
276  double sigmaY;
277 
278  ATensor2 th;
279  ATensor2 tH;
280  ATensor2 tP;
281  ATensor2 tSigma;
282  adouble phi;
283  adouble s2;
284  adouble f;
285  adouble energy;
286 
287  adouble detH;
288  adouble detF;
289  adouble trE;
290  adouble meanCauchy;
291  ATensor2 tInvH;
292  ATensor2 tF;
293  ATensor2 tInvF;
294  ATensor2 tC;
295  ATensor2 tE;
296  ATensor2 tS;
297  ATensor2 tCauchy;
298  ATensor2 tDevCauchy;
299 
300  ATensor2 tPhi_dDevCauchy;
301  ATensor2 tPhi_dCauchy;
302  ATensor2 tPhi_dSigma;
303 
304  ATensor2 tPulledP;
305  ATensor2 tPulledSigma;
306  ATensor2 tPulledPhi_dSigma;
307 
308  MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr) {
309  MoFEMFunctionBegin;
310  FTensor::Index<'i', 3> i;
311  FTensor::Index<'j', 3> j;
312  FTensor::Index<'I', 3> I;
313  FTensor::Index<'J', 3> J;
314  FTensor::Number<0> N0;
315  FTensor::Number<1> N1;
316  FTensor::Number<2> N2;
317 
318  CHKERR getOptions();
319 
320  auto ih = get_h();
321  auto iH = get_H();
322  if (t_h_ptr)
323  ih(i, j) = (*t_h_ptr)(i, j);
324  else {
325  ih(i, j) = 0;
326  for (auto ii : {0, 1, 2})
327  ih(ii, ii) = 1;
328  }
329 
330  iH(i, j) = 0;
331  for (auto ii : {0, 1, 2})
332  iH(ii, ii) = 1;
333 
334  auto r_P = get_P();
335  auto r_Sigma = get_Sigma();
336  double &r_phi = get_PhiRef();
337  auto r_Flow = get_Flow();
338 
339  enableMinMaxUsingAbs();
340  trace_on(tape);
341 
342  // Set active variables to ADOL-C
343  th(i, j) <<= get_h()(i, j);
344  tH(i, j) <<= get_H()(i, j);
345 
346  // Deformation gradient
347  CHKERR determinantTensor3by3(tH, detH);
348  CHKERR invertTensor3by3(tH, detH, tInvH);
349  tF(i, I) = th(i, J) * tInvH(J, I);
350  CHKERR determinantTensor3by3(tF, detF);
351  CHKERR invertTensor3by3(tF, detF, tInvF);
352 
353  // Deformation and strain
354  tC(I, J) = tF(i, I) * tF(i, J);
355  tE(I, J) = tC(I, J);
356  tE(N0, N0) -= 1;
357  tE(N1, N1) -= 1;
358  tE(N2, N2) -= 1;
359  tE(I, J) *= 0.5;
360 
361  // Energy
362  trE = tE(I, I);
363  energy = 0.5 * lambda * trE * trE + mu * (tE(I, J) * tE(I, J));
364 
365  // Stress Piola II
366  tS(I, J) = tE(I, J);
367  tS(I, J) *= 2 * mu;
368  tS(N0, N0) += lambda * trE;
369  tS(N1, N1) += lambda * trE;
370  tS(N2, N2) += lambda * trE;
371  // Stress Piola I
372  tP(i, J) = tF(i, I) * tS(I, J);
373  // Stress Cauchy
374  tCauchy(i, j) = tP(i, J) * tF(j, J);
375  tCauchy(i, j) *= (1 / detF);
376 
377  // Stress Eshelby
378  tSigma(I, J) = -tF(i, I) * tP(i, J);
379  tSigma(N0, N0) += energy;
380  tSigma(N1, N1) += energy;
381  tSigma(N2, N2) += energy;
382 
383  // Deviator Cauchy Stress
384  meanCauchy = third * tCauchy(i, i);
385  tDevCauchy(i, j) = tCauchy(i, j);
386  tDevCauchy(N0, N0) -= meanCauchy;
387  tDevCauchy(N1, N1) -= meanCauchy;
388  tDevCauchy(N2, N2) -= meanCauchy;
389 
390  // Plastic surface
391 
392  s2 = tDevCauchy(i, j) * tDevCauchy(i, j); // s:s
393  f = sqrt(1.5 * s2);
394  phi = f - sigmaY;
395 
396  // Flow
397  tPhi_dDevCauchy(i, j) = tDevCauchy(i, j);
398  tPhi_dDevCauchy(i, j) *= 1.5 / f;
399  tPhi_dCauchy(i, j) = tPhi_dDevCauchy(i, j);
400  tPhi_dCauchy(N0, N0) -= third * tPhi_dDevCauchy(i, i);
401  tPhi_dCauchy(N1, N1) -= third * tPhi_dDevCauchy(i, i);
402  tPhi_dCauchy(N2, N2) -= third * tPhi_dDevCauchy(i, i);
403  tPhi_dSigma(I, J) = tInvF(I, i) * (tPhi_dCauchy(i, j) * tF(j, J));
404  tPhi_dSigma(I, J) *= -(1 / detF);
405 
406  // Pull back
407  tPulledP(i, J) = tP(i, I) * tInvH(J, I);
408  tPulledP(i, J) *= detH;
409  tPulledSigma(i, J) = tSigma(i, I) * tInvH(J, I);
410  tPulledSigma(i, J) *= detH;
411  tPulledPhi_dSigma(i, J) = tPhi_dSigma(i, I) * tInvH(J, I);
412  tPulledPhi_dSigma(i, J) *= detH;
413 
414  // Set dependent variables to ADOL-C
415  tPulledP(i, j) >>= r_P(i, j);
416  tPulledSigma(i, j) >>= r_Sigma(i, j);
417  phi >>= r_phi;
418  tPulledPhi_dSigma(i, j) >>= r_Flow(i, j);
419 
420  trace_off();
421 
422  MoFEMFunctionReturn(0);
423  }
424 };
425 
426 struct HMHPMooneyRivlinWriggersEq63 : public PhysicalEquations {
427 
428  static constexpr int numberOfActiveVariables = 9 + 9;
429  static constexpr int numberOfDependentVariables = 9 + 9 + 1 + 9;
430 
431  HMHPMooneyRivlinWriggersEq63(const double alpha, const double beta,
432  const double lambda, const double sigma_y)
433  : PhysicalEquations(numberOfActiveVariables, numberOfDependentVariables),
434  alpha(alpha), beta(beta), lambda(lambda), epsilon(0), sigmaY(sigma_y) {}
435 
436  virtual OpJacobian *
437  returnOpJacobian(const std::string &field_name, const int tag,
438  const bool eval_rhs, const bool eval_lhs,
439  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
440  boost::shared_ptr<PhysicalEquations> &physics_ptr) {
441  return (
442  new OpHMHH(field_name, tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
443  }
444 
445  MoFEMErrorCode getOptions() {
446  MoFEMFunctionBegin;
447  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "mooneyrivlin_", "", "none");
448 
449  alpha = 1;
450  CHKERR PetscOptionsScalar("-alpha", "Alpha", "", alpha, &alpha, PETSC_NULL);
451  beta = 1;
452  CHKERR PetscOptionsScalar("-beta", "Beta", "", beta, &beta, PETSC_NULL);
453 
454  lambda = 1;
455  CHKERR PetscOptionsScalar("-lambda", "Lambda", "", lambda, &lambda,
456  PETSC_NULL);
457 
458  epsilon = 0;
459  CHKERR PetscOptionsScalar("-epsilon", "Epsilon", "", epsilon, &epsilon,
460  PETSC_NULL);
461 
462  CHKERR PetscOptionsScalar("-sigma_y", "plastic strain", "", sigmaY, &sigmaY,
463  PETSC_NULL);
464 
465  ierr = PetscOptionsEnd();
466  CHKERRG(ierr);
467  MoFEMFunctionReturn(0);
468  }
469 
470  double alpha;
471  double beta;
472  double lambda;
473  double epsilon;
474  double sigmaY;
475 
476  ATensor2 th;
477  ATensor2 tH;
478  ATensor2 tF;
479 
480  adouble detH;
481  adouble detF;
482  ATensor2 tInvH;
483  ATensor2 tInvF;
484 
485  ATensor2 tP;
486  ATensor2 tSigma;
487  ATensor2 tCauchy;
488  ATensor2 tDevCauchy;
489 
490  ATensor2 tCof;
491  // ATensor3 tCofT1;
492  // ATensor3 tCofT2;
493  ATensor2 tBF;
494  ATensor2 tBCof;
495  adouble tBj;
496 
497  adouble energy;
498  adouble meanCauchy;
499  adouble s2;
500  adouble f;
501  adouble phi;
502  adouble A;
503  adouble B;
504 
505  ATensor2 tPhi_dDevCauchy;
506  ATensor2 tPhi_dCauchy;
507  ATensor2 tPhi_dSigma;
508 
509  ATensor2 tPulledP;
510  ATensor2 tPulledSigma;
511  ATensor2 tPulledPhi_dSigma;
512 
513  MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr) {
514  MoFEMFunctionBegin;
515 
516  FTensor::Index<'i', 3> i;
517  FTensor::Index<'j', 3> j;
518  FTensor::Index<'k', 3> k;
519  FTensor::Index<'I', 3> I;
520  FTensor::Index<'J', 3> J;
521  FTensor::Index<'K', 3> K;
522  FTensor::Number<0> N0;
523  FTensor::Number<1> N1;
524  FTensor::Number<2> N2;
525 
526  CHKERR getOptions();
527 
528  auto ih = get_h();
529  auto iH = get_H();
530  if (t_h_ptr)
531  ih(i, j) = (*t_h_ptr)(i, j);
532  else {
533  ih(i, j) = 0;
534  for (auto ii : {0, 1, 2})
535  ih(ii, ii) = 1;
536  }
537 
538  iH(i, j) = 0;
539  for (auto ii : {0, 1, 2})
540  iH(ii, ii) = 1;
541 
542  auto r_P = get_P();
543  auto r_Sigma = get_Sigma();
544  double &r_phi = get_PhiRef();
545  auto r_Flow = get_Flow();
546 
547  enableMinMaxUsingAbs();
548  trace_on(tape);
549 
550  // Set active variables to ADOL-C
551  th(i, j) <<= get_h()(i, j);
552  tH(i, j) <<= get_H()(i, j);
553 
554  // Deformation gradient
555  CHKERR determinantTensor3by3(tH, detH);
556  CHKERR invertTensor3by3(tH, detH, tInvH);
557 
558  tF(i, I) = th(i, J) * tInvH(J, I);
559  CHKERR determinantTensor3by3(tF, detF);
560  CHKERR invertTensor3by3(tF, detF, tInvF);
561 
562  tCof(i, I) = detF * tInvF(I, i);
563 
564  A = tF(k, K) * tF(k, K);
565  B = tCof(k, K) * tCof(k, K);
566 
567  tBF(i, I) = 4 * alpha * (A * tF(i, I));
568  tBCof(i, I) = 4 * beta * (B * tCof(i, I));
569  tBj = (-12 * alpha - 24 * beta) / detF +
570  0.5 * (lambda / epsilon) *
571  (pow(detF, epsilon - 1) - pow(detF, -epsilon - 1));
572 
573  tP(i, I) = tBF(i, I);
574  tP(i, I) += (levi_civita(i, j, k) * tBCof(j, J)) *
575  (levi_civita(I, J, K) * tF(k, K));
576  tP(i, I) += tCof(i, I) * tBj;
577 
578  energy = alpha * pow(A, 2) + beta * pow(B, 2);
579  energy += (-12 * alpha - 24 * beta) * log(detF) +
580  (lambda / (2 * epsilon * epsilon)) *
581  (pow(detF, epsilon) + pow(detF, -epsilon));
582 
583  // Stress Eshelby
584  tSigma(I, J) = -tF(i, I) * tP(i, J);
585  tSigma(N0, N0) += energy;
586  tSigma(N1, N1) += energy;
587  tSigma(N2, N2) += energy;
588 
589  // Deviator Cauchy Stress
590  meanCauchy = third * tCauchy(i, i);
591  tDevCauchy(i, j) = tCauchy(i, j);
592  tDevCauchy(N0, N0) -= meanCauchy;
593  tDevCauchy(N1, N1) -= meanCauchy;
594  tDevCauchy(N2, N2) -= meanCauchy;
595 
596  // // Plastic surface
597 
598  s2 = tDevCauchy(i, j) * tDevCauchy(i, j); // s:s
599  f = sqrt(1.5 * s2);
600  phi = f - sigmaY;
601 
602  // Flow
603  tPhi_dDevCauchy(i, j) = tDevCauchy(i, j);
604  tPhi_dDevCauchy(i, j) *= 1.5 / f;
605  tPhi_dCauchy(i, j) = tPhi_dDevCauchy(i, j);
606  tPhi_dCauchy(N0, N0) -= third * tPhi_dDevCauchy(i, i);
607  tPhi_dCauchy(N1, N1) -= third * tPhi_dDevCauchy(i, i);
608  tPhi_dCauchy(N2, N2) -= third * tPhi_dDevCauchy(i, i);
609  tPhi_dSigma(I, J) = tInvF(I, i) * (tPhi_dCauchy(i, j) * tF(j, J));
610  tPhi_dSigma(I, J) *= -(1 / detF);
611 
612  // Pull back
613  tPulledP(i, J) = tP(i, I) * tInvH(J, I);
614  tPulledP(i, J) *= detH;
615  tPulledSigma(i, J) = tSigma(i, I) * tInvH(J, I);
616  tPulledSigma(i, J) *= detH;
617  tPulledPhi_dSigma(i, J) = tPhi_dSigma(i, I) * tInvH(J, I);
618  tPulledPhi_dSigma(i, J) *= detH;
619 
620  // Set dependent variables to ADOL-C
621  tPulledP(i, j) >>= r_P(i, j);
622  tPulledSigma(i, j) >>= r_Sigma(i, j);
623  phi >>= r_phi;
624  tPulledPhi_dSigma(i, j) >>= r_Flow(i, j);
625 
626  trace_off();
627 
628  MoFEMFunctionReturn(0);
629  }
630 };
631 
632 MoFEMErrorCode EshelbianCore::addMaterial_HMHHStVenantKirchhoff(
633  const int tape, const double lambda, const double mu,
634  const double sigma_y) {
635  MoFEMFunctionBegin;
636  physicalEquations = boost::shared_ptr<HMHStVenantKirchhoff>(
637  new HMHStVenantKirchhoff(lambda, mu, sigma_y));
638  CHKERR physicalEquations->recordTape(tape, nullptr);
639  MoFEMFunctionReturn(0);
640 }
641 
642 MoFEMErrorCode EshelbianCore::addMaterial_HMHMooneyRivlin(
643  const int tape, const double alpha, const double beta, const double lambda,
644  const double sigma_y) {
645  MoFEMFunctionBegin;
646  physicalEquations = boost::shared_ptr<HMHPMooneyRivlinWriggersEq63>(
647  new HMHPMooneyRivlinWriggersEq63(alpha, beta, lambda, sigma_y));
648  CHKERR physicalEquations->recordTape(tape, nullptr);
649  MoFEMFunctionReturn(0);
650 }
651 
652 }; // namespace EshelbianPlasticity
m
static Index< 'm', 3 > m
Definition: BasicFeTools.hpp:83
n
static Index< 'n', 3 > n
Definition: BasicFeTools.hpp:84
J
static Index< 'J', 3 > J
Definition: BasicFeTools.hpp:89
N2
static Number< 2 > N2
Definition: BasicFeTools.hpp:97
I
static Index< 'I', 3 > I
Definition: BasicFeTools.hpp:88
N1
static Number< 1 > N1
Definition: BasicFeTools.hpp:96
K
static Index< 'K', 3 > K
Definition: BasicFeTools.hpp:90
N0
static Number< 0 > N0
Definition: BasicFeTools.hpp:95
third
constexpr double third
Definition: EshelbianADOL-C.cpp:14
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:126
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
MU
#define MU(E, NU)
Definition: fem_tools.h:33
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:32
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
i
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EshelbianPlasticity::getFTensor3FromMat
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
Definition: EshelbianPlasticity.hpp:24
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
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:886
MoFEM::getFTensor0FromVec
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:143
MoFEM::determinantTensor3by3
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:905
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33
sigmaY
double sigmaY
Definition: plastic.cpp:106
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::returnOpJacobian
virtual OpJacobian * returnOpJacobian(const std::string &field_name, 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:437
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::recordTape
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
Definition: EshelbianADOL-C.cpp:513
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::HMHPMooneyRivlinWriggersEq63
HMHPMooneyRivlinWriggersEq63(const double alpha, const double beta, const double lambda, const double sigma_y)
Definition: EshelbianADOL-C.cpp:431
EshelbianPlasticity::HMHStVenantKirchhoff::HMHStVenantKirchhoff
HMHStVenantKirchhoff(const double lambda, const double mu, const double sigma_y)
Definition: EshelbianADOL-C.cpp:236
EshelbianPlasticity::HMHStVenantKirchhoff::returnOpJacobian
virtual OpJacobian * returnOpJacobian(const std::string &field_name, 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:266
EshelbianPlasticity::HMHStVenantKirchhoff::recordTape
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
Definition: EshelbianADOL-C.cpp:308
EshelbianPlasticity::OpHMHH::OpHMHH
OpHMHH(const std::string &field_name, 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
MoFEM::DataForcesAndSourcesCore::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: DataStructures.hpp:60
MoFEM::DataForcesAndSourcesCore::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: DataStructures.hpp:1288
Generated on Tue Feb 16 2021 20:34:55 for MoFEM by Doxygen 1.9.1 and hosted at University of Glasgow