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 for 3D beams
16 //
17 // =============================================================================
18
19 #include "chrono/physics/ChSystemSMC.h"
20 #include "chrono/physics/ChLinkMate.h"
21 #include "chrono/physics/ChBodyEasy.h"
22 #include "chrono/timestepper/ChTimestepper.h"
23 #include "chrono/solver/ChIterativeSolverLS.h"
24
25 #include "chrono/fea/ChElementBeamEuler.h"
26 #include "chrono/fea/ChBuilderBeam.h"
27 #include "chrono/fea/ChMesh.h"
28 #include "chrono/fea/ChVisualizationFEAmesh.h"
29 #include "chrono/fea/ChLinkPointFrame.h"
30 #include "chrono/fea/ChLinkDirFrame.h"
31
32 #include "chrono_irrlicht/ChIrrApp.h"
33
34 using namespace chrono;
35 using namespace chrono::fea;
36 using namespace chrono::irrlicht;
37
38 using namespace irr;
39
main(int argc,char * argv[])40 int main(int argc, char* argv[]) {
41 GetLog() << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << "\n\n";
42
43 // Create a Chrono::Engine physical system
44 ChSystemSMC my_system;
45
46 // Create the Irrlicht visualization (open the Irrlicht device,
47 // bind a simple user interface, etc. etc.)
48 ChIrrApp application(&my_system, L"Beams (SPACE for dynamics, F10 / F11 statics)", core::dimension2d<u32>(800, 600));
49
50 // Easy shortcuts to add camera, lights, logo and sky in Irrlicht scene:
51 application.AddTypicalLogo();
52 application.AddTypicalSky();
53 application.AddTypicalLights();
54 application.AddTypicalCamera(core::vector3df(-0.1f, 0.2f, -0.2f));
55
56 // Create a mesh, that is a container for groups
57 // of elements and their referenced nodes.
58 auto my_mesh = chrono_types::make_shared<ChMesh>();
59
60 // Create a section, i.e. thickness and material properties
61 // for beams. This will be shared among some beams.
62
63 auto msection = chrono_types::make_shared<ChBeamSectionEulerAdvanced>();
64
65 double beam_wy = 0.012;
66 double beam_wz = 0.025;
67 msection->SetAsRectangularSection(beam_wy, beam_wz);
68 msection->SetYoungModulus(0.01e9);
69 msection->SetGshearModulus(0.01e9 * 0.3);
70 msection->SetBeamRaleyghDamping(0.000);
71 // msection->SetCentroid(0,0.02);
72 // msection->SetShearCenter(0,0.1);
73 // msection->SetSectionRotation(45*CH_C_RAD_TO_DEG);
74
75 //
76 // Add some EULER-BERNOULLI BEAMS:
77 //
78
79 double beam_L = 0.1;
80
81 auto hnode1 = chrono_types::make_shared<ChNodeFEAxyzrot>(ChFrame<>(ChVector<>(0, 0, 0)));
82 auto hnode2 = chrono_types::make_shared<ChNodeFEAxyzrot>(ChFrame<>(ChVector<>(beam_L, 0, 0)));
83 auto hnode3 = chrono_types::make_shared<ChNodeFEAxyzrot>(ChFrame<>(ChVector<>(beam_L * 2, 0, 0)));
84
85 my_mesh->AddNode(hnode1);
86 my_mesh->AddNode(hnode2);
87 my_mesh->AddNode(hnode3);
88
89 auto belement1 = chrono_types::make_shared<ChElementBeamEuler>();
90
91 belement1->SetNodes(hnode1, hnode2);
92 belement1->SetSection(msection);
93
94 my_mesh->AddElement(belement1);
95
96 auto belement2 = chrono_types::make_shared<ChElementBeamEuler>();
97
98 belement2->SetNodes(hnode2, hnode3);
99 belement2->SetSection(msection);
100
101 my_mesh->AddElement(belement2);
102
103 // Apply a force or a torque to a node:
104 hnode2->SetForce(ChVector<>(4, 2, 0));
105 // hnode3->SetTorque( ChVector<>(0, -0.04, 0));
106
107 // Fix a node to ground:
108 // hnode1->SetFixed(true);
109 auto mtruss = chrono_types::make_shared<ChBody>();
110 mtruss->SetBodyFixed(true);
111 my_system.Add(mtruss);
112
113 auto constr_bc = chrono_types::make_shared<ChLinkMateGeneric>();
114 constr_bc->Initialize(hnode3, mtruss, false, hnode3->Frame(), hnode3->Frame());
115 my_system.Add(constr_bc);
116 constr_bc->SetConstrainedCoords(true, true, true, // x, y, z
117 true, true, true); // Rx, Ry, Rz
118
119 auto constr_d = chrono_types::make_shared<ChLinkMateGeneric>();
120 constr_d->Initialize(hnode1, mtruss, false, hnode1->Frame(), hnode1->Frame());
121 my_system.Add(constr_d);
122 constr_d->SetConstrainedCoords(false, true, true, // x, y, z
123 false, false, false); // Rx, Ry, Rz
124
125 //
126 // Add some EULER-BERNOULLI BEAMS (the fast way!)
127 //
128
129 // Shortcut!
130 // This ChBuilderBeamEuler helper object is very useful because it will
131 // subdivide 'beams' into sequences of finite elements of beam type, ex.
132 // one 'beam' could be made of 5 FEM elements of ChElementBeamEuler class.
133 // If new nodes are needed, it will create them for you.
134 ChBuilderBeamEuler builder;
135
136 // Now, simply use BuildBeam to create a beam from a point to another:
137 builder.BuildBeam(my_mesh, // the mesh where to put the created nodes and elements
138 msection, // the ChBeamSectionEulerAdvanced to use for the ChElementBeamEuler elements
139 5, // the number of ChElementBeamEuler to create
140 ChVector<>(0, 0, -0.1), // the 'A' point in space (beginning of beam)
141 ChVector<>(0.2, 0, -0.1), // the 'B' point in space (end of beam)
142 ChVector<>(0, 1, 0)); // the 'Y' up direction of the section for the beam
143
144 // After having used BuildBeam(), you can retrieve the nodes used for the beam,
145 // For example say you want to fix the A end and apply a force to the B end:
146 builder.GetLastBeamNodes().back()->SetFixed(true);
147 builder.GetLastBeamNodes().front()->SetForce(ChVector<>(0, -1, 0));
148
149 // Again, use BuildBeam for creating another beam, this time
150 // it uses one node (the last node created by the last beam) and one point:
151 builder.BuildBeam(my_mesh, msection, 5,
152 builder.GetLastBeamNodes().front(), // the 'A' node in space (beginning of beam)
153 ChVector<>(0.2, 0.1, -0.1), // the 'B' point in space (end of beam)
154 ChVector<>(0, 1, 0)); // the 'Y' up direction of the section for the beam
155
156 //
157 // Final touches..
158 //
159
160 // We do not want gravity effect on FEA elements in this demo
161 my_mesh->SetAutomaticGravity(false);
162
163 // Remember to add the mesh to the system!
164 my_system.Add(my_mesh);
165
166 // ==Asset== attach a visualization of the FEM mesh.
167 // This will automatically update a triangle mesh (a ChTriangleMeshShape
168 // asset that is internally managed) by setting proper
169 // coordinates and vertex colors as in the FEM elements.
170 // Such triangle mesh can be rendered by Irrlicht or POVray or whatever
171 // postprocessor that can handle a colored ChTriangleMeshShape).
172 // Do not forget AddAsset() at the end!
173
174 /*
175 auto mvisualizebeamA = chrono_types::make_shared<ChVisualizationFEAmesh>(*(my_mesh.get()));
176 mvisualizebeamA->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_SURFACE);
177 mvisualizebeamA->SetSmoothFaces(true);
178 my_mesh->AddAsset(mvisualizebeamA);
179 */
180
181 auto mvisualizebeamA = chrono_types::make_shared<ChVisualizationFEAmesh>(*(my_mesh.get()));
182 mvisualizebeamA->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_ELEM_BEAM_MZ);
183 mvisualizebeamA->SetColorscaleMinMax(-0.4, 0.4);
184 mvisualizebeamA->SetSmoothFaces(true);
185 mvisualizebeamA->SetWireframe(false);
186 my_mesh->AddAsset(mvisualizebeamA);
187
188 auto mvisualizebeamC = chrono_types::make_shared<ChVisualizationFEAmesh>(*(my_mesh.get()));
189 mvisualizebeamC->SetFEMglyphType(ChVisualizationFEAmesh::E_GLYPH_NODE_CSYS);
190 mvisualizebeamC->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_NONE);
191 mvisualizebeamC->SetSymbolsThickness(0.006);
192 mvisualizebeamC->SetSymbolsScale(0.01);
193 mvisualizebeamC->SetZbufferHide(false);
194 my_mesh->AddAsset(mvisualizebeamC);
195
196 // ==IMPORTANT!== Use this function for adding a ChIrrNodeAsset to all items
197 // in the system. These ChIrrNodeAsset assets are 'proxies' to the Irrlicht meshes.
198 // If you need a finer control on which item really needs a visualization proxy in
199 // Irrlicht, just use application.AssetBind(myitem); on a per-item basis.
200
201 application.AssetBindAll();
202
203 // ==IMPORTANT!== Use this function for 'converting' into Irrlicht meshes the assets
204 // that you added to the bodies into 3D shapes, they can be visualized by Irrlicht!
205
206 application.AssetUpdateAll();
207
208 // THE SIMULATION LOOP
209
210 auto solver = chrono_types::make_shared<ChSolverMINRES>();
211 my_system.SetSolver(solver);
212 solver->SetMaxIterations(500);
213 solver->SetTolerance(1e-14);
214 solver->EnableDiagonalPreconditioner(true);
215 solver->EnableWarmStart(true); // IMPORTANT for convergence when using EULER_IMPLICIT_LINEARIZED
216 solver->SetVerbose(false);
217
218 my_system.SetSolverForceTolerance(1e-13);
219
220 // Change type of integrator:
221 ////auto stepper = chrono_types::make_shared<ChTimestepperHHT>();
222 ////my_system.SetTimestepper(stepper);
223 ////stepper->SetAlpha(-0.2);
224 ////stepper->SetMaxiters(6);
225 ////stepper->SetAbsTolerances(1e-12);
226 ////stepper->SetVerbose(true);
227 ////stepper->SetStepControl(false);
228
229 my_system.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED);
230
231 application.SetTimestep(0.001);
232
233 GetLog() << "\n\n\n===========STATICS======== \n\n\n";
234
235 application.GetSystem()->DoStaticLinear();
236
237 GetLog() << "BEAM RESULTS (LINEAR STATIC ANALYSIS) \n\n";
238
239 ChVector<> F, M;
240 ChVectorDynamic<> displ;
241
242 belement1->GetStateBlock(displ);
243 GetLog() << displ;
244 for (double eta = -1; eta <= 1; eta += 0.4) {
245 belement1->EvaluateSectionForceTorque(eta, F, M);
246 GetLog() << " b1_at " << eta << " Mx=" << M.x() << " My=" << M.y() << " Mz=" << M.z() << " Tx=" << F.x()
247 << " Ty=" << F.y() << " Tz=" << F.z() << "\n";
248 }
249 GetLog() << "\n";
250 belement2->GetStateBlock(displ);
251 for (double eta = -1; eta <= 1; eta += 0.4) {
252 belement2->EvaluateSectionForceTorque(eta, F, M);
253 GetLog() << " b2_at " << eta << " Mx=" << M.x() << " My=" << M.y() << " Mz=" << M.z() << " Tx=" << F.x()
254 << " Ty=" << F.y() << " Tz=" << F.z() << "\n";
255 }
256
257 GetLog() << "Node 3 coordinate x= " << hnode3->Frame().GetPos().x() << " y=" << hnode3->Frame().GetPos().y()
258 << " z=" << hnode3->Frame().GetPos().z() << "\n\n";
259
260 GetLog() << "Press SPACE bar to start/stop dynamic simulation \n\n";
261 GetLog() << "Press F10 for nonlinear static solution \n\n";
262 GetLog() << "Press F11 for linear static solution \n\n";
263
264 while (application.GetDevice()->run()) {
265 application.BeginScene();
266
267 application.DrawAll();
268
269 application.DoStep();
270
271 application.EndScene();
272 }
273
274 return 0;
275 }
276