v0.13.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>
7using namespace MoFEM;
8
9#include <BasicFiniteElements.hpp>
10
11#include <EshelbianPlasticity.hpp>
12#include <boost/math/constants/constants.hpp>
13
14constexpr double third = boost::math::constants::third<double>();
15
16namespace EshelbianPlasticity {
17
18struct 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
29MoFEMErrorCode 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
93MoFEMErrorCode 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
231struct 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
426struct 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 // Look at equation 137 from Bonet Gil paper.
579 energy = alpha * pow(A, 2) + beta * pow(B, 2);
580 energy += (-12 * alpha - 24 * beta) * log(detF) +
581 (lambda / (2 * epsilon * epsilon)) *
582 (pow(detF, epsilon) + pow(detF, -epsilon));
583
584 // Stress Eshelby
585 tSigma(I, J) = -tF(i, I) * tP(i, J);
586 tSigma(N0, N0) += energy;
587 tSigma(N1, N1) += energy;
588 tSigma(N2, N2) += energy;
589
590 // Deviator Cauchy Stress
591 meanCauchy = third * tCauchy(i, i);
592 tDevCauchy(i, j) = tCauchy(i, j);
593 tDevCauchy(N0, N0) -= meanCauchy;
594 tDevCauchy(N1, N1) -= meanCauchy;
595 tDevCauchy(N2, N2) -= meanCauchy;
596
597 // // Plastic surface
598
599 s2 = tDevCauchy(i, j) * tDevCauchy(i, j); // s:s
600 f = sqrt(1.5 * s2);
601 phi = f - sigmaY;
602
603 // Flow
604 tPhi_dDevCauchy(i, j) = tDevCauchy(i, j);
605 tPhi_dDevCauchy(i, j) *= 1.5 / f;
606 tPhi_dCauchy(i, j) = tPhi_dDevCauchy(i, j);
607 tPhi_dCauchy(N0, N0) -= third * tPhi_dDevCauchy(i, i);
608 tPhi_dCauchy(N1, N1) -= third * tPhi_dDevCauchy(i, i);
609 tPhi_dCauchy(N2, N2) -= third * tPhi_dDevCauchy(i, i);
610 tPhi_dSigma(I, J) = tInvF(I, i) * (tPhi_dCauchy(i, j) * tF(j, J));
611 tPhi_dSigma(I, J) *= -(1 / detF);
612
613 // Pull back
614 tPulledP(i, J) = tP(i, I) * tInvH(J, I);
615 tPulledP(i, J) *= detH;
616 tPulledSigma(i, J) = tSigma(i, I) * tInvH(J, I);
617 tPulledSigma(i, J) *= detH;
618 tPulledPhi_dSigma(i, J) = tPhi_dSigma(i, I) * tInvH(J, I);
619 tPulledPhi_dSigma(i, J) *= detH;
620
621 // Set dependent variables to ADOL-C
622 tPulledP(i, j) >>= r_P(i, j);
623 tPulledSigma(i, j) >>= r_Sigma(i, j);
624 phi >>= r_phi;
625 tPulledPhi_dSigma(i, j) >>= r_Flow(i, j);
626
627 trace_off();
628
629 MoFEMFunctionReturn(0);
630 }
631};
632
633MoFEMErrorCode EshelbianCore::addMaterial_HMHHStVenantKirchhoff(
634 const int tape, const double lambda, const double mu,
635 const double sigma_y) {
636 MoFEMFunctionBegin;
637 physicalEquations = boost::shared_ptr<HMHStVenantKirchhoff>(
638 new HMHStVenantKirchhoff(lambda, mu, sigma_y));
639 CHKERR physicalEquations->recordTape(tape, nullptr);
640 MoFEMFunctionReturn(0);
641}
642
643MoFEMErrorCode EshelbianCore::addMaterial_HMHMooneyRivlin(
644 const int tape, const double alpha, const double beta, const double lambda,
645 const double sigma_y) {
646 MoFEMFunctionBegin;
647 physicalEquations = boost::shared_ptr<HMHPMooneyRivlinWriggersEq63>(
648 new HMHPMooneyRivlinWriggersEq63(alpha, beta, lambda, sigma_y));
649 CHKERR physicalEquations->recordTape(tape, nullptr);
650 MoFEMFunctionReturn(0);
651}
652
653}; // namespace EshelbianPlasticity
J
static Index< 'J', 3 > J
Definition: BasicFeTools.hpp:93
N2
static Number< 2 > N2
Definition: BasicFeTools.hpp:101
I
static Index< 'I', 3 > I
Definition: BasicFeTools.hpp:92
N1
static Number< 1 > N1
Definition: BasicFeTools.hpp:100
K
static Index< 'K', 3 > K
Definition: BasicFeTools.hpp:94
N0
static Number< 0 > N0
Definition: BasicFeTools.hpp:99
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:346
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MU
#define MU(E, NU)
Definition: fem_tools.h:23
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65
lambda
constexpr double lambda
Definition: free_surface.cpp:84
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EshelbianPlasticity::EntData
DataForcesAndSourcesCore::EntData EntData
Definition: EshelbianPlasticity.hpp:28
EshelbianPlasticity::getFTensor3FromMat
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
Definition: EshelbianPlasticity.hpp:19
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:76
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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:1204
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1193
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:21
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:17
mu
constexpr double mu
Definition: shallow_wave.cpp:64
EshelbianPlasticity::EshelbianCore::addMaterial_HMHHStVenantKirchhoff
MoFEMErrorCode addMaterial_HMHHStVenantKirchhoff(const int tape, const double lambda, const double mu, const double sigma_y)
Definition: EshelbianADOL-C.cpp:633
EshelbianPlasticity::EshelbianCore::addMaterial_HMHMooneyRivlin
MoFEMErrorCode addMaterial_HMHMooneyRivlin(const int tape, const double alpha, const double beta, const double lambda, const double sigma_y)
Definition: EshelbianADOL-C.cpp:643
EshelbianPlasticity::EshelbianCore::physicalEquations
boost::shared_ptr< PhysicalEquations > physicalEquations
Definition: EshelbianPlasticity.hpp:1171
EshelbianPlasticity::HMHPMooneyRivlinWriggersEq63::recordTape
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
Definition: EshelbianADOL-C.cpp:513
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::HMHPMooneyRivlinWriggersEq63
HMHPMooneyRivlinWriggersEq63(const double alpha, const double beta, const double lambda, const double sigma_y)
Definition: EshelbianADOL-C.cpp:431
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::HMHStVenantKirchhoff
HMHStVenantKirchhoff(const double lambda, const double mu, const double sigma_y)
Definition: EshelbianADOL-C.cpp:236
EshelbianPlasticity::HMHStVenantKirchhoff::recordTape
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
Definition: EshelbianADOL-C.cpp:308
EshelbianPlasticity::OpHMHH::evaluateLhs
MoFEMErrorCode evaluateLhs(EntData &data)
Definition: EshelbianADOL-C.cpp:93
EshelbianPlasticity::OpHMHH::evaluateRhs
MoFEMErrorCode evaluateRhs(EntData &data)
Definition: EshelbianADOL-C.cpp:29
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::OpJacobian::physicsPtr
boost::shared_ptr< PhysicalEquations > physicsPtr
material physical equations
Definition: EshelbianPlasticity.hpp:578
EshelbianPlasticity::OpJacobian::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
Definition: EshelbianPlasticity.hpp:576
Generated on Thu Aug 11 2022 23:36:16 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow