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 0 : }
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 0 : }
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 :
82 30 : std::pair<std::string, Set::Scalar> moduli[2];
83 :
84 48 : pp.forbid("lame","Use 'lambda' instead for lame constant");
85 48 : pp.forbid("shear","Use 'mu' instead for shear modulus");
86 42 : pp.forbid("bulk","Use 'kappa' instead for bulk modulus");
87 :
88 : // Specify exactly two of: lame constant \(\lambda\),
89 : // shear modulus \(\mu\), Young's modulus \(E\), Poisson's ratio \(\nu\),
90 : // bulk modulus \(\kappa\).
91 : // \(\mu\) and \(\lambda\) are how the final values are stored.
92 12 : pp.query_exactly<2>({"lambda","mu","E","nu","kappa"}, moduli);
93 :
94 6 : if (moduli[0].first == "lambda" && moduli[1].first == "mu")
95 : {
96 0 : lambda = moduli[0].second;
97 0 : mu = moduli[1].second;
98 : }
99 6 : else if (moduli[0].first == "lambda" && moduli[1].first == "E")
100 : {
101 0 : lambda = moduli[0].second;
102 0 : Set::Scalar E = moduli[1].second;
103 0 : mu = (E - 3.0 * lambda + sqrt(E * E + 9.0 * lambda * lambda + 2.0 * E * lambda)) / 4.0;
104 : }
105 6 : else if (moduli[0].first == "lambda" && moduli[1].first == "nu")
106 : {
107 0 : lambda = moduli[0].second;
108 0 : Set::Scalar nu = moduli[1].second;
109 0 : mu = lambda * (1.0 - 2.0 * nu) / 2.0 / nu;
110 : }
111 6 : else if (moduli[0].first == "mu" && moduli[1].first == "E")
112 : {
113 0 : mu = moduli[0].second;
114 0 : Set::Scalar E = moduli[1].second;
115 0 : lambda = mu * (E - 2.0 * mu) / (3.0 * mu - E);
116 : }
117 6 : else if (moduli[0].first == "mu" && moduli[1].first == "nu")
118 : {
119 0 : mu = moduli[0].second;
120 0 : Set::Scalar nu = moduli[1].second;
121 0 : lambda = 2.0 * mu * nu / (1.0 - 2.0 * nu);
122 : }
123 6 : else if (moduli[0].first == "mu" && moduli[1].first == "K")
124 : {
125 0 : mu = moduli[0].second;
126 0 : Set::Scalar K = moduli[1].second;
127 0 : lambda = K - (2.0 * mu / 3.0);
128 : }
129 6 : else if (moduli[0].first == "E" && moduli[1].first == "nu")
130 : {
131 6 : Set::Scalar E = moduli[0].second, nu = moduli[1].second;
132 6 : lambda = E * nu / (1.0 + nu) / (1.0 - 2.0 * nu);
133 6 : mu = E / 2.0 / (1.0 + nu);
134 : }
135 : else
136 : {
137 0 : Util::Exception(INFO,"Haven't implemented ",moduli[0].first," and ",moduli[1].first," (sorry!)");
138 : }
139 :
140 12 : if (pp.contains("F0"))
141 : {
142 30 : pp_queryarr("F0", F0); // Eigendeformation gradient
143 : }
144 6 : value.Define(mu, lambda, F0);
145 24 : }
146 : #define OP_CLASS Isotropic
147 : #define OP_VARS X(ddw) X(F0)
148 : #include "Model/Solid/InClassOperators.H"
149 : };
150 : #include "Model/Solid/ExtClassOperators.H"
151 : }
152 : }
153 : }
154 :
155 : template<>
156 : ALAMO_SINGLE_DEFINITION
157 226 : int Set::Field<Model::Solid::Affine::Isotropic>::NComp() const
158 : {
159 226 : return AMREX_SPACEDIM * AMREX_SPACEDIM + 2;
160 : }
161 : template<>
162 : ALAMO_SINGLE_DEFINITION
163 145 : std::string Set::Field<Model::Solid::Affine::Isotropic>::Name(int i) const
164 : {
165 145 : if (i == 0) return name + "_lambda";
166 130 : if (i == 1) return name + "_mu";
167 : #if AMREX_SPACEDIM==2
168 16 : if (i == 2) return name + "_F0_xx";
169 12 : if (i == 3) return name + "_F0_xy";
170 8 : if (i == 4) return name + "_F0_yx";
171 4 : if (i == 5) return name + "_F0_yy";
172 : #elif AMREX_SPACEDIM==3
173 99 : if (i == 2) return name + "_F0_xx";
174 88 : if (i == 3) return name + "_F0_xy";
175 77 : if (i == 4) return name + "_F0_xy";
176 66 : if (i == 5) return name + "_F0_yx";
177 55 : if (i == 6) return name + "_F0_yy";
178 44 : if (i == 7) return name + "_F0_yz";
179 33 : if (i == 8) return name + "_F0_zy";
180 22 : if (i == 9) return name + "_F0_zx";
181 11 : if (i == 10) return name + "_F0_zz";
182 : #endif
183 0 : return name;
184 : }
185 : template<>
186 : ALAMO_SINGLE_DEFINITION
187 51 : void Set::Field<Model::Solid::Affine::Isotropic>::Copy(int a_lev, amrex::MultiFab& a_dst, int a_dstcomp, int a_nghost) const
188 : {
189 334 : for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
190 : {
191 283 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
192 283 : if (bx.ok())
193 : {
194 283 : amrex::Array4<const Model::Solid::Affine::Isotropic> const& src = ((*this)[a_lev])->array(mfi);
195 283 : amrex::Array4<Set::Scalar> const& dst = a_dst.array(mfi);
196 1470 : for (int n = 0; n < AMREX_SPACEDIM * AMREX_SPACEDIM; n++)
197 : {
198 1187 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
199 901422 : dst(i, j, k, a_dstcomp + 0) = src(i, j, k).ddw.Lambda();
200 901422 : dst(i, j, k, a_dstcomp + 1) = src(i, j, k).ddw.Mu();
201 : #if AMREX_SPACEDIM==2
202 876672 : dst(i, j, k, a_dstcomp + 2) = src(i, j, k).F0(0, 0);
203 876672 : dst(i, j, k, a_dstcomp + 3) = src(i, j, k).F0(0, 1);
204 876672 : dst(i, j, k, a_dstcomp + 4) = src(i, j, k).F0(1, 0);
205 876672 : dst(i, j, k, a_dstcomp + 5) = src(i, j, k).F0(1, 1);
206 : #elif AMREX_SPACEDIM==3
207 24750 : dst(i, j, k, a_dstcomp + 2) = src(i, j, k).F0(0, 0);
208 24750 : dst(i, j, k, a_dstcomp + 3) = src(i, j, k).F0(0, 1);
209 24750 : dst(i, j, k, a_dstcomp + 4) = src(i, j, k).F0(0, 2);
210 24750 : dst(i, j, k, a_dstcomp + 5) = src(i, j, k).F0(1, 0);
211 24750 : dst(i, j, k, a_dstcomp + 6) = src(i, j, k).F0(1, 1);
212 24750 : dst(i, j, k, a_dstcomp + 7) = src(i, j, k).F0(1, 2);
213 24750 : dst(i, j, k, a_dstcomp + 8) = src(i, j, k).F0(2, 0);
214 24750 : dst(i, j, k, a_dstcomp + 9) = src(i, j, k).F0(2, 1);
215 24750 : dst(i, j, k, a_dstcomp + 10) = src(i, j, k).F0(2, 2);
216 : #endif
217 450711 : });
218 : }
219 : }
220 51 : }
221 51 : }
222 :
223 :
224 :
225 :
226 : #endif
|
