Line data Source code
1 : //
2 : // Solid models are used with the :ref:`Integrator::Mechanics` integrator, which
3 : // implements the :ref:`Solver::Nonlocal::Newton` elasticity solver.
4 : // All solid models inherit from the :code:`Model::Solid` base class, which requires
5 : // all of the necessary components to be used in a mechanics solve.
6 : // Model classes have basically two functions:
7 : //
8 : // #. Provide energy (W), stress (DW), and modulus (DDW) based on a kinematic variable
9 : // #. Evolve internal variables in time.
10 : //
11 :
12 : #ifndef MODEL_SOLID_H_
13 : #define MODEL_SOLID_H_
14 :
15 : #include <AMReX.H>
16 : #include <AMReX_REAL.H>
17 : #include <eigen3/Eigen/Core>
18 : #include <type_traits>
19 :
20 : #include "Set/Set.H"
21 : #include "Set/Matrix4.H"
22 : #include "Set/Matrix4_Diagonal.H"
23 : #include "Set/Matrix4_Full.H"
24 : #include "Set/Matrix4_Isotropic.H"
25 : #include "Set/Matrix4_Major.H"
26 : #include "Set/Matrix4_MajorMinor.H"
27 : #include "Util/Util.H"
28 :
29 : namespace Model
30 : {
31 : namespace Solid
32 : {
33 :
34 : enum KinematicVariable{gradu,epsilon,F};
35 :
36 : template<Set::Sym SYM>
37 : class Solid
38 : {
39 : public:
40 : Solid() = default;
41 :
42 : ~Solid() = default;
43 :
44 : static constexpr Set::Sym sym = SYM;
45 :
46 : Set::Scalar W(const Set::Matrix &) const
47 : {Util::Abort(INFO,"W not implemented"); return 0.0;};
48 : Set::Matrix DW(const Set::Matrix &) const
49 : {Util::Abort(INFO,"DW not implemented"); return Set::Matrix::Zero();};
50 : Set::Matrix4<AMREX_SPACEDIM,SYM> DDW(const Set::Matrix &) const
51 : {Util::Abort(INFO,"DDW not implemented"); return Set::Matrix4<AMREX_SPACEDIM,SYM>::Zero();};
52 :
53 1757488 : void Advance(Set::Scalar /*dt*/, Set::Matrix /*strain*/, Set::Matrix /*stress*/, Set::Scalar /*time*/) {}
54 :
55 : // Todo: implement for all models and make this pure virtual
56 : bool ContainsNan() {return true;};
57 :
58 :
59 : public:
60 : static const KinematicVariable kinvar = KinematicVariable::F;
61 :
62 :
63 0 : friend std::ostream& operator<<(std::ostream &out, const Solid &a)
64 : {
65 0 : a.Print(out);
66 0 : return out;
67 : }
68 0 : void Print(std::ostream &out) const
69 : {
70 0 : out << "No print function written for this model.";
71 0 : }
72 :
73 : public:
74 : template <class T>
75 24 : static int PODTest(int verbose = 0)
76 : {
77 24 : int ret = 0;
78 : if (!std::is_trivially_copyable<T>::value)
79 : {
80 36 : if (verbose) Util::Test::SubWarning("Model is not trivially copyable");
81 : }
82 : if (std::has_virtual_destructor<T>::value)
83 : {
84 : if (verbose) Util::Test::SubWarning("Model has a virtual destructor");
85 : }
86 : if (std::is_polymorphic<T>::value)
87 : {
88 : if (verbose) Util::Test::SubWarning("Model is polymorphic");
89 : }
90 : if (!std::is_standard_layout<T>::value)
91 : {
92 24 : if (verbose) Util::Test::SubWarning("Model does not have standard layout.");
93 : }
94 24 : return ret;
95 : }
96 :
97 :
98 : template <class T>
99 24 : static int ArithmeticTest(int verbose = 0)
100 : {
101 24 : const T z = T::Zero();
102 36 : if (!(z==z) && verbose)
103 : {
104 0 : Util::Message(INFO,"Zero() failed: z = \n",z); return 1;
105 : }
106 :
107 24 : const T a = T::Random();
108 36 : if (!(a==a) && verbose)
109 : {
110 0 : Util::Message(INFO,"Random() failed: a = \n",a); return 1;
111 : }
112 :
113 24 : T c = a;
114 36 : if (!(c==a) && verbose) {Util::Message(INFO); return 1;}
115 36 : c = 1.0 * a;
116 36 : if (!(c==a) && verbose) {Util::Message(INFO); return 1;}
117 36 : c = 0.0 * a;
118 36 : if (!(c==z) && verbose)
119 : {
120 0 : Util::Message(INFO,"0.0*a = \n",c);
121 0 : Util::Message(INFO,"Zero() returns \n",z);
122 0 : return 1;
123 : }
124 :
125 : // Todo: add some more arithmetic checks in here
126 :
127 24 : return 0;
128 : }
129 :
130 : template <class T>
131 24 : static int DerivativeTest1(int verbose = 0)
132 : {
133 264 : for (int iter = 0; iter < 10; iter++)
134 : {
135 240 : T model = T::Random();
136 :
137 240 : Set::Scalar dx = 1E-8, tol = 1E-6;
138 :
139 240 : Set::Matrix F = Set::Matrix::Random();
140 614 : while (F.determinant() < 0.1) F = Set::Matrix::Random(); // Ensure that F in GL(3)
141 :
142 240 : Set::Matrix dw_exact = model.DW(F);
143 240 : Set::Matrix dw_numeric = Set::Matrix::Zero();
144 840 : for (int i = 0; i < AMREX_SPACEDIM; i++)
145 2160 : for (int j = 0; j < AMREX_SPACEDIM; j++)
146 : {
147 1560 : Set::Matrix dF = Set::Matrix::Zero();
148 1560 : dF(i,j) = dx;
149 1560 : dw_numeric(i,j) = (model.W(F+dF) - model.W(F-dF)) / (2.0 * dx);
150 : }
151 240 : Set::Scalar relnorm = (dw_numeric-dw_exact).norm()/(dw_numeric.norm());
152 240 : if (relnorm > tol || std::isnan(relnorm) || std::isinf(relnorm))
153 : {
154 0 : if (verbose)
155 : {
156 0 : Util::Message(INFO,"F \n",F);
157 0 : Util::Message(INFO,"det(F) = ",F.determinant());
158 0 : Util::Message(INFO,"dw exact \n",dw_exact);
159 0 : Util::Message(INFO,"dw numeric \n",dw_numeric);
160 0 : Util::Message(INFO,"error norm \n",relnorm);
161 : }
162 0 : return 1;
163 : }
164 : }
165 24 : return 0;
166 : }
167 :
168 : template <class T>
169 24 : static int DerivativeTest2(int verbose = 0)
170 : {
171 264 : for (int iter = 0; iter < 10; iter++)
172 : {
173 240 : T model = T::Random();
174 :
175 240 : Set::Scalar dx = 1E-8, tol = 1E-4;
176 240 : Set::Matrix F = Set::Matrix::Random();
177 641 : while (F.determinant() < 0.1) F = Set::Matrix::Random(); // Ensure that F in GL(3)
178 :
179 240 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Major> ddw_exact = model.DDW(F);
180 240 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Major> ddw_numeric = Set::Matrix4<AMREX_SPACEDIM,T::sym>::Zero();
181 840 : for (int i = 0; i < AMREX_SPACEDIM; i++)
182 2160 : for (int j = 0; j < AMREX_SPACEDIM; j++)
183 5760 : for (int k = 0; k < AMREX_SPACEDIM; k++)
184 15840 : for (int l = 0; l < AMREX_SPACEDIM; l++)
185 : {
186 11640 : Set::Matrix dF = Set::Matrix::Zero();
187 11640 : dF(k,l) = dx;
188 23280 : ddw_numeric(i,j,k,l) = (model.DW(F+dF) - model.DW(F-dF))(i,j) / (2.0 * dx);
189 : }
190 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Major> error = ddw_numeric-ddw_exact;
191 240 : Set::Scalar relnorm = error.Norm()/ddw_numeric.Norm();
192 240 : if (relnorm > tol || std::isnan(relnorm) || std::isinf(relnorm))
193 : {
194 0 : if (verbose)
195 : {
196 0 : Util::Message(INFO,"F \n",F);
197 0 : Util::Message(INFO,"det(F) = ",F.determinant());
198 0 : Util::Message(INFO,"ddw exact \n",ddw_exact);
199 0 : Util::Message(INFO,"ddw numeric \n",ddw_numeric);
200 0 : Util::Message(INFO,"ddw difference \n",ddw_exact - ddw_numeric);
201 0 : Util::Message(INFO,"error norm \n",relnorm);
202 : }
203 0 : return 1;
204 : }
205 : }
206 24 : return 0;
207 : }
208 :
209 : template <class T>
210 8 : static int MaterialFrameIndifference(int verbose = 0)
211 : {
212 : if (T::kinvar != Model::Solid::KinematicVariable::F)
213 : Util::Abort(INFO,"Attempting to test material frame indifference in a non-finite-kinematics model will fail");
214 :
215 88 : for (int iter = 0; iter < 10; iter++)
216 : {
217 80 : T model = T::Random();
218 :
219 80 : Set::Scalar tol = 1E-4;
220 80 : Set::Matrix F = Set::Matrix::Random();
221 80 : while (
222 717 : F.determinant() < 0.5 || // skip if determinant negative or too small
223 80 : F.determinant() > 2 // also skip if determinant too large
224 : )
225 557 : F = Set::Matrix::Random(); // Ensure that F is reasonably well-behaved.
226 :
227 80 : Set::Scalar W_unrotated = model.W(F);
228 :
229 880 : for (int jter = 0; jter < 10; jter++)
230 : {
231 : #if AMREX_SPACEDIM == 2
232 400 : Set::Scalar theta = 2.0 * Set::Constant::Pi * Util::Random();
233 400 : Set::Matrix R;
234 400 : R(0,0) = cos(theta);
235 400 : R(0,1) = sin(theta);
236 400 : R(1,0) = -sin(theta);
237 400 : R(1,1) = cos(theta);
238 : #elif AMREX_SPACEDIM == 3
239 400 : Set::Scalar phi1 = 2.0 * Set::Constant::Pi *(Util::Random() - 0.5);
240 400 : Set::Scalar Phi = Set::Constant::Pi * Util::Random();
241 400 : Set::Scalar phi2 = 2.0 * Set::Constant::Pi *(Util::Random() - 0.5);
242 400 : Eigen::Matrix3d R;
243 400 : R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
244 400 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
245 400 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
246 : #endif
247 :
248 5600 : Util::Assert(INFO,TEST( fabs(R.determinant() - 1.0) < tol ));
249 5600 : Util::Assert(INFO,TEST( fabs((R.transpose() * R - Set::Matrix::Identity()).norm()) < tol ));
250 :
251 800 : Set::Scalar W_rotated = model.W(R*F);
252 800 : Set::Scalar error = fabs(W_rotated - W_unrotated) / fabs(W_rotated + W_unrotated + tol);
253 :
254 800 : if (error > tol)
255 : {
256 0 : if (verbose)
257 : {
258 0 : Util::Message(INFO,"F = \n",F);
259 0 : Util::Message(INFO,"R = \n",R);
260 0 : Util::Message(INFO,"W unrotated = ",W_unrotated);
261 0 : Util::Message(INFO,"W rotated = ",W_rotated);
262 : }
263 0 : return 1;
264 : }
265 : }
266 : }
267 8 : return 0;
268 : }
269 :
270 : };
271 :
272 : }
273 : }
274 :
275 :
276 : #endif
|
