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 thin shells of Kirchhoff-Love type, with BST triangle finite elements
16 //
17 // =============================================================================
18
19 #include <vector>
20
21 #include "chrono/physics/ChBodyEasy.h"
22 #include "chrono/physics/ChLinkMate.h"
23 #include "chrono/physics/ChSystemSMC.h"
24 #include "chrono/solver/ChIterativeSolverLS.h"
25 #include "chrono/timestepper/ChTimestepper.h"
26
27 #include "chrono/fea/ChElementShellBST.h"
28 #include "chrono/fea/ChLinkPointFrame.h"
29 #include "chrono/fea/ChMesh.h"
30 #include "chrono/fea/ChVisualizationFEAmesh.h"
31 #include "chrono/fea/ChMeshFileLoader.h"
32
33 #include "chrono_irrlicht/ChIrrApp.h"
34
35 #include "chrono_pardisomkl/ChSolverPardisoMKL.h"
36
37 #include "chrono_postprocess/ChGnuPlot.h"
38
39 #include "chrono_thirdparty/filesystem/path.h"
40
41 // Remember to use the namespace 'chrono' because all classes
42 // of Chrono::Engine belong to this namespace and its children...
43
44 using namespace chrono;
45 using namespace chrono::fea;
46 using namespace chrono::irrlicht;
47 using namespace chrono::postprocess;
48 using namespace irr;
49
50 // Output directory
51 const std::string out_dir = GetChronoOutputPath() + "FEA_SHELLS";
52
main(int argc,char * argv[])53 int main(int argc, char* argv[]) {
54 GetLog() << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << "\n\n";
55
56 // Create (if needed) output directory
57 if (!filesystem::create_directory(filesystem::path(out_dir))) {
58 std::cout << "Error creating directory " << out_dir << std::endl;
59 return 1;
60 }
61
62 // Create a Chrono::Engine physical system
63 ChSystemSMC my_system;
64
65 // Create a mesh, that is a container for groups
66 // of elements and their referenced nodes.
67 auto my_mesh = chrono_types::make_shared<ChMesh>();
68
69 // Remember to add the mesh to the system!
70 my_system.Add(my_mesh);
71
72 // my_system.Set_G_acc(VNULL); // to remove gravity effect, or:
73 // my_mesh->SetAutomaticGravity(false);
74
75 std::shared_ptr<ChNodeFEAxyz> nodePlotA;
76 std::shared_ptr<ChNodeFEAxyz> nodePlotB;
77 std::vector<std::shared_ptr<ChNodeFEAxyz>> nodesLoad;
78
79 ChFunction_Recorder ref_X;
80 ChFunction_Recorder ref_Y;
81
82 ChVector<> load_force;
83
84 //
85 // BENCHMARK n.1
86 //
87 // Add a single BST element:
88 //
89
90 if (false) // set as 'true' to execute this
91 {
92 // Create a material
93
94 double density = 0.1;
95 double E = 1.2e3;
96 double nu = 0.3;
97 double thickness = 0.001;
98
99 auto melasticity = chrono_types::make_shared<ChElasticityKirchhoffIsothropic>(E, nu);
100 auto material = chrono_types::make_shared<ChMaterialShellKirchhoff>(melasticity);
101 material->SetDensity(density);
102
103 // Create nodes
104 double L = 1.0;
105
106 ChVector<> p0(0, 0, 0);
107 ChVector<> p1(L, 0, 0);
108 ChVector<> p2(0, L, 0);
109 ChVector<> p3(L, L, 0);
110 ChVector<> p4(-L, L, 0);
111 ChVector<> p5(L, -L, 0);
112
113 auto mnode0 = chrono_types::make_shared<ChNodeFEAxyz>(p0);
114 auto mnode1 = chrono_types::make_shared<ChNodeFEAxyz>(p1);
115 auto mnode2 = chrono_types::make_shared<ChNodeFEAxyz>(p2);
116 auto mnode3 = chrono_types::make_shared<ChNodeFEAxyz>(p3);
117 auto mnode4 = chrono_types::make_shared<ChNodeFEAxyz>(p4);
118 auto mnode5 = chrono_types::make_shared<ChNodeFEAxyz>(p5);
119 my_mesh->AddNode(mnode0);
120 my_mesh->AddNode(mnode1);
121 my_mesh->AddNode(mnode2);
122 my_mesh->AddNode(mnode3);
123 my_mesh->AddNode(mnode4);
124 my_mesh->AddNode(mnode5);
125
126 // Create element
127
128 auto melement = chrono_types::make_shared<ChElementShellBST>();
129 my_mesh->AddElement(melement);
130
131 melement->SetNodes(mnode0, mnode1, mnode2, nullptr, nullptr, nullptr); // mnode3, mnode4, mnode5);
132
133 melement->AddLayer(thickness, 0 * CH_C_DEG_TO_RAD, material);
134
135 // TEST
136 my_system.Setup();
137 my_system.Update();
138 GetLog() << "BST initial: \n"
139 << "Area: " << melement->area << "\n"
140 << "l0: " << melement->l0 << "\n"
141 << "phi0: " << melement->phi0 << "\n"
142 << "k0: " << melement->k0 << "\n"
143 << "e0: " << melement->e0 << "\n";
144
145 mnode1->SetPos(mnode1->GetPos() + ChVector<>(0.1, 0, 0));
146
147 my_system.Update();
148 ChVectorDynamic<double> Fi(melement->GetNdofs());
149 melement->ComputeInternalForces(Fi);
150 GetLog() << "BST updated: \n"
151 << "phi: " << melement->phi << "\n"
152 << "k: " << melement->k << "\n"
153 << "e: " << melement->e << "\n"
154 << "m: " << melement->m << "\n"
155 << "n: " << melement->n << "\n";
156 ChVector<> resultant = VNULL;
157 GetLog() << "Fi= \n";
158 for (int i = 0; i < Fi.size(); ++i) {
159 if (i % 3 == 0) {
160 GetLog() << "-------" << i / 3 << "\n";
161 ChVector<> fi = Fi.segment(i, 3);
162 resultant += fi;
163 }
164 GetLog() << Fi(i) << "\n";
165 }
166 GetLog() << "resultant: " << resultant << "\n";
167 // system("pause");
168 }
169
170 //
171 // BENCHMARK n.2
172 //
173 // Add a rectangular mesh of BST elements:
174 //
175
176 std::shared_ptr<ChNodeFEAxyz> mnodemonitor;
177 std::shared_ptr<ChElementShellBST> melementmonitor;
178
179 if (true) // set as 'true' to execute this
180 {
181 // Create a material
182
183 double density = 100;
184 double E = 6e4;
185 double nu = 0.0;
186 double thickness = 0.01;
187
188 auto melasticity = chrono_types::make_shared<ChElasticityKirchhoffIsothropic>(E, nu);
189 auto material = chrono_types::make_shared<ChMaterialShellKirchhoff>(melasticity);
190 material->SetDensity(density);
191
192 double L_x = 1;
193 size_t nsections_x = 40;
194 double L_z = 1;
195 size_t nsections_z = 40;
196
197 // Create nodes
198 std::vector<std::shared_ptr<ChNodeFEAxyz>> mynodes; // for future loop when adding elements
199 for (size_t iz = 0; iz <= nsections_z; ++iz) {
200 for (size_t ix = 0; ix <= nsections_x; ++ix) {
201 ChVector<> p(ix * (L_x / nsections_x), 0, iz * (L_z / nsections_z));
202
203 auto mnode = chrono_types::make_shared<ChNodeFEAxyz>(p);
204
205 my_mesh->AddNode(mnode);
206
207 mynodes.push_back(mnode);
208 }
209 }
210 // Create elements
211 for (size_t iz = 0; iz < nsections_z; ++iz) {
212 for (size_t ix = 0; ix < nsections_x; ++ix) {
213 auto melementA = chrono_types::make_shared<ChElementShellBST>();
214 my_mesh->AddElement(melementA);
215
216 if (iz == 0 && ix == 1)
217 melementmonitor = melementA;
218
219 std::shared_ptr<ChNodeFEAxyz> boundary_1;
220 std::shared_ptr<ChNodeFEAxyz> boundary_2;
221 std::shared_ptr<ChNodeFEAxyz> boundary_3;
222
223 boundary_1 = mynodes[(iz + 1) * (nsections_x + 1) + ix + 1];
224 if (ix > 0)
225 boundary_2 = mynodes[(iz + 1) * (nsections_x + 1) + ix - 1];
226 else
227 boundary_2 = nullptr;
228 if (iz > 0)
229 boundary_3 = mynodes[(iz - 1) * (nsections_x + 1) + ix + 1];
230 else
231 boundary_3 = nullptr;
232
233 melementA->SetNodes(mynodes[(iz) * (nsections_x + 1) + ix], mynodes[(iz) * (nsections_x + 1) + ix + 1],
234 mynodes[(iz + 1) * (nsections_x + 1) + ix], boundary_1, boundary_2, boundary_3);
235
236 melementA->AddLayer(thickness, 0 * CH_C_DEG_TO_RAD, material);
237
238 auto melementB = chrono_types::make_shared<ChElementShellBST>();
239 my_mesh->AddElement(melementB);
240
241 boundary_1 = mynodes[(iz) * (nsections_x + 1) + ix];
242 if (ix < nsections_x - 1)
243 boundary_2 = mynodes[(iz) * (nsections_x + 1) + ix + 2];
244 else
245 boundary_2 = nullptr;
246 if (iz < nsections_z - 1)
247 boundary_3 = mynodes[(iz + 2) * (nsections_x + 1) + ix];
248 else
249 boundary_3 = nullptr;
250
251 melementB->SetNodes(mynodes[(iz + 1) * (nsections_x + 1) + ix + 1],
252 mynodes[(iz + 1) * (nsections_x + 1) + ix],
253 mynodes[(iz) * (nsections_x + 1) + ix + 1], boundary_1, boundary_2, boundary_3);
254
255 melementB->AddLayer(thickness, 0 * CH_C_DEG_TO_RAD, material);
256 }
257 }
258
259 // fix some nodes
260 for (int j = 0; j < 30; ++j) {
261 for (int k = 0; k < 30; ++k) {
262 mynodes[j * (nsections_x + 1) + k]->SetFixed(true);
263 }
264 }
265 /*
266 mynodes[0*(nsections_x+1)+1]->SetFixed(true);
267 mynodes[1*(nsections_x+1)+1]->SetFixed(true);
268
269 mynodes[0*(nsections_x+1)+nsections_x-1]->SetFixed(true);
270 mynodes[1*(nsections_x+1)+nsections_x-1]->SetFixed(true);
271
272 mynodes[0*(nsections_x+1)+nsections_x]->SetFixed(true);
273 mynodes[1*(nsections_x+1)+nsections_x]->SetFixed(true);
274 */
275 }
276
277 //
278 // BENCHMARK n.3
279 //
280 // Load and create shells from a .obj file containing a triangle mesh surface
281 //
282
283 if (false) {
284 double density = 100;
285 double E = 6e5;
286 double nu = 0.0;
287 double thickness = 0.01;
288
289 auto melasticity = chrono_types::make_shared<ChElasticityKirchhoffIsothropic>(E, nu);
290 auto material = chrono_types::make_shared<ChMaterialShellKirchhoff>(melasticity);
291 material->SetDensity(density);
292
293 ChMeshFileLoader::BSTShellFromObjFile(my_mesh, GetChronoDataFile("models/cube.obj").c_str(), material, thickness);
294 }
295
296 // ==Asset== attach a visualization of the FEM mesh.
297 // This will automatically update a triangle mesh (a ChTriangleMeshShape
298 // asset that is internally managed) by setting proper
299 // coordinates and vertex colors as in the FEM elements.
300 // Such triangle mesh can be rendered by Irrlicht or POVray or whatever
301 // postprocessor that can handle a colored ChTriangleMeshShape).
302 // Do not forget AddAsset() at the end!
303
304 auto mvisualizeshellA = chrono_types::make_shared<ChVisualizationFEAmesh>(*(my_mesh.get()));
305 // mvisualizeshellA->SetSmoothFaces(true);
306 // mvisualizeshellA->SetWireframe(true);
307 mvisualizeshellA->SetShellResolution(2);
308 // mvisualizeshellA->SetBackfaceCull(true);
309 my_mesh->AddAsset(mvisualizeshellA);
310
311 auto mvisualizeshellB = chrono_types::make_shared<ChVisualizationFEAmesh>(*(my_mesh.get()));
312 mvisualizeshellB->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_NONE);
313 mvisualizeshellB->SetFEMglyphType(ChVisualizationFEAmesh::E_GLYPH_NODE_DOT_POS);
314 mvisualizeshellB->SetSymbolsThickness(0.006);
315 my_mesh->AddAsset(mvisualizeshellB);
316
317 //
318 // VISUALIZATION
319 //
320
321 // Create the Irrlicht visualization (open the Irrlicht device,
322 // bind a simple user interface, etc. etc.)
323 ChIrrApp application(&my_system, L"Shells FEA test: triangle BST elements", core::dimension2d<u32>(1024, 768));
324
325 // Easy shortcuts to add camera, lights, logo and sky in Irrlicht scene:
326 application.AddTypicalLogo();
327 application.AddTypicalSky();
328 application.AddLightWithShadow(irr::core::vector3df(2, 2, 2), irr::core::vector3df(0, 0, 0), 6, 0.2, 6, 50);
329 application.AddLight(irr::core::vector3df(-2, -2, 0), 6, irr::video::SColorf(0.6f, 1.0f, 1.0f, 1.0f));
330 application.AddLight(irr::core::vector3df(0, -2, -2), 6, irr::video::SColorf(0.6f, 1.0f, 1.0f, 1.0f));
331 // application.AddTypicalLights();
332 application.AddTypicalCamera(core::vector3df(1.f, 0.3f, 1.3f), core::vector3df(0.5f, -0.3f, 0.5f));
333
334 // ==IMPORTANT!== Use this function for adding a ChIrrNodeAsset to all items
335 // in the system. These ChIrrNodeAsset assets are 'proxies' to the Irrlicht meshes.
336 // If you need a finer control on which item really needs a visualization proxy in
337 // Irrlicht, just use application.AssetBind(myitem); on a per-item basis.
338
339 application.AssetBindAll();
340
341 // ==IMPORTANT!== Use this function for 'converting' into Irrlicht meshes the assets
342 // that you added to the bodies into 3D shapes, they can be visualized by Irrlicht!
343
344 application.AssetUpdateAll();
345
346 application.AddShadowAll();
347 //
348 // THE SOFT-REAL-TIME CYCLE
349 //
350
351 // Change solver to PardisoMKL
352 auto mkl_solver = chrono_types::make_shared<ChSolverPardisoMKL>();
353 mkl_solver->LockSparsityPattern(true);
354 my_system.SetSolver(mkl_solver);
355
356 // auto gmres_solver = chrono_types::make_shared<ChSolverGMRES>();
357 // gmres_solver->SetMaxIterations(50);
358 // my_system.SetSolver(gmres_solver);
359
360 double timestep = 0.005;
361 application.SetTimestep(timestep);
362 my_system.Setup();
363 my_system.Update();
364
365 ChFunction_Recorder rec_X;
366 ChFunction_Recorder rec_Y;
367
368 while (application.GetDevice()->run()) {
369 application.BeginScene();
370
371 application.DrawAll();
372
373 application.DoStep();
374
375 application.EndScene();
376 }
377
378 return 0;
379 }
380