1 /** @file benchmark_UniaxialTension.cpp
2
3 @brief Uniaxial Tension Test benchmark
4
5 e.g. Section 8.1. from Kiendl et al 2015
6 Kiendl, J., Hsu, M.-C., Wu, M. C. H., & Reali, A. (2015). Isogeometric Kirchhoff–Love shell formulations for general hyperelastic materials. Computer Methods in Applied Mechanics and Engineering, 291, 280–303. https://doi.org/10.1016/J.CMA.2015.03.010
7
8 This file is part of the G+Smo library.
9
10 This Source Code Form is subject to the terms of the Mozilla Public
11 License, v. 2.0. If a copy of the MPL was not distributed with this
12 file, You can obtain one at http://mozilla.org/MPL/2.0/.
13
14 Author(s): H.M. Verhelst (2019-..., TU Delft)
15 */
16
17 #include <gismo.h>
18
19 #include <gsKLShell/gsThinShellAssembler.h>
20 #include <gsKLShell/getMaterialMatrix.h>
21
22 #include <gsStructuralAnalysis/gsALMBase.h>
23 #include <gsStructuralAnalysis/gsALMLoadControl.h>
24 #include <gsStructuralAnalysis/gsALMRiks.h>
25 #include <gsStructuralAnalysis/gsALMCrisfield.h>
26
27 using namespace gismo;
28
main(int argc,char ** argv)29 int main (int argc, char** argv)
30 {
31 // Input options
32 int numElevate = 1;
33 int numRefine = 1;
34 bool plot = false;
35 bool stress = false;
36 bool quasiNewton = false;
37 int quasiNewtonInt= -1;
38 bool adaptive = false;
39 int step = 10;
40 int method = 2; // (0: Load control; 1: Riks' method; 2: Crisfield's method; 3: consistent crisfield method; 4: extended iterations)
41
42 index_t Compressibility = 0;
43 index_t material = 1;
44 bool composite = false;
45 index_t impl = 1; // 1= analytical, 2= generalized, 3= spectral
46
47 real_t relax = 1.0;
48 int result = 0;
49 index_t maxit = 20;
50 // Arc length method options
51 real_t dL = 2e0; // General arc length
52 real_t tol = 1e-7;
53 real_t tolU = 1e-7;
54 real_t tolF = 1e-4;
55
56 std::string wn("data.csv");
57
58 gsCmdLine cmd("Uniaxially stretched sheet.");
59
60 cmd.addInt("r","hRefine", "Number of dyadic h-refinement (bisection) steps to perform before solving", numRefine);
61 cmd.addInt("e","degreeElevation", "Number of degree elevation steps to perform on the Geometry's basis before solving", numElevate);
62 cmd.addInt( "M", "Material", "Material law", material );
63 cmd.addInt( "c", "Compressibility", "1: compressible, 0: incompressible", Compressibility );
64 cmd.addInt( "I", "Implementation", "Implementation: 1= analytical, 2= generalized, 3= spectral", impl );
65 cmd.addSwitch("composite", "Composite material", composite);
66
67 cmd.addInt("m","Method", "Arc length method; 1: Crisfield's method; 2: RIks' method.", method);
68 cmd.addReal("L","dL", "arc length", dL);
69 cmd.addReal("A","relaxation", "Relaxation factor for arc length method", relax);
70
71 cmd.addInt("q","QuasiNewtonInt","Use the Quasi Newton method every INT iterations",quasiNewtonInt);
72 cmd.addInt("N", "maxsteps", "Maximum number of steps", step);
73
74 cmd.addSwitch("adaptive", "Adaptive length ", adaptive);
75 cmd.addSwitch("quasi", "Use the Quasi Newton method", quasiNewton);
76 cmd.addSwitch("plot", "Plot result in ParaView format", plot);
77 cmd.addSwitch("stress", "Plot stress in ParaView format", stress);
78
79 try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
80
81 /*
82 Uniaxial tension of a square plate --- Validation settings: -L 1eX -l 1eX -M 14 -N 500 -r X -e X
83 (bottom boundary fixed in Y, left boundary fixed in X, right boundary normal load)
84 */
85 real_t mu = 1.5e6;
86 real_t thickness = 0.001;
87 real_t PoissonRatio;
88 if (!Compressibility)
89 {
90 if (material==0)
91 PoissonRatio = 0.499;
92 else
93 PoissonRatio = 0.5;
94 }
95 else
96 PoissonRatio = 0.45;
97 real_t E_modulus = 2*mu*(1+PoissonRatio);
98 real_t Density = 1e0;
99 real_t Ratio = 7.0;
100
101 gsMultiPatch<> mp,mp_def;
102 mp.addPatch( gsNurbsCreator<>::BSplineSquare(1) ); // degree
103 mp.addAutoBoundaries();
104
105 // p-refine
106 if (numElevate!=0)
107 mp.degreeElevate(numElevate);
108
109 // h-refine
110 for (int r =0; r < numRefine; ++r)
111 mp.uniformRefine();
112
113 gsInfo<<"mu = "<<E_modulus / (2 * (1 + PoissonRatio))<<"\n";
114
115 gsMultiBasis<> dbasis(mp);
116 gsInfo<<"Basis (patch 0): "<< mp.patch(0).basis() << "\n";
117 mp_def = mp;
118
119 // Boundary conditions
120 gsBoundaryConditions<> BCs;
121 BCs.setGeoMap(mp);
122 gsPointLoads<real_t> pLoads = gsPointLoads<real_t>();
123
124 BCs.addCondition(boundary::west, condition_type::dirichlet, 0, 0, false, 0 );
125 BCs.addCondition(boundary::east, condition_type::collapsed, 0, 0, false, 0 );
126 BCs.addCondition(boundary::south, condition_type::dirichlet, 0, 0, false, 1 );
127
128 BCs.setGeoMap(mp);
129
130 real_t Load = 1e0;
131 gsVector<> point(2);
132 gsVector<> load (2);
133 point<< 1.0, 0.5 ;
134 load << Load,0.0;
135 pLoads.addLoad(point, load, 0 );
136
137 std::string dirname = "ArcLengthResults";
138 std::string output = dirname + "/UniaxialTension";
139 output = "solution";
140
141 std::string commands = "mkdir -p " + dirname;
142 const char *command = commands.c_str();
143 int systemRet = system(command);
144 GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument");
145
146
147 // plot geometry
148 if (plot)
149 gsWriteParaview(mp,dirname + "/" + "mp",1000,true);
150
151 // Linear isotropic material model
152 gsFunctionExpr<> force("0","0",2);
153 gsConstantFunction<> t(thickness,2);
154 gsConstantFunction<> E(E_modulus,2);
155 gsConstantFunction<> nu(PoissonRatio,2);
156 gsConstantFunction<> rho(Density,2);
157 gsConstantFunction<> ratio(Ratio,2);
158
159 gsConstantFunction<> alpha1(1.3,2);
160 gsConstantFunction<> mu1(6.3e5/4.225e5*mu,2);
161 gsConstantFunction<> alpha2(5.0,2);
162 gsConstantFunction<> mu2(0.012e5/4.225e5*mu,2);
163 gsConstantFunction<> alpha3(-2.0,2);
164 gsConstantFunction<> mu3(-0.1e5/4.225e5*mu,2);
165
166 index_t kmax = 1;
167
168 std::vector<gsFunctionSet<> * > Gs(kmax);
169 std::vector<gsFunctionSet<> * > Ts(kmax);
170 std::vector<gsFunctionSet<> * > Phis(kmax);
171
172 gsMatrix<> Gmat = gsCompositeMatrix(E_modulus,E_modulus,0.5 * E_modulus / (1+PoissonRatio),PoissonRatio,PoissonRatio);
173 Gmat.resize(Gmat.rows()*Gmat.cols(),1);
174 gsConstantFunction<> Gfun(Gmat,3);
175 Gs[0] = &Gfun;
176
177 gsConstantFunction<> phi;
178 phi.setValue(0,3);
179
180 Phis[0] = φ
181
182 gsConstantFunction<> thicks(thickness/kmax,3);
183 Ts[0] = &thicks;
184
185 std::vector<gsFunction<>*> parameters;
186 if (material==0) // SvK & Composites
187 {
188 parameters.resize(2);
189 parameters[0] = &E;
190 parameters[1] = ν
191 }
192 else if (material==1 || material==2) // NH & NH_ext
193 {
194 parameters.resize(2);
195 parameters[0] = &E;
196 parameters[1] = ν
197 }
198 else if (material==3) // MR
199 {
200 parameters.resize(3);
201 parameters[0] = &E;
202 parameters[1] = ν
203 parameters[2] = ∶
204 }
205 else if (material==4) // OG
206 {
207 parameters.resize(8);
208 parameters[0] = &E;
209 parameters[1] = ν
210 parameters[2] = &mu1;
211 parameters[3] = &alpha1;
212 parameters[4] = &mu2;
213 parameters[5] = &alpha2;
214 parameters[6] = &mu3;
215 parameters[7] = &alpha3;
216 }
217
218 gsMaterialMatrixBase<real_t>* materialMatrix;
219
220 gsOptionList options;
221 if (material==0 && impl==1)
222 {
223 if (composite)
224 {
225 materialMatrix = new gsMaterialMatrixComposite<3,real_t>(mp,Ts,Gs,Phis);
226 }
227 else
228 {
229 options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",0);
230 options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",1);
231 materialMatrix = getMaterialMatrix<2,real_t>(mp,t,parameters,rho,options);
232 }
233 }
234 else
235 {
236 options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",material);
237 options.addSwitch("Compressibility","Compressibility: (false): Imcompressible | (true): Compressible",Compressibility);
238 options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",impl);
239 materialMatrix = getMaterialMatrix<2,real_t>(mp,t,parameters,rho,options);
240 }
241
242 gsThinShellAssemblerBase<real_t>* assembler;
243 assembler = new gsThinShellAssembler<2,real_t,false >(mp,dbasis,BCs,force,materialMatrix);
244
245 // Construct assembler object
246 assembler->setPointLoads(pLoads);
247
248 gsStopwatch stopwatch;
249 real_t time = 0.0;
250
251 typedef std::function<gsSparseMatrix<real_t> (gsVector<real_t> const &)> Jacobian_t;
252 typedef std::function<gsVector<real_t> (gsVector<real_t> const &, real_t, gsVector<real_t> const &) > ALResidual_t;
253 // Function for the Jacobian
254 Jacobian_t Jacobian = [&time,&stopwatch,&assembler,&mp_def](gsVector<real_t> const &x)
255 {
256 stopwatch.restart();
257 assembler->constructSolution(x,mp_def);
258 assembler->assembleMatrix(mp_def);
259 time += stopwatch.stop();
260
261 gsSparseMatrix<real_t> m = assembler->matrix();
262 return m;
263 };
264 // Function for the Residual
265 ALResidual_t ALResidual = [&time,&stopwatch,&assembler,&mp_def](gsVector<real_t> const &x, real_t lam, gsVector<real_t> const &force)
266 {
267 stopwatch.restart();
268 assembler->constructSolution(x,mp_def);
269 assembler->assembleVector(mp_def);
270 gsVector<real_t> Fint = -(assembler->rhs() - force);
271 gsVector<real_t> result = Fint - lam * force;
272 time += stopwatch.stop();
273 return result; // - lam * force;
274 };
275 // Assemble linear system to obtain the force vector
276 assembler->assemble();
277 gsVector<> Force = assembler->rhs();
278
279 gsALMBase<real_t> * arcLength;
280 if (method==0)
281 arcLength = new gsALMLoadControl<real_t>(Jacobian, ALResidual, Force);
282 else if (method==1)
283 arcLength = new gsALMRiks<real_t>(Jacobian, ALResidual, Force);
284 else if (method==2)
285 arcLength = new gsALMCrisfield<real_t>(Jacobian, ALResidual, Force);
286 else
287 GISMO_ERROR("Method "<<method<<" unknown");
288
289 arcLength->options().setInt("Solver",0); // CG solver
290 arcLength->options().setInt("BifurcationMethod",0); // 0: determinant, 1: eigenvalue
291 arcLength->options().setReal("Length",dL);
292 arcLength->options().setInt("AngleMethod",0); // 0: step, 1: iteration
293 arcLength->options().setSwitch("AdaptiveLength",adaptive);
294 arcLength->options().setInt("AdaptiveIterations",5);
295 arcLength->options().setReal("Scaling",0.0);
296 arcLength->options().setReal("Tol",tol);
297 arcLength->options().setReal("TolU",tolU);
298 arcLength->options().setReal("TolF",tolF);
299 arcLength->options().setInt("MaxIter",maxit);
300 arcLength->options().setSwitch("Verbose",true);
301 arcLength->options().setReal("Relaxation",relax);
302 if (quasiNewtonInt>0)
303 {
304 quasiNewton = true;
305 arcLength->options().setInt("QuasiIterations",quasiNewtonInt);
306 }
307 arcLength->options().setSwitch("Quasi",quasiNewton);
308
309
310 gsDebug<<arcLength->options();
311 arcLength->applyOptions();
312 arcLength->initialize();
313
314
315 gsParaviewCollection collection(dirname + "/" + output);
316 gsParaviewCollection Smembrane(dirname + "/" + "membrane");
317 gsParaviewCollection Sflexural(dirname + "/" + "flexural");
318 gsParaviewCollection Smembrane_p(dirname + "/" + "membrane_p");
319 gsMultiPatch<> deformation = mp;
320
321 // Make objects for previous solutions
322 real_t Lold = 0;
323 gsMatrix<> Uold = Force;
324 Uold.setZero();
325
326 gsMatrix<> solVector;
327 real_t indicator = 0.0;
328 arcLength->setIndicator(indicator); // RESET INDICATOR
329
330 gsMatrix<> lambdas(3,step);
331 gsVector<> S(step);
332 gsVector<> San(step);
333 for (index_t k=0; k<step; k++)
334 {
335 gsInfo<<"Load step "<< k<<"\n";
336 arcLength->step();
337
338 if (!(arcLength->converged()))
339 GISMO_ERROR("Loop terminated, arc length method did not converge.\n");
340
341 solVector = arcLength->solutionU();
342 Uold = solVector;
343 Lold = arcLength->solutionL();
344
345 assembler->constructSolution(solVector,mp_def);
346
347 gsMatrix<> pts(2,1);
348 pts<<0.5,0.5;
349
350 lambdas.col(k) = assembler->computePrincipalStretches(pts,mp_def,0);
351 S.at(k) = Lold / 1e-3 / lambdas(0,k) / lambdas(2,k);
352 San.at(k) = mu * (math::pow(lambdas(1,k),2)-1/lambdas(1,k));
353
354 deformation = mp_def;
355 deformation.patch(0).coefs() -= mp.patch(0).coefs();// assuming 1 patch here
356
357 gsInfo<<"Total ellapsed assembly time: "<<time<<" s\n";
358
359 if (plot)
360 {
361 gsField<> solField;
362 solField= gsField<>(mp,deformation);
363
364 std::string fileName = dirname + "/" + output + util::to_string(k);
365 gsWriteParaview<>(solField, fileName, 1000,true);
366 fileName = output + util::to_string(k) + "0";
367 collection.addTimestep(fileName,k,".vts");
368 collection.addTimestep(fileName,k,"_mesh.vtp");
369 }
370 if (stress)
371 {
372 std::string fileName;
373
374 gsField<> membraneStress, flexuralStress, membraneStress_p;
375
376 gsPiecewiseFunction<> membraneStresses;
377 assembler->constructStress(mp_def,membraneStresses,stress_type::membrane);
378 membraneStress = gsField<>(mp,membraneStresses,true);
379
380 fileName = dirname + "/" + "membrane" + util::to_string(k);
381 gsWriteParaview( membraneStress, fileName, 1000);
382 fileName = "membrane" + util::to_string(k) + "0";
383 Smembrane.addTimestep(fileName,k,".vts");
384
385 gsPiecewiseFunction<> flexuralStresses;
386 assembler->constructStress(mp_def,flexuralStresses,stress_type::flexural);
387 flexuralStress = gsField<>(mp,flexuralStresses, true);
388
389 fileName = dirname + "/" + "flexural" + util::to_string(k);
390 gsWriteParaview( flexuralStress, fileName, 1000);
391 fileName = "flexural" + util::to_string(k) + "0";
392 Sflexural.addTimestep(fileName,k,".vts");
393
394 if (impl==3)
395 {
396 gsPiecewiseFunction<> membraneStresses_p;
397 assembler->constructStress(mp_def,membraneStresses_p,stress_type::principal_stress_membrane);
398 membraneStress_p = gsField<>(mp,membraneStresses_p, true);
399
400 fileName = dirname + "/" + "membrane_p" + util::to_string(k);
401 gsWriteParaview( membraneStress_p, fileName, 1000);
402 fileName = "membrane_p" + util::to_string(k) + "0";
403 Smembrane_p.addTimestep(fileName,k,".vts");
404 }
405
406 }
407 }
408
409 gsInfo<<"Lambdas:\n"<<lambdas<<"\n";
410 gsInfo<<"S\t:\n"<<S<<"\n";
411 gsInfo<<"San\t:\n"<<San<<"\n";
412
413 if (plot)
414 {
415 collection.save();
416 Smembrane.save();
417 Sflexural.save();
418 Smembrane_p.save();
419 }
420
421 delete materialMatrix;
422 delete assembler;
423 delete arcLength;
424
425 return result;
426 }
427