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