Hi dev team/community, Recently I've been trying to do some benchmarking of Chrono and other libraries as FEA frameworks for use at a marine engineering company. In particular, we are interested in modeling algae as long cable-like elements.
I've created a test-model where numerous multi-element lines hang from a single "backbone" line with one free end (see model6 in the attached .h file). Briefly, I create a backbone line with fixed ends using 1m long ChSpringElements to link the nodes; I then compute nodal positions for each of the hanging lines, each with 10 1m ShSpringElement segments per-line. I've also set up this simulation in MoorDyn (an open source mooring model), MSC Marc, and Orcaflex. Unfortunately, Chrono appears to be much slower than the other options: running a 16s simulation of a system with 100 hanging lines, for instance takes 3058 seconds in Chrono using the HHT timestepper and MINRES solver, 533 seconds in the alpha version of MoorDyn (using an RK2 stepper) and just 288 seconds in Orcaflex. With pending improvements to MoorDyn (including RK4/RKF45 integrators) I've been able to knock that simulation time down to just 31 seconds for MoorDyn. Given this large discrepancy in runtimes, I was curious if: 1. there were obvious bugs with my model formulation, 2. if there are any standard steps that can be taken to speed up the FEA integration in chrono, and 3. Whether anyone has ideas as to why these massive discrepancies might exist. Some basic profiling suggests that the solver is the computational bottleneck - I've tried swapping out all of the default solvers, but haven't seen huge increases in speed and am not familiar enough with the math of the solvers to delve deeply into the implementations. I've attached my .h file defining my models (model6 is the model in question) and the wrapper program I am using to run my simulation (my_example.cpp). Any feedback/ideas are greatly appreciated. Thanks, David -- You received this message because you are subscribed to the Google Groups "ProjectChrono" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/projectchrono/2920668d-0033-4dd3-aea1-9c65a5e3191bn%40googlegroups.com.
// ============================================================================= // PROJECT CHRONO - http://projectchrono.org // // Copyright (c) 2014 projectchrono.org // All rights reserved. // // Use of this source code is governed by a BSD-style license that can be found // in the LICENSE file at the top level of the distribution and at // http://projectchrono.org/license-chrono.txt. // // ============================================================================= // Authors: Alessandro Tasora, Radu Serban // ============================================================================= // // FEA for 3D beams of 'cable' type (ANCF gradient-deficient beams) // // ============================================================================= #define WRITE_OUTFILE 1 #include <ctime> #include "chrono/physics/ChSystemSMC.h" #include "chrono/solver/ChDirectSolverLS.h" #include "chrono/solver/ChIterativeSolverLS.h" #include "chrono/timestepper/ChTimestepper.h" #include "chrono_irrlicht/ChVisualSystemIrrlicht.h" #include "myFEAcables.h" using namespace chrono; using namespace chrono::fea; using namespace chrono::irrlicht; using namespace irr; // Select solver type (SPARSE_QR, SPARSE_LU, or MINRES). ChSolver::Type solver_type = ChSolver::Type::MINRES; int main(int argc, char* argv[]) { GetLog() << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << "\n\n"; // Settings bool do_animate = true; double current_time = 0; double end_time = 16; double time_step = 0.0002; int number_of_chains = 10; int segments_per_chain = 10; // Create a Chrono::Engine physical system ChSystemSMC sys; // Create a mesh, that is a container for groups of elements and // their referenced nodes. auto mesh = chrono_types::make_shared<ChMesh>(); // Create one of the available models (defined in FEAcables.h) //auto model = Model1(sys, mesh); auto model = Model6(sys, mesh, number_of_chains, segments_per_chain); //auto model = Model5(sys, mesh, number_of_chains, segments_per_chain); //auto model = Model3(sys, mesh, number_of_chains); // Remember to add the mesh to the system! sys.Add(mesh); // Visualization of the FEM mesh. // This will automatically update a triangle mesh (a ChTriangleMeshShape asset that is internally managed) by // setting proper coordinates and vertex colors as in the FEM elements. Such triangle mesh can be rendered by // Irrlicht or POVray or whatever postprocessor that can handle a colored ChTriangleMeshShape). //auto vis_beam_A = chrono_types::make_shared<ChVisualShapeFEA>(mesh); //vis_beam_A->SetFEMdataType(ChVisualShapeFEA::DataType::ELEM_BEAM_MZ); //vis_beam_A->SetColorscaleMinMax(-0.4, 0.4); //vis_beam_A->SetSmoothFaces(true); //vis_beam_A->SetWireframe(true); //mesh->AddVisualShapeFEA(vis_beam_A); //auto vis_beam_B = chrono_types::make_shared<ChVisualShapeFEA>(mesh); //vis_beam_B->SetFEMglyphType(ChVisualShapeFEA::GlyphType::NODE_DOT_POS); //vis_beam_B->SetFEMdataType(ChVisualShapeFEA::DataType::NONE); //vis_beam_B->SetSymbolsThickness(0.06); //vis_beam_B->SetSymbolsScale(0.01); //vis_beam_B->SetZbufferHide(false); //mesh->AddVisualShapeFEA(vis_beam_B); // Set solver and solver settings switch (solver_type) { case ChSolver::Type::SPARSE_QR: { std::cout << "Using SparseQR solver" << std::endl; auto solver = chrono_types::make_shared<ChSolverSparseQR>(); sys.SetSolver(solver); solver->UseSparsityPatternLearner(true); solver->LockSparsityPattern(true); solver->SetVerbose(false); break; } case ChSolver::Type::SPARSE_LU: { std::cout << "Using SparseLU solver" << std::endl; auto solver = chrono_types::make_shared<ChSolverSparseLU>(); sys.SetSolver(solver); solver->UseSparsityPatternLearner(true); solver->LockSparsityPattern(true); solver->SetVerbose(false); break; } case ChSolver::Type::MINRES: { std::cout << "Using MINRES solver" << std::endl; auto solver = chrono_types::make_shared<ChSolverMINRES>(); sys.SetSolver(solver); solver->SetMaxIterations(200); solver->SetTolerance(1e-10); solver->EnableDiagonalPreconditioner(true); solver->EnableWarmStart(true); // IMPORTANT for convergence when using EULER_IMPLICIT_LINEARIZED solver->SetVerbose(false); break; } default: { std::cout << "Solver type not supported." << std::endl; break; } } // Create the Irrlicht visualization system //auto vis = chrono_types::make_shared<ChVisualSystemIrrlicht>(); //vis->SetWindowSize(800, 600); //vis->SetWindowTitle("Cables FEM"); //vis->Initialize(); //vis->AddLogo(); //vis->AddSkyBox(); //vis->AddTypicalLights(); ////vis->AddCamera(ChVector<>(0, 0.6, -1.0)); //vis->AddCamera(ChVector<>(0, -1.0, 0.0), ChVector<>(0,0,0)); //sys.SetVisualSystem(vis); // Set integrator //sys.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED); sys.SetTimestepperType(ChTimestepper::Type::RUNGEKUTTA45); //sys.SetTimestepperType(ChTimestepper::Type::NEWMARK); //if (auto mystepper = std::dynamic_pointer_cast<ChTimestepperNewmark>(sys.GetTimestepper())) { // mystepper->SetGammaBeta(0.5, 0.5); //} //sys.SetTimestepperType(ChTimestepper::Type::HHT); //if (auto mystepper = std::dynamic_pointer_cast<ChTimestepperHHT>(sys.GetTimestepper())) { // mystepper->SetAlpha(-0.2); // mystepper->SetMinStepSize(0.00000001); // mystepper->SetStepIncreaseFactor(1.2); // mystepper->SetStepDecreaseFactor(0.8); //} sys.SetSolverForceTolerance(1e-13); std::cout << "Num threads: " << sys.GetNumThreadsChrono() << std::endl; sys.SetNumThreads(36); std::cout << "Num threads: " << sys.GetNumThreadsChrono() << std::endl; //std::system("pause"); time_t start = time(nullptr); ofstream outfile; int step = 0; int frames_per_sec = 30; int output_step = (end_time / time_step) / (end_time * frames_per_sec); // Output so as to have frames_per_sec outputs while (current_time < end_time) { //vis->BeginScene(); //vis->DrawAll(); //vis->EndScene(); sys.DoStepDynamics(time_step); current_time += time_step; std::cout << "Simulation Time " << current_time << '\n'; #ifdef WRITE_OUTFILE // Print node positions: if (step % output_step == 0) { outfile.open("longline.out", std::ios_base::app); outfile << current_time << "\t"; for (unsigned int inode = 0; inode < mesh->GetNnodes(); ++inode) { if (auto mnode = std::dynamic_pointer_cast<ChNodeFEAxyz>(mesh->GetNode(inode))) { /*if (mnode->GetPos().x() < 0.01) { GetLog() << "Node at y=" << mnode->GetPos().y() << " has T=" << mnode->GetP() << "\n"; }*/ outfile << "(" << mnode->GetPos().x() << ", " << mnode->GetPos().y() << ", " << mnode->GetPos().z() << ")\t"; } } outfile << "\n"; outfile.close(); } #endif ++step; //std::cout << "Wall Clock time: " << clock() / CLOCKS_PER_SEC << " seconds" << std::endl; } time_t end = time(nullptr); std::cout << "Finished time: " << end - start << " seconds" << std::endl; std::system("pause"); return 0; }
// ============================================================================= // PROJECT CHRONO - http://projectchrono.org // // Copyright (c) 2014 projectchrono.org // All rights reserved. // // Use of this source code is governed by a BSD-style license that can be found // in the LICENSE file at the top level of the distribution and at // http://projectchrono.org/license-chrono.txt. // // ============================================================================= // // Models using the ANCF gradient-deficient cable element // // ============================================================================= #include "chrono/physics/ChSystem.h" #include "chrono/physics/ChBodyEasy.h" #include "chrono/fea/ChElementCableANCF.h" #include "chrono/fea/ChBuilderBeam.h" #include "chrono/fea/ChMesh.h" #include "chrono/assets/ChVisualShapeFEA.h" #include "chrono/fea/ChLinkPointFrame.h" #include "chrono/fea/ChLinkDirFrame.h" #include "chrono/fea/ChElementSpring.h" #include "chrono/fea/ChElementBar.h" #include <iostream> #include <fstream> #include <ctime> #include <cmath> using namespace chrono; using namespace chrono::fea; using namespace chrono::irrlicht; using namespace irr; using namespace std; const int N_SEGMENTS = 10; // ---------------------------------------------------------------------------- // Model1: A beam composed of a single ANCF cable element, with one end fixed // and a constant force applied at the other node. A rigid body is // connected to node2. // ---------------------------------------------------------------------------- class Model1 { public: Model1(ChSystem& system, std::shared_ptr<ChMesh> mesh) { double beam_L = 0.1; double beam_diameter = 0.015; // Create a section, i.e. thickness and material properties // for beams. This will be shared among some beams. auto msection_cable = chrono_types::make_shared<ChBeamSectionCable>(); msection_cable->SetDiameter(beam_diameter); msection_cable->SetYoungModulus(0.01e9); msection_cable->SetBeamRaleyghDamping(0.000); // Create the nodes auto hnodeancf1 = chrono_types::make_shared<ChNodeFEAxyzD>(ChVector<>(0, 0, -0.2), ChVector<>(1, 0, 0)); auto hnodeancf2 = chrono_types::make_shared<ChNodeFEAxyzD>(ChVector<>(beam_L, 0, -0.2), ChVector<>(1, 0, 0)); mesh->AddNode(hnodeancf1); mesh->AddNode(hnodeancf2); // Create the element auto belementancf1 = chrono_types::make_shared<ChElementCableANCF>(); belementancf1->SetNodes(hnodeancf1, hnodeancf2); belementancf1->SetSection(msection_cable); mesh->AddElement(belementancf1); // Apply a force or a torque to a node: hnodeancf2->SetForce(ChVector<>(0, 3, 0)); hnodeancf1->SetFixed(true); // Add a rigid body connected to the end of the beam: body = chrono_types::make_shared<ChBodyEasyBox>(0.1, 0.02, 0.02, 1000); body->SetPos(hnodeancf2->GetPos() + ChVector<>(0.05, 0, 0)); system.Add(body); auto constraint_pos = chrono_types::make_shared<ChLinkPointFrame>(); constraint_pos->Initialize(hnodeancf2, body); system.Add(constraint_pos); auto constraint_dir = chrono_types::make_shared<ChLinkDirFrame>(); constraint_dir->Initialize(hnodeancf2, body); constraint_dir->SetDirectionInAbsoluteCoords(ChVector<>(1, 0, 0)); system.Add(constraint_dir); } // Print position of end body void PrintBodyPosition() { std::cout << "Time: " << body->GetChTime() << std::endl; std::cout << " " << body->GetPos() << std::endl; } private: std::shared_ptr<ChBodyEasyBox> body; }; // ---------------------------------------------------------------------------- // Model2: A beam composed of 10 ANCF beam element, with one end hinged to // ground, moving under gravity alone. // This model demonstrates the use of the utility class ChBuilderCableANCF. // ---------------------------------------------------------------------------- class Model2 { public: Model2(ChSystem& system, std::shared_ptr<ChMesh> mesh) { // Create a section, i.e. thickness and material properties // for beams. This will be shared among some beams. auto msection_cable2 = chrono_types::make_shared<ChBeamSectionCable>(); msection_cable2->SetDiameter(0.015); msection_cable2->SetYoungModulus(0.01e9); msection_cable2->SetBeamRaleyghDamping(0.000); // This ChBuilderCableANCF helper object is very useful because it will // subdivide 'beams' into sequences of finite elements of beam type, ex. // one 'beam' could be made of 5 FEM elements of ChElementBeamANCF_3333 class. // If new nodes are needed, it will create them for you. ChBuilderCableANCF builder; // Now, simply use BuildBeam to create a beam from a point to another: builder.BuildBeam(mesh, // the mesh where to put the created nodes and elements msection_cable2, // the ChBeamSectionCable to use for the ChElementBeamANCF_3333 elements 10, // the number of ChElementBeamANCF_3333 to create ChVector<>(0, 0, -0.1), // the 'A' point in space (beginning of beam) ChVector<>(0.5, 0, -0.1)); // the 'B' point in space (end of beam) // After having used BuildBeam(), you can retrieve the nodes used for the beam, // For example say you want to fix both pos and dir of A end and apply a force to the B end: // builder.GetLastBeamNodes().back()->SetFixed(true); builder.GetLastBeamNodes().front()->SetForce(ChVector<>(0, -0.2, 0)); // For instance, now retrieve the A end and add a constraint to // block the position only of that node: auto mtruss = chrono_types::make_shared<ChBody>(); mtruss->SetBodyFixed(true); auto constraint_hinge = chrono_types::make_shared<ChLinkPointFrame>(); constraint_hinge->Initialize(builder.GetLastBeamNodes().back(), mtruss); system.Add(constraint_hinge); } }; // ---------------------------------------------------------------------------- // Model3: A set of beam elements with connected bodies, each with different // number of ANCF cable elements. // ---------------------------------------------------------------------------- class Model3 { public: Model3(ChSystem& system, std::shared_ptr<ChMesh> mesh, int n_chains = 10) : bodies(n_chains) { auto msection_cable2 = chrono_types::make_shared<ChBeamSectionCable>(); msection_cable2->SetArea(1.40e-4); msection_cable2->SetI(2e-9); msection_cable2->SetDiameter(13.332e-3); msection_cable2->SetYoungModulus(53.7e9); msection_cable2->SetDensity(7880.0); msection_cable2->SetBeamRaleyghDamping(0.000); msection_cable2->SetDrawCircularRadius(0.08); auto mtruss = chrono_types::make_shared<ChBody>(); mtruss->SetBodyFixed(true); for (int j = 0; j < n_chains; ++j) { ChBuilderCableANCF builder; // Now, simply use BuildBeam to create a beam from a point to another: builder.BuildBeam(mesh, // the mesh where to put the created nodes and elements msection_cable2, // ChBeamSectionCable to use for the ChElementBeamANCF_3333 elements N_SEGMENTS, // number of ChElementBeamANCF_3333 to create ChVector<>(5.0 * j, -16.6, 10.0), // point A (beginning of beam) ChVector<>(5.0 * j, -25.26025404, 15.0) // point B (end of beam) ); builder.GetLastBeamNodes().back()->SetForce(ChVector<>(0, -0.2, 0)); auto constraint_hinge = chrono_types::make_shared<ChLinkPointFrame>(); constraint_hinge->Initialize(builder.GetLastBeamNodes().front(), mtruss); system.Add(constraint_hinge); auto msphere = chrono_types::make_shared<ChSphereShape>(); msphere->GetSphereGeometry().rad = 0.02; constraint_hinge->AddVisualShape(msphere); } } // Print positions of end bodies in each chain void PrintBodyPositions() { std::cout << "Time: " << bodies[0]->GetChTime() << std::endl; for (auto body : bodies) { std::cout << " " << body->GetPos() << std::endl; } } private: std::vector<std::shared_ptr<ChBodyEasyBox>> bodies; // end bodies }; // ---------------------------------------------------------------------------- // Model4: A set of beam elements with connected bodies, each with different // number of ANCF cable elements. // // A swinging chain for Kelson Marine AM_OpenTools speed test // Chain composed of ANCF beam elements // [X] Mass // [X] Diameter // [X] Young's modulus // [ ] Fluid drag // [ ] Fluid inertia/added mass // [X] Length // [X] Gravity // This model demonstrates the use of the utility class ChBuilderCableANCF. // ---------------------------------------------------------------------------- class Model4 { public: Model4(ChSystem& system, std::shared_ptr<ChMesh> mesh, int n_chains = 6, int elements = 10) : bodies(n_chains) { auto msection_cable2 = chrono_types::make_shared<ChBeamSectionCable>(); msection_cable2->SetArea(1.40e-4); msection_cable2->SetI(2e-9); msection_cable2->SetDiameter(13.332e-3); msection_cable2->SetDensity(7880); msection_cable2->SetYoungModulus(53.7e9); msection_cable2->SetBeamRaleyghDamping(0.000); auto mtruss = chrono_types::make_shared<ChBody>(); mtruss->SetBodyFixed(true); for (int j = 0; j < n_chains; ++j) { ChBuilderCableANCF builder; // Now, simply use BuildBeam to create a beam from a point to another: builder.BuildBeam(mesh, // the mesh where to put the created nodes and elements msection_cable2, // ChBeamSectionCable to use for the ChElementBeamANCF_3333 elements elements, // number of ChElementBeamANCF_3333 to create ChVector<>(5.0 * j, -16.6, 10.0), // point A (beginning of beam) ChVector<>(5.0 * j, -25.26025404, 15.0) // point B (end of beam) ); //// Now, simply use BuildBeam to create a beam from a point to another: //builder.BuildBeam(mesh, // the mesh where to put the created nodes and elements // msection_cable2, // the ChBeamSectionCable to use for the ChElementBeamANCF_3333 elements // 10, // the number of ChElementBeamANCF_3333 to create // ChVector<>(0, 0, -0.1), // the 'A' point in space (beginning of beam) // ChVector<>(0.5, 0, -0.1)); // the 'B' point in space (end of beam) builder.GetLastBeamNodes().back()->SetForce(ChVector<>(0, -107.91, 0)); // Add a hinged constraint to the top of each cable - the cable will swing from this point auto constraint_hinge = chrono_types::make_shared<ChLinkPointFrame>(); constraint_hinge->Initialize(builder.GetLastBeamNodes().front(), mtruss); system.Add(constraint_hinge); // Visualize the hinge with a sphere auto msphere = chrono_types::make_shared<ChSphereShape>(); msphere->GetSphereGeometry().rad = 0.02; constraint_hinge->AddVisualShape(msphere); //// make a box and connect it //bodies[j] = chrono_types::make_shared<ChBodyEasyBox>(0.04, 0.04, 0.04, 1); //bodies[j]->SetPos(builder.GetLastBeamNodes().back()->GetPos() + ChVector<>(0.1, 0, 0)); //system.Add(bodies[j]); //auto constraint_pos3 = chrono_types::make_shared<ChLinkPointFrame>(); //constraint_pos3->Initialize(builder.GetLastBeamNodes().back(), bodies[j]); //system.Add(constraint_pos3); //auto constraint_dir3 = chrono_types::make_shared<ChLinkDirFrame>(); //constraint_dir3->Initialize(builder.GetLastBeamNodes().back(), bodies[j]); //constraint_dir3->SetDirectionInAbsoluteCoords(ChVector<>(1, 0, 0)); //system.Add(constraint_dir3); } } // Print positions of end bodies in each chain double PrintBodyPositions() { std::cout << "Time: " << bodies[0]->GetChTime() << std::endl; ofstream myfile("out.csv", std::ios_base::app); if (myfile.is_open()) { myfile << bodies[0]->GetChTime(); myfile << clock() / CLOCKS_PER_SEC; for (auto body : bodies) { myfile << " " << body->GetPos(); } myfile << std::endl; myfile.close(); } else cout << "Unable to open file"; return bodies[0]->GetChTime(); } private: std::vector<std::shared_ptr<ChBodyEasyBox>> bodies; // end bodies }; // ---------------------------------------------------------------------------- // ---------------------------------------------------------------------------- // Model5: A set of spring elements with connected bodies // // A swinging chain for Kelson Marine AM_OpenTools speed test // Chain composed of spring elements // [X] Mass // [X] Diameter // [X] Young's modulus // [ ] Fluid drag // [ ] Fluid inertia/added mass // [X] Length // [X] Gravity // // Reference for spring-node implementation: // https://api.projectchrono.org/tutorial_demo__f_e_a_basic.html // ---------------------------------------------------------------------------- class Model5 { public: Model5(ChSystem& system, std::shared_ptr<ChMesh> mesh, int n_chains = 6, int elements = 10) : nodes(elements + 1) { // Global constants double PI = 3.14159; // Material inputs // Note for this simple test I am ignoring bending stiffness double chain_length = 10; // m double chain_mass_length = 1.1; // kg/m double chain_density = 7880; // kg/m^3 double chain_diameter = 13.332e-3; // m double chain_youngs_modulus = 53.7e9; // Pa // Material calculations double chain_cross_section_area = pow(chain_diameter, 2) * PI / 4; double element_length = chain_length / elements; double node_mass = chain_length / (elements + 1) * chain_mass_length; double element_stiffness = chain_youngs_modulus * chain_cross_section_area / element_length; auto mtruss = chrono_types::make_shared<ChBody>(); mtruss->SetBodyFixed(true); // Create the nodes for (int j = 0; j < n_chains; ++j) { for (int i = 0; i <= elements; ++i) { // We need to specify each nodes xyz in this formulation (inconvenient) // To make this easy, we say that for the purposes of the test lines // hang from the x-axis and nodes are rotated about the x-axis using the R_X // rotation matrix (see: https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations) double x{ 5.0 * j }; double y{ -1.0 * i * cos(0.523599) }; // 0.523599 == the angle of rotation from normal (radians) double z{ -1.0 * i * sin(0.523599) }; // This is really just SOHCAHTOA std::cout << "(" << x << ", " << y << ", " << z << ")\n"; nodes[i] = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(x, y, z)); nodes[i]->SetMass(node_mass); mesh->AddNode(nodes[i]); if (i == 0) { //// Add a hinged constraint to the top of each cable - the cable will swing from this point auto constraint_hinge = chrono_types::make_shared<ChLinkPointFrame>(); constraint_hinge->Initialize(nodes[i], mtruss); system.Add(constraint_hinge); // Visualize the hinge with a sphere auto msphere = chrono_types::make_shared<ChSphereShape>(); msphere->GetSphereGeometry().rad = 0.02; constraint_hinge->AddVisualShape(msphere); } else { //// Add a bar element to connect current node to previous one & add linear stiffness //auto melement = chrono_types::make_shared<ChElementBar>(); //melement->SetNodes(nodes[i], nodes[i-1]); //melement->SetBarArea(chain_cross_section_area); //melement->SetBarDensity(chain_density); //melement->SetBarYoungModulus(chain_youngs_modulus); //melement->SetBarRaleyghDamping(0); //mesh->AddElement(melement); // Add a spring connecting the current node to the previous one auto melement = chrono_types::make_shared<ChElementSpring>(); melement->SetNodes(nodes[i], nodes[i - 1]); melement->SetSpringK(element_stiffness); mesh->AddElement(melement); } }; } } // Print positions of end bodies in each chain double PrintBodyPositions() { std::cout << "test" << std::endl; nodes[0]->GetMass(); //std::cout << " " << position[0] << std::endl; //std::cout << " " << nodes[0]->GetX0()[0]; //std::cout << "Time: " << nodes[0]->GetChTime() << std::endl; //ofstream myfile("out.csv", std::ios_base::app); //if (myfile.is_open()) //{ // myfile << bodies[0]->GetChTime(); // myfile << clock() / CLOCKS_PER_SEC; // for (auto body : bodies) { // myfile << " " << body->GetPos(); // } // myfile << std::endl; // myfile.close(); //} //else cout << "Unable to open file"; return 0; } private: std::vector<std::shared_ptr<ChBodyEasyBox>> bodies; // end bodies std::vector<std::shared_ptr<ChNodeFEAxyz>> nodes; // end nodes }; class Model6 { public: Model6(ChSystem& system, std::shared_ptr<ChMesh> mesh, int n_chains = 6, int elements = 10) : nodes(elements + 1) { // Global constants double PI = 3.14159; // Material inputs // Note for this simple test I am ignoring bending stiffness double chain_length = 10; // m double chain_mass_length = 1.1; // kg/m double chain_density = 7880; // kg/m^3 double chain_diameter = 13.332e-3; // m double chain_youngs_modulus = 53.7e9; // Pa // Material calculations double chain_cross_section_area = pow(chain_diameter, 2) * PI / 4; double element_length = chain_length / elements; double node_mass = chain_length / (elements + 1) * chain_mass_length; double element_stiffness = chain_youngs_modulus * chain_cross_section_area / element_length; auto mtruss = chrono_types::make_shared<ChBody>(); mtruss->SetBodyFixed(true); // Create the nodes for the backbone and link them together for (int j = 0; j < n_chains; ++j) { for (int i = 0; i < elements; ++i) { double x{ double(j * elements + i) }; double y{ 0 }; double z{ 0 }; std::cout << "Backbone: (" << x << ", " << y << ", " << z << ")\n"; int node_idx = j * elements + i; backbone_nodes.push_back(chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(x, y, z))); backbone_nodes[node_idx]->SetMass(node_mass); mesh->AddNode(backbone_nodes[node_idx]); // We need to fix the ends: if (j == 0 && i == 0 || (j == n_chains - 1 && i == elements - 1)) { auto constraint_hinge = chrono_types::make_shared<ChLinkPointFrame>(); constraint_hinge->Initialize(backbone_nodes[node_idx], mtruss); system.Add(constraint_hinge); // Visualize the hinge with a sphere auto msphere = chrono_types::make_shared<ChSphereShape>(); msphere->GetSphereGeometry().rad = 0.02; constraint_hinge->AddVisualShape(msphere); } // Connect all nodes to the prior node, except the first one (where there is no prior node): if (!(node_idx == 0)) { // Add a spring connecting the current node to the previous one auto melement = chrono_types::make_shared<ChElementSpring>(); melement->SetNodes(backbone_nodes[node_idx], backbone_nodes[node_idx - 1]); melement->SetSpringK(element_stiffness); mesh->AddElement(melement); } } } // Create the nodes for the hanging chains and link them together for (int j = 0; j < n_chains; ++j) { for (int i = 0; i < elements; ++i) { // We need to specify each nodes xyz in this formulation (inconvenient) // To make this easy, we say that for the purposes of the test lines // hang from the x-axis and nodes are rotated about the x-axis using the R_X // rotation matrix (see: https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations) double x{ double(j * elements) }; double y{ -1.0 * i * cos(0.523599) }; // 0.523599 == the angle of rotation from normal (radians) double z{ -1.0 * i * sin(0.523599) }; // This is really just SOHCAHTOA //std::cout << "(" << x << ", " << y << ", " << z << ")\n"; nodes.push_back(chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(x, y, z))); int node_idx = int(nodes.size()) - 1; int backbone_idx = j * elements; std::cout << backbone_idx << "\n"; //// Add a bar element to connect current node to previous one & add linear stiffness //auto melement = chrono_types::make_shared<ChElementBar>(); //melement->SetNodes(nodes[i], nodes[i-1]); //melement->SetBarArea(chain_cross_section_area); //melement->SetBarDensity(chain_density); //melement->SetBarYoungModulus(chain_youngs_modulus); //melement->SetBarRaleyghDamping(0); //mesh->AddElement(melement); // Add a spring connecting the current node to the previous one auto melement = chrono_types::make_shared<ChElementSpring>(); if (i == 0) { // Do nothing - we're going to use the zeroth node from the backbone } else if (i == 1) { // Set the first element to the appropriate backbone node nodes[node_idx]->SetMass(node_mass); mesh->AddNode(nodes[node_idx]); melement->SetNodes(nodes[node_idx], backbone_nodes[j * elements]); std::cout << "Connecting node at (" << nodes[node_idx]->GetPos().x() << ", " << nodes[node_idx]->GetPos().y() << ", " << nodes[node_idx]->GetPos().z(); std::cout << ") to (" << backbone_nodes[backbone_idx]->GetPos().x() << ", " << backbone_nodes[backbone_idx]->GetPos().y() << ", " << backbone_nodes[backbone_idx]->GetPos().z() << ")\n"; melement->SetSpringK(element_stiffness); mesh->AddElement(melement); } else { nodes[node_idx]->SetMass(node_mass); mesh->AddNode(nodes[node_idx]); melement->SetNodes(nodes[node_idx], nodes[node_idx - 1]); melement->SetSpringK(element_stiffness); mesh->AddElement(melement); } }; } } // Print positions of end bodies in each chain double PrintBodyPositions(ChSystem& system) { //std::cout << "Time: " << system.GetChTime() << std::endl; ofstream myfile("out.csv", std::ios_base::app); if (myfile.is_open()) { //myfile << system.GetChTime(); //myfile << clock() / CLOCKS_PER_SEC; for (auto node : nodes) { std::cout << node->GetPos().x() << std::endl; myfile << " " << node->GetPos(); } myfile << std::endl; myfile.close(); } else cout << "Unable to open file"; return bodies[0]->GetChTime(); return 0; } private: std::vector<std::shared_ptr<ChBodyEasyBox>> bodies; // end bodies std::vector<std::shared_ptr<ChNodeFEAxyz>> nodes; // hanging line nodes std::vector<std::shared_ptr<ChNodeFEAxyz>> backbone_nodes; // backbone nodes };
