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