v0.15.0
Loading...
Searching...
No Matches
HMHPMooneyRivlinWriggersEq63.cpp
Go to the documentation of this file.
1/**
2 * @file HMHPMooneyRivlinWriggersEq63.cpp
3 * @brief Implement of MooneyRivlinWriggersEq63
4 *
5 * @copyright Copyright (c) 2024
6 *
7 */
8
9namespace EshelbianPlasticity {
10
12
13 static constexpr int numberOfActiveVariables = 9;
14 static constexpr int numberOfDependentVariables = 9;
15
20
21 MoFEMErrorCode getOptions() {
23 PetscOptionsBegin(PETSC_COMM_WORLD, "mooneyrivlin_", "", "none");
24 alpha = 1;
25 CHKERR PetscOptionsScalar("-alpha", "Alpha", "", alpha, &alpha, PETSC_NULLPTR);
26 beta = 1;
27 CHKERR PetscOptionsScalar("-beta", "Beta", "", beta, &beta, PETSC_NULLPTR);
28
29 lambda = 1;
30 CHKERR PetscOptionsScalar("-lambda", "Lambda", "", lambda, &lambda,
31 PETSC_NULLPTR);
32
33 epsilon = 0;
34 CHKERR PetscOptionsScalar("-epsilon", "Epsilon", "", epsilon, &epsilon,
35 PETSC_NULLPTR);
36
37 PetscOptionsEnd();
39 }
40
41 double alpha;
42 double beta;
43 double lambda;
44 double epsilon;
45
49
54
57
62
67
69
70 MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr) {
72
73 FTensor::Index<'i', 3> i;
74 FTensor::Index<'j', 3> j;
75 FTensor::Index<'k', 3> k;
76 FTensor::Index<'I', 3> I;
77 FTensor::Index<'J', 3> J;
78 FTensor::Index<'K', 3> K;
82
84
86
87 auto ih = get_h();
88 if (t_h_ptr)
89 ih(i, j) = (*t_h_ptr)(i, j);
90 else {
91 ih(i, j) = t_kd(i, j);
92 }
93
94 auto r_P = get_P();
95
96 enableMinMaxUsingAbs();
97 trace_on(tape);
98
99 // Set active variables to ADOL-C
100 th(i, j) <<= get_h()(i, j);
101
102 tF(i, I) = th(i, I);
103 CHKERR determinantTensor3by3(tF, detF);
104 CHKERR invertTensor3by3(tF, detF, tInvF);
105
106 tCof(i, I) = detF * tInvF(I, i);
107
108 A = tF(k, K) * tF(k, K);
109 B = tCof(k, K) * tCof(k, K);
110
111 tBF(i, I) = 4 * alpha * (A * tF(i, I));
112 tBCof(i, I) = 4 * beta * (B * tCof(i, I));
113 tBj = (-12 * alpha - 24 * beta) / detF +
114 0.5 * (lambda / epsilon) *
115 (pow(detF, epsilon - 1) - pow(detF, -epsilon - 1));
116
117 tP(i, I) = tBF(i, I);
118 tP(i, I) += (levi_civita(i, j, k) * tBCof(j, J)) *
119 (levi_civita(I, J, K) * tF(k, K));
120 tP(i, I) += tCof(i, I) * tBj;
121
122 // Set dependent variables to ADOL-C
123 tP(i, j) >>= r_P(i, j);
124
125 trace_off();
126
128 }
129};
130
131
132
133} // namespace EshelbianPlasticity
Kronecker Delta class symmetric.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr IntegrationType I
HMHPMooneyRivlinWriggersEq63(const double alpha, const double beta, const double lambda)
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)