Line data Source code
1 : //
2 : // This model implements an isotropic linear elastic material.
3 : // See `this link <https://en.wikipedia.org/wiki/Linear_elasticity#(An)isotropic_(in)homogeneous_media>`_
4 : // for more information about the theory.
5 : //
6 : // Free energy for a linear material is defined as
7 : //
8 : // .. math::
9 : //
10 : // W(\nabla\mathbf{u}) =
11 : // \frac{1}{2}\nabla\mathbf{u}\cdot\mathbb{C}\,\nabla\mathbf{u}
12 : //
13 : // For an isotropic material, stress and strain are related through
14 : //
15 : // .. math::
16 : //
17 : // \mathbb{C}_{ijkl} = \lambda \delta_{ij}\varepsilon_{kk} + 2\mu\varepsilon_{ij}
18 : //
19 : // where :math:`\lambda` and :math:`\mu` are the Lame constant and shear modulus, respectively.
20 : // Users can specify these either through (:code:`lame` and :code:`shear`)
21 : // OR (:code:`lambda` and :code:`mu`) OR (:code:`E` and :code:`nu`).
22 : //
23 : // Class methods:
24 : //
25 : // #. :code:`Isotropic()`:
26 : // Basic constructor. Does nothing, and leaves all values initiated as NAN.
27 : // #. :code:`Isotropic(Solid<Set::Sym::Isotropic> base)`
28 : // Basic constructor. Does nothing gut allows for inheritance.
29 : // #. :code:`Isotropic(Set::Scalar a_mu, Set::Scalar a_lambda)`
30 : // BAD old-fashioned constructor. Do not use!
31 : // #. :code:`~Isotropic()`
32 : // Simple destructor. Don't need to change it.
33 : // #. :code:`void Define(Set::Scalar a_mu, Set::Scalar a_lambda)`
34 : // BAD old-fashioned way of doing things. Use :code:`Parse` instead.
35 : // #. :code:`Set::Scalar W(const Set::Matrix & gradu) const override`
36 : // Returns elastic free energy density
37 : // #. :code:`Set::Matrix DW(const Set::Matrix & gradu) const override`
38 : // Returns first derivative of free energy, the stress tensor
39 : // #. :code:`Set::Matrix4<...> DDW(const Set::Matrix & ) const override`
40 : // Returns second derivative of free energy, the modulus tensor
41 : // #. :code:`virtual void Print(std::ostream &out) const override`
42 : // Prints the modulus tensor object to output stream (usually the terminal)
43 : // #. :code:`static Isotropic Random()`
44 : // Static method that generates a random yet acceptable model.
45 : // #. :code:`static Isotropic Zero()`
46 : // Static method that generates a "zero" element (so that there is no effect under addition)
47 : // #. :code:`static void Parse(Isotropic & value, IO::ParmParse & pp)`
48 : // Parser where all the IO occurs
49 : //
50 : //
51 :
52 : #ifndef MODEL_SOLID_LINEAR_ISOTROPIC_H_
53 : #define MODEL_SOLID_LINEAR_ISOTROPIC_H_
54 :
55 : #include "Model/Solid/Solid.H"
56 : #include "IO/ParmParse.H"
57 : #include "Set/Matrix4_Isotropic.H"
58 :
59 : namespace Model
60 : {
61 : namespace Solid
62 : {
63 : namespace Linear
64 : {
65 : class Isotropic : public Solid<Set::Sym::Isotropic>
66 : {
67 : public:
68 :
69 207892 : Isotropic() = default;
70 : Isotropic(Solid<Set::Sym::Isotropic> base) : Solid<Set::Sym::Isotropic>(base) {};
71 : Isotropic(Set::Scalar a_mu, Set::Scalar a_lambda)
72 : {
73 : Define(a_mu,a_lambda);
74 : };
75 : ~Isotropic() = default;
76 :
77 49904 : void Define(Set::Scalar a_mu, Set::Scalar a_lambda)
78 : {
79 49904 : ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Isotropic>(a_lambda,a_mu);
80 49904 : }
81 :
82 260 : Set::Scalar W(const Set::Matrix & gradu) const
83 : {
84 520 : return ( 0.5 * gradu.transpose() * (ddw*gradu) ).trace();
85 : }
86 134932 : Set::Matrix DW(const Set::Matrix & gradu) const
87 : {
88 269864 : return ddw*gradu;
89 : }
90 66506 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Isotropic> DDW(const Set::Matrix & /*gradu*/) const
91 : {
92 66506 : return ddw;
93 : }
94 : void Print(std::ostream &out) const
95 : {
96 : out << ddw;
97 : }
98 :
99 : public:
100 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Isotropic> ddw;
101 : static const KinematicVariable kinvar = KinematicVariable::gradu;
102 :
103 : public:
104 42 : static Isotropic Random()
105 : {
106 42 : Isotropic ret;
107 42 : ret.Define(Util::Random(),Util::Random());
108 42 : return ret;
109 : }
110 49855 : static Isotropic Zero()
111 : {
112 49855 : Isotropic ret;
113 49855 : ret.Define(0.,0.);
114 49855 : return ret;
115 : }
116 7 : static void Parse(Isotropic & value, IO::ParmParse & pp)
117 : {
118 7 : Set::Scalar mu = NAN, lambda = NAN;
119 7 : bool planestress = false;
120 42 : std::pair<std::string, Set::Scalar> moduli[2];
121 :
122 28 : pp.forbid("lame","Use 'lambda' instead for lame constant");
123 28 : pp.forbid("shear","Use 'mu' instead for shear modulus");
124 21 : pp.forbid("bulk","Use 'kappa' instead for bulk modulus");
125 :
126 : // Specify exactly two of: lame constant \(\lambda\),
127 : // shear modulus \(\mu\), Young's modulus \(E\), Poisson's ratio \(\nu\),
128 : // bulk modulus \(\kappa\).
129 : // \(\mu\) and \(\lambda\) are how the final values are stored.
130 21 : pp.query_exactly<2>({"lambda","mu","E","nu","kappa"}, moduli);
131 :
132 7 : if (moduli[0].first == "lambda" && moduli[1].first == "mu")
133 : {
134 0 : lambda = moduli[0].second;
135 0 : mu = moduli[1].second;
136 : }
137 7 : else if (moduli[0].first == "lambda" && moduli[1].first == "E")
138 : {
139 0 : lambda = moduli[0].second;
140 0 : Set::Scalar E = moduli[1].second;
141 0 : mu = (E - 3.0 * lambda + sqrt(E * E + 9.0 * lambda * lambda + 2.0 * E * lambda)) / 4.0;
142 : }
143 7 : else if (moduli[0].first == "lambda" && moduli[1].first == "nu")
144 : {
145 0 : lambda = moduli[0].second;
146 0 : Set::Scalar nu = moduli[1].second;
147 0 : mu = lambda * (1.0 - 2.0 * nu) / 2.0 / nu;
148 : }
149 7 : else if (moduli[0].first == "mu" && moduli[1].first == "E")
150 : {
151 0 : mu = moduli[0].second;
152 0 : Set::Scalar E = moduli[1].second;
153 0 : lambda = mu * (E - 2.0 * mu) / (3.0 * mu - E);
154 : }
155 7 : else if (moduli[0].first == "mu" && moduli[1].first == "nu")
156 : {
157 0 : mu = moduli[0].second;
158 0 : Set::Scalar nu = moduli[1].second;
159 0 : lambda = 2.0 * mu * nu / (1.0 - 2.0 * nu);
160 : }
161 7 : else if (moduli[0].first == "mu" && moduli[1].first == "K")
162 : {
163 0 : mu = moduli[0].second;
164 0 : Set::Scalar K = moduli[1].second;
165 0 : lambda = K - (2.0 * mu / 3.0);
166 : }
167 7 : else if (moduli[0].first == "E" && moduli[1].first == "nu")
168 : {
169 7 : Set::Scalar E = moduli[0].second, nu = moduli[1].second;
170 7 : lambda = E * nu / (1.0 + nu) / (1.0 - 2.0 * nu);
171 7 : mu = E / 2.0 / (1.0 + nu);
172 : }
173 : else
174 : {
175 0 : Util::Exception(INFO,"Haven't implemented ",moduli[0].first," and ",moduli[1].first," (sorry!)");
176 : }
177 :
178 : // Whether or not to use the
179 : // `plane stress <https://en.wikipedia.org/wiki/Plane_stress>`_
180 : // approximation.
181 8 : pp.query_default("planestress",planestress, false);
182 :
183 6 : if (AMREX_SPACEDIM==2 && planestress)
184 5 : value.Define(mu,lambda*(1.0 - lambda/(2.*mu + lambda)));
185 : else
186 2 : value.Define(mu,lambda);
187 28 : }
188 :
189 : #define OP_CLASS Isotropic
190 : #define OP_VARS X(ddw)
191 : #include "Model/Solid/InClassOperators.H"
192 : };
193 : #include "Model/Solid/ExtClassOperators.H"
194 :
195 :
196 :
197 : }
198 : }
199 : }
200 :
201 : template<>
202 : ALAMO_SINGLE_DEFINITION
203 184 : int Set::Field<Model::Solid::Linear::Isotropic>::NComp() const
204 : {
205 184 : return 2;
206 : }
207 : template<>
208 : ALAMO_SINGLE_DEFINITION
209 48 : std::string Set::Field<Model::Solid::Linear::Isotropic>::Name(int i) const
210 : {
211 48 : if (i == 0) return name + "_lambda";
212 24 : if (i == 1) return name + "_mu";
213 0 : return name;
214 : }
215 : template<>
216 : ALAMO_SINGLE_DEFINITION
217 88 : void Set::Field<Model::Solid::Linear::Isotropic>::Copy(int a_lev, amrex::MultiFab& a_dst, int a_dstcomp, int a_nghost) const
218 : {
219 461 : for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
220 : {
221 373 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
222 373 : if (bx.ok())
223 : {
224 373 : amrex::Array4<const Model::Solid::Linear::Isotropic> const& src = ((*this)[a_lev])->array(mfi);
225 373 : amrex::Array4<Set::Scalar> const& dst = a_dst.array(mfi);
226 1920 : for (int n = 0; n < AMREX_SPACEDIM * AMREX_SPACEDIM; n++)
227 : {
228 1547 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
229 1241054 : dst(i, j, k, a_dstcomp + 0) = src(i, j, k).ddw.Lambda();
230 1241054 : dst(i, j, k, a_dstcomp + 1) = src(i, j, k).ddw.Mu();
231 620527 : });
232 : }
233 : }
234 88 : }
235 88 : }
236 :
237 : #endif
|
