1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Alessandro Tasora
13 // =============================================================================
14 //
15 // FEA (introduction to dynamics)
16 //
17 // =============================================================================
18
19 #include "chrono/physics/ChSystemSMC.h"
20 #include "chrono/solver/ChIterativeSolverLS.h"
21
22 #include "chrono/fea/ChElementSpring.h"
23 #include "chrono/fea/ChElementShellANCF_3423.h"
24 #include "chrono/fea/ChElementHexaANCF_3813.h"
25 #include "chrono/fea/ChElementBar.h"
26 #include "chrono/fea/ChElementTetraCorot_4.h"
27 #include "chrono/fea/ChElementTetraCorot_10.h"
28 #include "chrono/fea/ChElementHexaCorot_8.h"
29 #include "chrono/fea/ChElementHexaCorot_20.h"
30 #include "chrono/fea/ChMesh.h"
31 #include "chrono/fea/ChLinkPointFrame.h"
32 #include "chrono/fea/ChLinkDirFrame.h"
33
34 // Remember to use the namespace 'chrono' because all classes
35 // of Chrono::Engine belong to this namespace and its children...
36
37 using namespace chrono;
38 using namespace fea;
39
test_1()40 void test_1() {
41 GetLog() << "\n-------------------------------------------------\n";
42 GetLog() << "TEST: spring FEM dynamics, implicit integration \n\n";
43
44 // The physical system: it contains all physical objects.
45 ChSystemSMC my_system;
46
47 // Create a mesh, that is a container for groups
48 // of elements and their referenced nodes.
49 auto my_mesh = chrono_types::make_shared<ChMesh>();
50
51 // Create some nodes. These are the classical point-like
52 // nodes with x,y,z degrees of freedom, that can be used
53 // for many types of FEM elements in space.
54 auto mnodeA = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 0, 0));
55 auto mnodeB = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 1, 0));
56
57 // Default mass for FEM nodes is zero, so set point-like
58 // masses because the ChElementSpring FEM element that we
59 // are going to use won't add any mass:
60 mnodeA->SetMass(100.0);
61 mnodeB->SetMass(100.0);
62
63 // For example, set an applied force to a node:
64 mnodeB->SetForce(ChVector<>(0, 5, 0));
65
66 // For example, set an initial displacement to a node:
67 mnodeB->SetPos(mnodeB->GetX0() + ChVector<>(0, 0.1, 0));
68
69 // Remember to add nodes and elements to the mesh!
70 my_mesh->AddNode(mnodeA);
71 my_mesh->AddNode(mnodeB);
72
73 // Create some elements of 'spring-damper' type, each connecting
74 // two 3D nodes:
75 auto melementA = chrono_types::make_shared<ChElementSpring>();
76 melementA->SetNodes(mnodeA, mnodeB);
77 melementA->SetSpringK(2000);
78 melementA->SetDamperR(0);
79
80 // Remember to add elements to the mesh!
81 my_mesh->AddElement(melementA);
82
83 // Remember to add the mesh to the system!
84 my_system.Add(my_mesh);
85
86 // Create also a truss
87 auto truss = chrono_types::make_shared<ChBody>();
88 truss->SetBodyFixed(true);
89 my_system.Add(truss);
90
91 // Create a constraint between a node and the truss
92 auto constraintA = chrono_types::make_shared<ChLinkPointFrame>();
93
94 constraintA->Initialize(mnodeA, // node
95 truss); // body to be connected to
96
97 my_system.Add(constraintA);
98
99 // Set no gravity
100 // my_system.Set_G_acc(VNULL);
101
102 // Perform a dynamic time integration:
103 auto solver = chrono_types::make_shared<ChSolverMINRES>();
104 my_system.SetSolver(solver);
105 solver->SetMaxIterations(40);
106
107 my_system.SetSolverForceTolerance(1e-10);
108
109 my_system.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED);
110
111 double timestep = 0.01;
112 while (my_system.GetChTime() < 2) {
113 my_system.DoStepDynamics(timestep);
114
115 GetLog() << " t=" << my_system.GetChTime() << " nodeB pos.y()=" << mnodeB->GetPos().y() << " \n";
116 }
117 }
118
test_2()119 void test_2() {
120 GetLog() << "\n-------------------------------------------------\n";
121 GetLog() << "TEST: bar FEM dynamics, implicit integration \n\n";
122
123 // The physical system: it contains all physical objects.
124 ChSystemSMC my_system;
125
126 // Create a mesh, that is a container for groups
127 // of elements and their referenced nodes.
128 auto my_mesh = chrono_types::make_shared<ChMesh>();
129
130 // Create some nodes. These are the classical point-like
131 // nodes with x,y,z degrees of freedom, that can be used
132 // for many types of FEM elements in space.
133 auto mnodeA = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 0, 0));
134 auto mnodeB = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 1, 0));
135
136 // Set no point-like masses because mass is already in bar element:
137 mnodeA->SetMass(0.0);
138 mnodeB->SetMass(0.0);
139
140 // For example, set an applied force to a node:
141 mnodeB->SetForce(ChVector<>(0, 5, 0));
142
143 // For example, set an initial displacement to a node:
144 mnodeB->SetPos(mnodeB->GetX0() + ChVector<>(0, 0.1, 0));
145
146 // Remember to add nodes and elements to the mesh!
147 my_mesh->AddNode(mnodeA);
148 my_mesh->AddNode(mnodeB);
149
150 // Create some elements of 'bar' type, each connecting
151 // two 3D nodes:
152 auto melementA = chrono_types::make_shared<ChElementBar>();
153 melementA->SetNodes(mnodeA, mnodeB);
154 melementA->SetBarArea(0.1 * 0.02);
155 melementA->SetBarYoungModulus(0.01e9); // rubber 0.01e9, steel 200e9
156 melementA->SetBarRaleyghDamping(0.01);
157 melementA->SetBarDensity(2. * 0.1 / (melementA->GetBarArea() * 1.0));
158 // melementA->SetBarDensity(0);
159
160 // Remember to add elements to the mesh!
161 my_mesh->AddElement(melementA);
162
163 // Remember to add the mesh to the system!
164 my_system.Add(my_mesh);
165
166 // Create also a truss
167 auto truss = chrono_types::make_shared<ChBody>();
168 truss->SetBodyFixed(true);
169 my_system.Add(truss);
170
171 // Create a constraint between a node and the truss
172 auto constraintA = chrono_types::make_shared<ChLinkPointFrame>();
173
174 constraintA->Initialize(mnodeA, // node
175 truss); // body to be connected to
176
177 my_system.Add(constraintA);
178
179 // Set no gravity
180 // my_system.Set_G_acc(VNULL);
181
182 // Perform a dynamic time integration:
183
184 auto solver = chrono_types::make_shared<ChSolverMINRES>();
185 my_system.SetSolver(solver);
186 solver->SetMaxIterations(100);
187 solver->SetTolerance(1e-8);
188 solver->EnableDiagonalPreconditioner(true);
189
190 my_system.SetSolverForceTolerance(1e-10);
191
192 my_system.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT);
193
194 double timestep = 0.001;
195 while (my_system.GetChTime() < 0.2) {
196 my_system.DoStepDynamics(timestep);
197
198 GetLog() << " t=" << my_system.GetChTime() << " nodeB pos.y()=" << mnodeB->GetPos().y() << " \n";
199 }
200
201 GetLog() << " Bar mass = " << melementA->GetMass() << " restlength = " << melementA->GetRestLength() << "\n";
202 }
203
test_2b()204 void test_2b() {
205 GetLog() << "\n-------------------------------------------------\n";
206 GetLog() << "TEST: spring FEM dynamics compare to bar \n\n";
207
208 // The physical system: it contains all physical objects.
209 ChSystemSMC my_system;
210
211 // Create a mesh, that is a container for groups
212 // of elements and their referenced nodes.
213 auto my_mesh = chrono_types::make_shared<ChMesh>();
214
215 // Create some nodes. These are the classical point-like
216 // nodes with x,y,z degrees of freedom, that can be used
217 // for many types of FEM elements in space.
218 auto mnodeA = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 0, 0));
219 auto mnodeB = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 1, 0));
220
221 // Default mass for FEM nodes is zero, so set point-like
222 // masses because the ChElementSpring FEM element that we
223 // are going to use won't add any mass:
224 mnodeA->SetMass(0.1);
225 mnodeB->SetMass(0.1);
226
227 // For example, set an applied force to a node:
228 mnodeB->SetForce(ChVector<>(0, 5, 0));
229
230 // For example, set an initial displacement to a node:
231 mnodeB->SetPos(mnodeB->GetX0() + ChVector<>(0, 0.1, 0));
232
233 // Remember to add nodes and elements to the mesh!
234 my_mesh->AddNode(mnodeA);
235 my_mesh->AddNode(mnodeB);
236
237 // Create some elements of 'spring-damper' type, each connecting
238 // two 3D nodes:
239 auto melementA = chrono_types::make_shared<ChElementSpring>();
240 melementA->SetNodes(mnodeA, mnodeB);
241 melementA->SetSpringK(20000);
242 melementA->SetDamperR(200);
243
244 // Remember to add elements to the mesh!
245 my_mesh->AddElement(melementA);
246
247 // Remember to add the mesh to the system!
248 my_system.Add(my_mesh);
249
250 // Create also a truss
251 auto truss = chrono_types::make_shared<ChBody>();
252 truss->SetBodyFixed(true);
253 my_system.Add(truss);
254
255 // Create a constraint between a node and the truss
256 auto constraintA = chrono_types::make_shared<ChLinkPointFrame>();
257
258 constraintA->Initialize(mnodeA, // node
259 truss); // body to be connected to
260
261 my_system.Add(constraintA);
262
263 // Set no gravity
264 // my_system.Set_G_acc(VNULL);
265
266 // Perform a dynamic time integration:
267
268 auto solver = chrono_types::make_shared<ChSolverMINRES>();
269 my_system.SetSolver(solver);
270 solver->SetMaxIterations(200);
271 solver->SetTolerance(1e-12);
272
273 my_system.SetSolverForceTolerance(1e-10);
274
275 my_system.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT);
276
277 double timestep = 0.001;
278 while (my_system.GetChTime() < 0.2) {
279 my_system.DoStepDynamics(timestep);
280
281 GetLog() << " t=" << my_system.GetChTime() << " nodeB pos.y()=" << mnodeB->GetPos().y() << " \n";
282 }
283 }
284
test_3()285 void test_3() {
286 GetLog() << "\n-------------------------------------------------\n";
287 GetLog() << "TEST: tetrahedron FEM dynamics, implicit integration \n\n";
288
289 // The physical system: it contains all physical objects.
290 ChSystemSMC my_system;
291
292 // Create a mesh, that is a container for groups
293 // of elements and their referenced nodes.
294 auto my_mesh = chrono_types::make_shared<ChMesh>();
295
296 // Create a material, that must be assigned to each element,
297 // and set its parameters
298 auto mmaterial = chrono_types::make_shared<ChContinuumElastic>();
299 mmaterial->Set_E(0.01e9); // rubber 0.01e9, steel 200e9
300 mmaterial->Set_v(0.3);
301 mmaterial->Set_RayleighDampingK(0.01);
302 mmaterial->Set_density(1000);
303
304 // Create some nodes. These are the classical point-like
305 // nodes with x,y,z degrees of freedom, that can be used
306 // for many types of FEM elements in space.
307 auto mnode1 = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 0, 0));
308 auto mnode2 = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 0, 1));
309 auto mnode3 = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 1, 0));
310 auto mnode4 = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(1, 0, 0));
311
312 // For example, set a point-like mass at a node:
313 mnode1->SetMass(200);
314 mnode2->SetMass(200);
315 mnode3->SetMass(200);
316 mnode4->SetMass(200);
317 // For example, set an initial displacement to a node:
318 mnode3->SetPos(mnode3->GetX0() + ChVector<>(0, 0.01, 0));
319
320 // Remember to add nodes and elements to the mesh!
321 my_mesh->AddNode(mnode1);
322 my_mesh->AddNode(mnode2);
323 my_mesh->AddNode(mnode3);
324 my_mesh->AddNode(mnode4);
325
326 // Create the tetrahedron element, and assign
327 // nodes and material
328 auto melement1 = chrono_types::make_shared<ChElementTetraCorot_4>();
329 melement1->SetNodes(mnode1, mnode2, mnode3, mnode4);
330 melement1->SetMaterial(mmaterial);
331
332 // Remember to add elements to the mesh!
333 my_mesh->AddElement(melement1);
334
335 // Remember to add the mesh to the system!
336 my_system.Add(my_mesh);
337
338 // Create also a truss
339 auto truss = chrono_types::make_shared<ChBody>();
340 truss->SetBodyFixed(true);
341 my_system.Add(truss);
342
343 // Create a constraint between a node and the truss
344 auto constraint1 = chrono_types::make_shared<ChLinkPointFrame>();
345 auto constraint2 = chrono_types::make_shared<ChLinkPointFrame>();
346 auto constraint3 = chrono_types::make_shared<ChLinkPointFrame>();
347
348 constraint1->Initialize(mnode1, // node
349 truss); // body to be connected to
350
351 constraint2->Initialize(mnode2, // node
352 truss); // body to be connected to
353
354 constraint3->Initialize(mnode4, // node
355 truss); // body to be connected to
356
357 my_system.Add(constraint1);
358 my_system.Add(constraint2);
359 my_system.Add(constraint3);
360
361 // Perform a dynamic time integration:
362
363 auto solver = chrono_types::make_shared<ChSolverMINRES>();
364 my_system.SetSolver(solver);
365 solver->SetMaxIterations(40);
366 solver->SetTolerance(1e-8);
367
368 my_system.SetSolverForceTolerance(1e-10);
369
370 my_system.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED);
371
372 double timestep = 0.001;
373 while (my_system.GetChTime() < 0.1) {
374 my_system.DoStepDynamics(timestep);
375
376 GetLog() << " t =" << my_system.GetChTime() << " mnode3 pos.y()=" << mnode3->GetPos().y() << " \n";
377 }
378 }
379
test_4()380 void test_4() {
381 GetLog() << "\n-------------------------------------------------\n";
382 GetLog() << "TEST: bar FEM dynamics (2 elements), implicit integration \n\n";
383
384 // The physical system: it contains all physical objects.
385 ChSystemSMC my_system;
386
387 // Create a mesh, that is a container for groups
388 // of elements and their referenced nodes.
389 auto my_mesh = chrono_types::make_shared<ChMesh>();
390
391 // Create some nodes. These are the classical point-like
392 // nodes with x,y,z degrees of freedom, that can be used
393 // for many types of FEM elements in space.
394 auto mnodeA = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 0, 0));
395 auto mnodeB = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 1, 0));
396 auto mnodeC = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(0, 2, 0));
397
398 // Set no point-like masses because mass is already in bar element:
399 mnodeA->SetMass(0.0);
400 mnodeB->SetMass(0.0);
401 mnodeC->SetMass(0.0);
402
403 // For example, set an applied force to a node:
404 mnodeB->SetForce(ChVector<>(0, 5, 0));
405 // For example, set an applied force to a node:
406 mnodeC->SetForce(ChVector<>(0, 2, 0));
407
408 // For example, set an initial displacement to a node:
409 mnodeC->SetPos(mnodeC->GetX0() + ChVector<>(0, 0.1, 0));
410
411 // Remember to add nodes and elements to the mesh!
412 my_mesh->AddNode(mnodeA);
413 my_mesh->AddNode(mnodeB);
414 my_mesh->AddNode(mnodeC);
415
416 // Create some elements of 'bar' type, each connecting
417 // two 3D nodes:
418 auto melementA = chrono_types::make_shared<ChElementBar>();
419 melementA->SetNodes(mnodeA, mnodeB);
420 melementA->SetBarArea(0.1 * 0.02);
421 melementA->SetBarYoungModulus(0.01e9); // rubber 0.01e9, steel 200e9
422 melementA->SetBarRaleyghDamping(0.01);
423 melementA->SetBarDensity(2. * 0.1 / (melementA->GetBarArea() * 1.0));
424
425 auto melementB = chrono_types::make_shared<ChElementBar>();
426 melementB->SetNodes(mnodeB, mnodeC);
427 melementB->SetBarArea(0.1 * 0.02);
428 melementB->SetBarYoungModulus(0.01e9); // rubber 0.01e9, steel 200e9
429 melementB->SetBarRaleyghDamping(0.01);
430 melementB->SetBarDensity(2. * 0.1 / (melementB->GetBarArea() * 1.0));
431
432 // Remember to add elements to the mesh!
433 my_mesh->AddElement(melementA);
434 my_mesh->AddElement(melementB);
435
436 // Remember to add the mesh to the system!
437 my_system.Add(my_mesh);
438
439 // Create also a truss
440 auto truss = chrono_types::make_shared<ChBody>();
441 truss->SetBodyFixed(true);
442 my_system.Add(truss);
443
444 // Create a constraint between a node and the truss
445 auto constraintA = chrono_types::make_shared<ChLinkPointFrame>();
446
447 constraintA->Initialize(mnodeA, // node
448 truss); // body to be connected to
449
450 my_system.Add(constraintA);
451
452 // Set no gravity
453 // my_system.Set_G_acc(VNULL);
454
455 // Perform a dynamic time integration:
456
457 auto solver = chrono_types::make_shared<ChSolverMINRES>();
458 my_system.SetSolver(solver);
459 solver->SetMaxIterations(100);
460 solver->SetTolerance(1e-8);
461 solver->EnableDiagonalPreconditioner(true);
462
463 my_system.SetSolverForceTolerance(1e-10);
464
465 my_system.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT);
466
467 double timestep = 0.001;
468 while (my_system.GetChTime() < 0.2) {
469 my_system.DoStepDynamics(timestep);
470
471 GetLog() << " t=" << my_system.GetChTime() << " nodeB pos.y()=" << mnodeB->GetPos().y()
472 << " nodeC pos.y()=" << mnodeC->GetPos().y() << " \n";
473 }
474 }
475
main(int argc,char * argv[])476 int main(int argc, char* argv[]) {
477 GetLog() << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << "\n\n";
478
479 // test_1();
480 // test_2();
481 test_3();
482 // test_4();
483
484 return 0;
485 }
486