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