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 :
19 : #include "Set/Set.H"
20 :
21 : namespace Model
22 : {
23 : namespace Solid
24 : {
25 :
26 : enum KinematicVariable{gradu,epsilon,F};
27 :
28 : template<Set::Sym SYM>
29 : class Solid
30 : {
31 : public:
32 497603 : Solid() {} ;
33 :
34 697074 : virtual ~Solid() {};
35 :
36 : static constexpr Set::Sym sym = SYM;
37 :
38 0 : virtual Set::Scalar W(const Set::Matrix &) const
39 0 : {Util::Abort(INFO,"W not implemented"); return 0.0;};
40 0 : virtual Set::Matrix DW(const Set::Matrix &) const
41 0 : {Util::Abort(INFO,"DW not implemented"); return Set::Matrix::Zero();};
42 0 : virtual Set::Matrix4<AMREX_SPACEDIM,SYM> DDW(const Set::Matrix &) const
43 0 : {Util::Abort(INFO,"DDW not implemented"); return Set::Matrix4<AMREX_SPACEDIM,SYM>::Zero();};
44 :
45 1757488 : virtual void Advance(Set::Scalar /*dt*/, Set::Matrix /*strain*/, Set::Matrix /*stress*/, Set::Scalar /*time*/) {}
46 :
47 : // Todo: implement for all models and make this pure virtual
48 0 : virtual bool ContainsNan() {return true;};
49 :
50 :
51 : public:
52 : static const KinematicVariable kinvar = KinematicVariable::F;
53 :
54 :
55 0 : friend std::ostream& operator<<(std::ostream &out, const Solid &a)
56 : {
57 0 : a.Print(out);
58 0 : return out;
59 : }
60 0 : virtual void Print(std::ostream &out) const
61 : {
62 0 : out << "No print function written for this model.";
63 0 : }
64 :
65 : public:
66 : template <class T>
67 22 : static int ArithmeticTest(int verbose = 0)
68 : {
69 22 : const T z = T::Zero();
70 44 : if (!(z==z) && verbose)
71 : {
72 0 : Util::Message(INFO,"Zero() failed: z = \n",z); return 1;
73 : }
74 :
75 22 : const T a = T::Random();
76 44 : if (!(a==a) && verbose)
77 : {
78 0 : Util::Message(INFO,"Random() failed: a = \n",a); return 1;
79 : }
80 :
81 22 : const T b = T::Random();
82 :
83 22 : T c = a;
84 44 : if (!(c==a) && verbose) {Util::Message(INFO); return 1;}
85 44 : c = 1.0 * a;
86 44 : if (!(c==a) && verbose) {Util::Message(INFO); return 1;}
87 44 : c = 0.0 * a;
88 44 : if (!(c==z) && verbose)
89 : {
90 0 : Util::Message(INFO,"0.0*a = \n",c);
91 0 : Util::Message(INFO,"Zero() returns \n",z);
92 0 : return 1;
93 : }
94 :
95 : // Todo: add some more arithmetic checks in here
96 :
97 22 : return 0;
98 22 : }
99 :
100 : template <class T>
101 22 : static int DerivativeTest1(int verbose = 0)
102 : {
103 462 : for (int iter = 0; iter < 10; iter++)
104 : {
105 220 : T model = T::Random();
106 :
107 220 : Set::Scalar dx = 1E-8, tol = 1E-6;
108 :
109 220 : Set::Matrix F = Set::Matrix::Random();
110 619 : while (F.determinant() < 0.1) F = Set::Matrix::Random(); // Ensure that F in GL(3)
111 :
112 220 : Set::Matrix dw_exact = model.DW(F);
113 220 : Set::Matrix dw_numeric = Set::Matrix::Zero();
114 770 : for (int i = 0; i < AMREX_SPACEDIM; i++)
115 1980 : for (int j = 0; j < AMREX_SPACEDIM; j++)
116 : {
117 1430 : Set::Matrix dF = Set::Matrix::Zero();
118 1430 : dF(i,j) = dx;
119 1430 : dw_numeric(i,j) = (model.W(F+dF) - model.W(F-dF)) / (2.0 * dx);
120 : }
121 220 : Set::Scalar relnorm = (dw_numeric-dw_exact).norm()/(dw_numeric.norm());
122 220 : if (relnorm > tol || std::isnan(relnorm) || std::isinf(relnorm))
123 : {
124 0 : if (verbose)
125 : {
126 0 : Util::Message(INFO,"F \n",F);
127 0 : Util::Message(INFO,"det(F) = ",F.determinant());
128 0 : Util::Message(INFO,"dw exact \n",dw_exact);
129 0 : Util::Message(INFO,"dw numeric \n",dw_numeric);
130 0 : Util::Message(INFO,"error norm \n",relnorm);
131 : }
132 0 : return 1;
133 : }
134 : }
135 22 : return 0;
136 : }
137 :
138 : template <class T>
139 22 : static int DerivativeTest2(int verbose = 0)
140 : {
141 462 : for (int iter = 0; iter < 10; iter++)
142 : {
143 220 : T model = T::Random();
144 :
145 220 : Set::Scalar dx = 1E-8, tol = 1E-4;
146 220 : Set::Matrix F = Set::Matrix::Random();
147 551 : while (F.determinant() < 0.1) F = Set::Matrix::Random(); // Ensure that F in GL(3)
148 :
149 220 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Major> ddw_exact = model.DDW(F);
150 220 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Major> ddw_numeric = Set::Matrix4<AMREX_SPACEDIM,T::sym>::Zero();
151 770 : for (int i = 0; i < AMREX_SPACEDIM; i++)
152 1980 : for (int j = 0; j < AMREX_SPACEDIM; j++)
153 5280 : for (int k = 0; k < AMREX_SPACEDIM; k++)
154 14520 : for (int l = 0; l < AMREX_SPACEDIM; l++)
155 : {
156 10670 : Set::Matrix dF = Set::Matrix::Zero();
157 10670 : dF(k,l) = dx;
158 21340 : ddw_numeric(i,j,k,l) = (model.DW(F+dF) - model.DW(F-dF))(i,j) / (2.0 * dx);
159 : }
160 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Major> error = ddw_numeric-ddw_exact;
161 220 : Set::Scalar relnorm = error.Norm()/ddw_numeric.Norm();
162 220 : if (relnorm > tol || std::isnan(relnorm) || std::isinf(relnorm))
163 : {
164 0 : if (verbose)
165 : {
166 0 : Util::Message(INFO,"F \n",F);
167 0 : Util::Message(INFO,"det(F) = ",F.determinant());
168 0 : Util::Message(INFO,"ddw exact \n",ddw_exact);
169 0 : Util::Message(INFO,"ddw numeric \n",ddw_numeric);
170 0 : Util::Message(INFO,"ddw difference \n",ddw_exact - ddw_numeric);
171 0 : Util::Message(INFO,"error norm \n",relnorm);
172 : }
173 0 : return 1;
174 : }
175 : }
176 22 : return 0;
177 : }
178 :
179 : template <class T>
180 8 : static int MaterialFrameIndifference(int verbose = 0)
181 : {
182 : if (T::kinvar != Model::Solid::KinematicVariable::F)
183 : Util::Abort(INFO,"Attempting to test material frame indifference in a non-finite-kinematics model will fail");
184 :
185 168 : for (int iter = 0; iter < 10; iter++)
186 : {
187 80 : T model = T::Random();
188 :
189 80 : Set::Scalar tol = 1E-4;
190 80 : Set::Matrix F = Set::Matrix::Random();
191 80 : while (
192 641 : F.determinant() < 0.5 || // skip if determinant negative or too small
193 80 : F.determinant() > 2 // also skip if determinant too large
194 : )
195 481 : F = Set::Matrix::Random(); // Ensure that F is reasonably well-behaved.
196 :
197 80 : Set::Scalar W_unrotated = model.W(F);
198 :
199 880 : for (int jter = 0; jter < 10; jter++)
200 : {
201 : #if AMREX_SPACEDIM == 2
202 400 : Set::Scalar theta = 2.0 * Set::Constant::Pi * Util::Random();
203 400 : Set::Matrix R;
204 400 : R(0,0) = cos(theta);
205 400 : R(0,1) = sin(theta);
206 400 : R(1,0) = -sin(theta);
207 400 : R(1,1) = cos(theta);
208 : #elif AMREX_SPACEDIM == 3
209 400 : Set::Scalar phi1 = 2.0 * Set::Constant::Pi *(Util::Random() - 0.5);
210 400 : Set::Scalar Phi = Set::Constant::Pi * Util::Random();
211 400 : Set::Scalar phi2 = 2.0 * Set::Constant::Pi *(Util::Random() - 0.5);
212 400 : Eigen::Matrix3d R;
213 400 : R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
214 400 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
215 400 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
216 : #endif
217 :
218 5600 : Util::Assert(INFO,TEST( fabs(R.determinant() - 1.0) < tol ));
219 5600 : Util::Assert(INFO,TEST( fabs((R.transpose() * R - Set::Matrix::Identity()).norm()) < tol ));
220 :
221 800 : Set::Scalar W_rotated = model.W(R*F);
222 800 : Set::Scalar error = fabs(W_rotated - W_unrotated) / fabs(W_rotated + W_unrotated + tol);
223 :
224 800 : if (error > tol)
225 : {
226 0 : if (verbose)
227 : {
228 0 : Util::Message(INFO,"F = \n",F);
229 0 : Util::Message(INFO,"R = \n",R);
230 0 : Util::Message(INFO,"W unrotated = ",W_unrotated);
231 0 : Util::Message(INFO,"W rotated = ",W_rotated);
232 : }
233 0 : return 1;
234 : }
235 : }
236 : }
237 8 : return 0;
238 : }
239 :
240 : };
241 :
242 : }
243 : }
244 :
245 :
246 : #endif
|
