Line data Source code
1 : #ifndef MODEL_SOLID_FINITE_NEOHOOKEAN_H_
2 : #define MODEL_SOLID_FINITE_NEOHOOKEAN_H_
3 :
4 : #include "IO/ParmParse.H"
5 : #include "Model/Solid/Solid.H"
6 :
7 : namespace Model
8 : {
9 : namespace Solid
10 : {
11 : namespace Finite
12 : {
13 : class NeoHookean : public Solid<Set::Sym::Major>
14 : {
15 : public:
16 39755 : NeoHookean() {};
17 : NeoHookean(Solid<Set::Sym::Major> base) : Solid<Set::Sym::Major>(base) {};
18 56969 : virtual ~NeoHookean() {};
19 :
20 960 : Set::Scalar W(const Set::Matrix& a_F) const override
21 : {
22 : #if AMREX_SPACEDIM==2
23 380 : Eigen::Matrix3d F = Eigen::Matrix3d::Identity();
24 380 : F(0, 0) = a_F(0, 0);
25 380 : F(0, 1) = a_F(0, 1);
26 380 : F(1, 0) = a_F(1, 0);
27 380 : F(1, 1) = a_F(1, 1);
28 : #elif AMREX_SPACEDIM==3
29 580 : Eigen::Matrix3d F = a_F;
30 : #endif
31 :
32 960 : Set::Scalar J = F.determinant();
33 960 : Set::Scalar J23 = std::pow(fabs(J), 2. / 3.);
34 960 : Set::Scalar w = 0.0;
35 960 : w += 0.5 * mu * ((F * F.transpose()).trace() / J23 - 3.);
36 960 : w += 0.5 * kappa * (J - 1.0) * (J - 1.0);
37 960 : return w;
38 : }
39 123488 : Set::Matrix DW(const Set::Matrix& a_F) const override
40 : {
41 : #if AMREX_SPACEDIM==2
42 120228 : Eigen::Matrix3d F = Eigen::Matrix3d::Identity();
43 120228 : F(0, 0) = a_F(0, 0);
44 120228 : F(0, 1) = a_F(0, 1);
45 120228 : F(1, 0) = a_F(1, 0);
46 120228 : F(1, 1) = a_F(1, 1);
47 : #elif AMREX_SPACEDIM==3
48 3260 : Eigen::Matrix3d F = a_F;
49 : #endif
50 :
51 123488 : Set::Scalar J = F.determinant();
52 123488 : Set::Scalar J23 = std::pow(fabs(J), 2. / 3.);
53 123488 : Eigen::Matrix3d FinvT = F.inverse().transpose();
54 :
55 123488 : Eigen::Matrix3d dw = Eigen::Matrix3d::Zero();
56 :
57 123488 : dw += mu * (F / J23 - (F * F.transpose()).trace() * FinvT / (3. * J23));
58 123488 : dw += kappa * (J - 1) * J * FinvT;
59 :
60 : #if AMREX_SPACEDIM==2
61 120228 : Set::Matrix r_dw;
62 120228 : r_dw(0, 0) = dw(0, 0);
63 120228 : r_dw(0, 1) = dw(0, 1);
64 120228 : r_dw(1, 0) = dw(1, 0);
65 120228 : r_dw(1, 1) = dw(1, 1);
66 240456 : return r_dw;
67 : #elif AMREX_SPACEDIM==3
68 6520 : return dw;
69 : #endif
70 : }
71 71074 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Major> DDW(const Set::Matrix& a_F) const override
72 : {
73 : #if AMREX_SPACEDIM==2
74 71054 : Eigen::Matrix3d F = Eigen::Matrix3d::Identity();
75 71054 : F(0, 0) = a_F(0, 0);
76 71054 : F(0, 1) = a_F(0, 1);
77 71054 : F(1, 0) = a_F(1, 0);
78 71054 : F(1, 1) = a_F(1, 1);
79 : #elif AMREX_SPACEDIM==3
80 20 : Eigen::Matrix3d F = a_F;
81 : #endif
82 :
83 71074 : Set::Matrix4<3, Set::Sym::Major> ddw;
84 71074 : Set::Scalar J = F.determinant();
85 71074 : Set::Scalar J23 = std::pow(fabs(J), 2. / 3.);
86 71074 : Eigen::Matrix3d FinvT = F.inverse().transpose();
87 284296 : for (int i = 0; i < 3; i++)
88 852888 : for (int j = 0; j < 3; j++)
89 2558660 : for (int k = 0; k < 3; k++)
90 7675990 : for (int l = 0; l < 3; l++)
91 : {
92 5756990 : ddw(i, j, k, l) = 0.0;
93 :
94 5756990 : Set::Scalar t1 = 0.0, t2 = 0.0;
95 :
96 5756990 : if (i == k && j == l) t1 += 1.0;
97 5756990 : t1 -= (2. / 3.) * F(i, j) * FinvT(k, l);
98 5756990 : t1 -= (2. / 3.) * FinvT(i, j) * F(k, l);
99 5756990 : t1 += (2. / 9.) * (F * F.transpose()).trace() * FinvT(i, j) * FinvT(k, l);
100 5756990 : t1 += (1. / 3.) * (F * F.transpose()).trace() * FinvT(i, l) * FinvT(k, j);
101 :
102 5756990 : t2 += (2. * J - 1.) * FinvT(i, j) * FinvT(k, l);
103 5756990 : t2 += (1. - J) * FinvT(i, l) * FinvT(k, j);
104 :
105 11513940 : ddw(i, j, k, l) = (mu / J23) * t1 + kappa * J * t2;
106 : }
107 : #if AMREX_SPACEDIM==2
108 71054 : Set::Matrix4<2, Set::Sym::Major> r_ddw;
109 213162 : for (int i = 0; i < 2; i++)
110 426324 : for (int j = 0; j < 2; j++)
111 852648 : for (int k = 0; k < 2; k++)
112 1705300 : for (int l = 0; l < 2; l++)
113 2273730 : r_ddw(i, j, k, l) = ddw(i, j, k, l);
114 142108 : return r_ddw;
115 : #elif AMREX_SPACEDIM==3
116 40 : return ddw;
117 : #endif
118 : }
119 0 : virtual void Print(std::ostream& out) const override
120 : {
121 0 : out << "mu = " << mu << " kappa = " << kappa;
122 0 : }
123 :
124 : public:
125 : Set::Scalar mu = NAN, kappa = NAN;
126 : static constexpr KinematicVariable kinvar = KinematicVariable::F;
127 :
128 : public:
129 8387 : static NeoHookean Zero()
130 : {
131 8387 : NeoHookean ret;
132 8387 : ret.mu = 0.0;
133 8387 : ret.kappa = 0.0;
134 8387 : return ret;
135 : }
136 128 : static NeoHookean Random()
137 : {
138 128 : NeoHookean ret;
139 128 : ret.mu = Util::Random();
140 128 : ret.kappa = Util::Random();
141 128 : return ret;
142 : }
143 4 : static void Parse(NeoHookean& value, IO::ParmParse& pp)
144 : {
145 : //Set::Scalar mu = NAN, kappa = NAN;
146 4 : if (pp.contains("shear") && pp.contains("kappa")) {
147 0 : pp_query("shear", value.mu); // Shear modulus
148 0 : pp_query("kappa", value.kappa); // Bulk modulus
149 : }
150 4 : else if (pp.contains("mu") && pp.contains("kappa")) {
151 4 : pp_query("mu", value.mu); // Alternative input for shear modulus
152 4 : pp_query("kappa", value.kappa); // Bulk modulus
153 : }
154 0 : else if (pp.contains("lame") && pp.contains("shear")) {
155 : Set::Scalar lame;
156 0 : pp_query("shear", value.mu); // Shear modulus
157 0 : pp_query("lame", lame); // Lame parameter
158 0 : value.kappa = lame + (2.0 * value.mu) / 3.0;
159 : }
160 0 : else if (pp.contains("E") && pp.contains("nu")) {
161 : Set::Scalar E, nu;
162 0 : pp_query("E", E); // Young's modulus
163 0 : pp_query("nu", nu); // Poisson's ratio
164 0 : value.kappa = E / (3.0 - 6.0 * nu);
165 0 : value.mu = E / (2.0 + 2.0 * nu);
166 : }
167 4 : }
168 :
169 : #define OP_CLASS NeoHookean
170 : #define OP_VARS X(kappa) X(mu)
171 : #include "Model/Solid/InClassOperators.H"
172 : };
173 : #include "Model/Solid/ExtClassOperators.H"
174 :
175 : }
176 : }
177 : }
178 :
179 : #endif
|
