Line data Source code
1 : #ifndef MODEL_SOLID_AFFINE_ISOTROPIC_H_
2 : #define MODEL_SOLID_AFFINE_ISOTROPIC_H_
3 :
4 : #include "AMReX.H"
5 : #include "IO/ParmParse.H"
6 : #include "Model/Solid/Affine/Affine.H"
7 :
8 : namespace Model
9 : {
10 : namespace Solid
11 : {
12 : namespace Affine
13 : {
14 : class Isotropic : public Affine<Set::Sym::Isotropic>
15 : {
16 : public:
17 :
18 222063 : Isotropic() {};
19 : Isotropic(Affine<Set::Sym::Isotropic> base) : Affine<Set::Sym::Isotropic>(base) {};
20 : Isotropic(Set::Scalar a_mu, Set::Scalar a_lambda, Set::Matrix a_F0 = Set::Matrix::Zero())
21 : {
22 : Define(a_mu, a_lambda, a_F0);
23 : };
24 :
25 38399 : void Define(Set::Scalar a_mu, Set::Scalar a_lambda, Set::Matrix a_F0)
26 : {
27 38399 : F0 = a_F0;
28 38399 : ddw = Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>(a_lambda, a_mu);
29 38399 : }
30 :
31 : void Define(Set::Scalar a_mu, Set::Scalar a_lambda)
32 : {
33 : ddw = Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>(a_lambda, a_mu);
34 : }
35 :
36 260 : Set::Scalar W(const Set::Matrix& F) const override
37 : {
38 520 : return 0.5 * ((F - F0).transpose() * (ddw * ((F - F0)))).trace();
39 : }
40 111902 : Set::Matrix DW(const Set::Matrix& F) const override
41 : {
42 223804 : return ddw * (F - F0);
43 : }
44 51294 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> DDW(const Set::Matrix& /*F*/) const override
45 : {
46 51294 : return ddw;
47 : }
48 0 : virtual void Print(std::ostream& out) const override
49 : {
50 0 : out << "lambda=" << ddw.Lambda() << ", mu=" << ddw.Mu() << ", F0=\n" << F0;
51 0 : }
52 0 : virtual bool ContainsNan() override
53 : {
54 0 : if (F0.hasNaN()) return true;
55 0 : if (ddw.contains_nan()) return true;
56 0 : else return false;
57 : }
58 : public:
59 38349 : static Isotropic Zero()
60 : {
61 38349 : Isotropic ret;
62 38349 : Set::Scalar mu = 0.0;
63 38349 : Set::Scalar lambda = 0.0;
64 38349 : Set::Matrix F0 = Set::Matrix::Zero();
65 38349 : ret.Define(mu, lambda, F0);
66 76698 : return ret;
67 : }
68 44 : static Isotropic Random()
69 : {
70 44 : Isotropic ret;
71 44 : Set::Scalar mu = Util::Random();
72 44 : Set::Scalar lambda = Util::Random();
73 44 : Set::Matrix F0 = Set::Matrix::Random();
74 44 : ret.Define(mu, lambda, F0);
75 88 : return ret;
76 : }
77 6 : static void Parse(Isotropic& value, IO::ParmParse& pp)
78 : {
79 6 : Set::Scalar mu = NAN, lambda = NAN;
80 6 : Set::Matrix F0 = Set::Matrix::Zero();
81 6 : if (pp.contains("lame") && pp.contains("shear"))
82 : {
83 0 : pp_query("lame", lambda); // Lame modulus
84 0 : pp_query("shear", mu); // Shear modulus
85 : }
86 6 : else if (pp.contains("E") && pp.contains("nu"))
87 : {
88 : Set::Scalar E, nu;
89 6 : pp_query("E", E); // Elastic modulus
90 6 : pp_query("nu", nu); // Poisson's ratio
91 6 : lambda = E * nu / (1.0 + nu) / (1.0 - 2.0 * nu);
92 6 : mu = E / 2.0 / (1.0 + nu);
93 : }
94 0 : else if (pp.contains("E") && pp.contains("shear"))
95 : {
96 : Set::Scalar E;
97 0 : pp_query("E", E); // Young's modulus
98 0 : pp_query("shear", mu); // Shear modulus
99 0 : lambda = mu * (E - 2.0 * mu) / (3.0 * mu - E);
100 : }
101 0 : else if (pp.contains("E") && pp.contains("lame"))
102 : {
103 : Set::Scalar E;
104 0 : pp_query("E", E); // Young's modulus
105 0 : pp_query("lame", lambda); // Lame parameter
106 0 : mu = (E - 3.0 * lambda + sqrt(E * E + 9.0 * lambda * lambda + 2.0 * E * lambda)) / 4.0;
107 : }
108 0 : else if (pp.contains("nu") && pp.contains("shear"))
109 : {
110 : Set::Scalar nu;
111 0 : pp_query("nu", nu); // Poisson's ratio
112 0 : pp_query("shear", mu); // Shear modulus
113 0 : lambda = 2.0 * mu * nu / (1.0 - 2.0 * nu);
114 : }
115 0 : else if (pp.contains("lambda") && pp.contains("nu"))
116 : {
117 : Set::Scalar nu;
118 0 : pp_query("lambda", lambda); // Lame parameter
119 0 : pp_query("nu", nu); // Poisson's ratio
120 0 : mu = lambda * (1.0 - 2.0 * nu) / 2.0 / nu;
121 : }
122 0 : else if (pp.contains("bulk") && pp.contains("shear"))
123 : {
124 : Set::Scalar K;
125 0 : pp_query("bulk", K); // Bulk modulus
126 0 : pp_query("shear", mu); // Shear modulus
127 0 : lambda = K - (2.0 * mu / 3.0);
128 : }
129 : else
130 : {
131 0 : Util::Exception(INFO, "Model parameters not specified with either (lame, shear), or (E, nu)");
132 : }
133 6 : if (pp.contains("F0"))
134 : {
135 5 : pp_queryarr("F0", F0); // Eigendeformation gradient
136 : }
137 6 : value.Define(mu, lambda, F0);
138 6 : }
139 : #define OP_CLASS Isotropic
140 : #define OP_VARS X(ddw) X(F0)
141 : #include "Model/Solid/InClassOperators.H"
142 : };
143 : #include "Model/Solid/ExtClassOperators.H"
144 : }
145 : }
146 : }
147 :
148 : template<>
149 : ALAMO_SINGLE_DEFINITION
150 226 : int Set::Field<Model::Solid::Affine::Isotropic>::NComp() const
151 : {
152 226 : return AMREX_SPACEDIM * AMREX_SPACEDIM + 2;
153 : }
154 : template<>
155 : ALAMO_SINGLE_DEFINITION
156 145 : std::string Set::Field<Model::Solid::Affine::Isotropic>::Name(int i) const
157 : {
158 145 : if (i == 0) return name + "_lambda";
159 130 : if (i == 1) return name + "_mu";
160 : #if AMREX_SPACEDIM==2
161 16 : if (i == 2) return name + "_F0_xx";
162 12 : if (i == 3) return name + "_F0_xy";
163 8 : if (i == 4) return name + "_F0_yx";
164 4 : if (i == 5) return name + "_F0_yy";
165 : #elif AMREX_SPACEDIM==3
166 99 : if (i == 2) return name + "_F0_xx";
167 88 : if (i == 3) return name + "_F0_xy";
168 77 : if (i == 4) return name + "_F0_xy";
169 66 : if (i == 5) return name + "_F0_yx";
170 55 : if (i == 6) return name + "_F0_yy";
171 44 : if (i == 7) return name + "_F0_yz";
172 33 : if (i == 8) return name + "_F0_zy";
173 22 : if (i == 9) return name + "_F0_zx";
174 11 : if (i == 10) return name + "_F0_zz";
175 : #endif
176 0 : return name;
177 : }
178 : template<>
179 : ALAMO_SINGLE_DEFINITION
180 51 : void Set::Field<Model::Solid::Affine::Isotropic>::Copy(int a_lev, amrex::MultiFab& a_dst, int a_dstcomp, int a_nghost) const
181 : {
182 334 : for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
183 : {
184 283 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
185 283 : if (bx.ok())
186 : {
187 283 : amrex::Array4<const Model::Solid::Affine::Isotropic> const& src = ((*this)[a_lev])->array(mfi);
188 283 : amrex::Array4<Set::Scalar> const& dst = a_dst.array(mfi);
189 1470 : for (int n = 0; n < AMREX_SPACEDIM * AMREX_SPACEDIM; n++)
190 : {
191 1187 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
192 901422 : dst(i, j, k, a_dstcomp + 0) = src(i, j, k).ddw.Lambda();
193 901422 : dst(i, j, k, a_dstcomp + 1) = src(i, j, k).ddw.Mu();
194 : #if AMREX_SPACEDIM==2
195 876672 : dst(i, j, k, a_dstcomp + 2) = src(i, j, k).F0(0, 0);
196 876672 : dst(i, j, k, a_dstcomp + 3) = src(i, j, k).F0(0, 1);
197 876672 : dst(i, j, k, a_dstcomp + 4) = src(i, j, k).F0(1, 0);
198 876672 : dst(i, j, k, a_dstcomp + 5) = src(i, j, k).F0(1, 1);
199 : #elif AMREX_SPACEDIM==3
200 24750 : dst(i, j, k, a_dstcomp + 2) = src(i, j, k).F0(0, 0);
201 24750 : dst(i, j, k, a_dstcomp + 3) = src(i, j, k).F0(0, 1);
202 24750 : dst(i, j, k, a_dstcomp + 4) = src(i, j, k).F0(0, 2);
203 24750 : dst(i, j, k, a_dstcomp + 5) = src(i, j, k).F0(1, 0);
204 24750 : dst(i, j, k, a_dstcomp + 6) = src(i, j, k).F0(1, 1);
205 24750 : dst(i, j, k, a_dstcomp + 7) = src(i, j, k).F0(1, 2);
206 24750 : dst(i, j, k, a_dstcomp + 8) = src(i, j, k).F0(2, 0);
207 24750 : dst(i, j, k, a_dstcomp + 9) = src(i, j, k).F0(2, 1);
208 24750 : dst(i, j, k, a_dstcomp + 10) = src(i, j, k).F0(2, 2);
209 : #endif
210 450711 : });
211 : }
212 : }
213 : }
214 51 : }
215 :
216 :
217 :
218 :
219 : #endif
|
