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, Radu Serban
13 // =============================================================================
14 
15 #ifndef CHELEMENTBAR_H
16 #define CHELEMENTBAR_H
17 
18 #include "chrono/fea/ChElementGeneric.h"
19 #include "chrono/fea/ChNodeFEAxyz.h"
20 
21 namespace chrono {
22 namespace fea {
23 
24 /// @addtogroup fea_elements
25 /// @{
26 
27 /// Simple finite element with two nodes and a bar that connects them.
28 /// No bending and torsion stiffness, just like a bar with two spherical joints.
29 /// In practical terms, this element works a bit like the class ChElementSpring,
30 /// but also adds mass along the element, hence point-like mass in the two nodes
31 /// is not needed.
32 class ChApi ChElementBar : public ChElementGeneric {
33   public:
34     ChElementBar();
35     ~ChElementBar();
36 
GetNnodes()37     virtual int GetNnodes() override { return 2; }
GetNdofs()38     virtual int GetNdofs() override { return 2 * 3; }
GetNodeNdofs(int n)39     virtual int GetNodeNdofs(int n) override { return 3; }
40 
GetNodeN(int n)41     virtual std::shared_ptr<ChNodeFEAbase> GetNodeN(int n) override { return nodes[n]; }
42 
43     virtual void SetNodes(std::shared_ptr<ChNodeFEAxyz> nodeA, std::shared_ptr<ChNodeFEAxyz> nodeB);
44 
45     //
46     // FEM functions
47     //
48 
49     /// Fills the D vector with the current field values at the nodes of the element, with proper ordering.
50     /// If the D vector size is not this->GetNdofs(), it will be resized.
51     virtual void GetStateBlock(ChVectorDynamic<>& mD) override;
52 
53     /// Sets H as the global stiffness matrix K, scaled  by Kfactor. Optionally, also
54     /// superimposes global damping matrix R, scaled by Rfactor, and global mass matrix M multiplied by Mfactor.
55     /// (For the spring matrix there is no need to corotate local matrices: we already know a closed form expression.)
56     virtual void ComputeKRMmatricesGlobal(ChMatrixRef H,
57                                           double Kfactor,
58                                           double Rfactor = 0,
59                                           double Mfactor = 0) override;
60 
61     /// Computes the internal forces (ex. the actual position of nodes is not in relaxed reference position) and set
62     /// values in the Fi vector.
63     virtual void ComputeInternalForces(ChVectorDynamic<>& Fi) override;
64 
65     //
66     // Custom properties functions
67     //
68 
69     /// Set the cross sectional area of the bar (m^2) (also changes stiffness keeping same E modulus)
SetBarArea(double ma)70     void SetBarArea(double ma) { this->area = ma; }
GetBarArea()71     double GetBarArea() { return this->area; }
72 
73     /// Set the density of the bar (kg/m^3)
SetBarDensity(double md)74     void SetBarDensity(double md) { this->density = md; }
GetBarDensity()75     double GetBarDensity() { return this->density; }
76 
77     /// Set the Young elastic modulus (N/m^2) (also sets stiffness)
SetBarYoungModulus(double mE)78     void SetBarYoungModulus(double mE) { this->E = mE; }
GetBarYoungModulus()79     double GetBarYoungModulus() { return this->E; }
80 
81     /// Set the Rayleigh damping ratio r (as in: R = r * K )
SetBarRaleyghDamping(double mr)82     void SetBarRaleyghDamping(double mr) { this->rdamping = mr; }
GetBarRaleyghDamping()83     double GetBarRaleyghDamping() { return this->rdamping; }
84 
85     /// The full mass of the bar
GetMass()86     double GetMass() { return this->mass; }
87 
88     /// The rest length of the bar
GetRestLength()89     double GetRestLength() { return this->length; }
90 
91     /// The current length of the bar (might be after deformation)
GetCurrentLength()92     double GetCurrentLength() { return (nodes[1]->GetPos() - nodes[0]->GetPos()).Length(); }
93 
94     /// Get the strain (epsilon), after deformation.
GetStrain()95     double GetStrain() { return (GetCurrentLength() - GetRestLength()) / GetRestLength(); }
96 
97     /// Get the elastic stress (sigma), after deformation.
GetStress()98     double GetStress() { return GetBarYoungModulus() * GetStrain(); }
99 
100 	/// Get the current force transmitted along the bar direction,
101 	/// including the effect of the damper. Positive if pulled. (N)
102 	virtual double GetCurrentForce();
103 
104     //
105     // Functions for interfacing to the solver
106     //            (***not needed, thank to bookkeeping in parent class ChElementGeneric)
107 
108   private:
109     /// Initial setup. Precompute mass and matrices that do not change during the simulation, such as the local tangent
110     /// stiffness Kl of each element, if needed, etc.
111     virtual void SetupInitial(ChSystem* system) override;
112 
113     std::vector<std::shared_ptr<ChNodeFEAxyz> > nodes;
114     double area;
115     double density;
116     double E;
117     double rdamping;
118     double mass;
119     double length;
120 };
121 
122 /// @} fea_elements
123 
124 }  // end namespace fea
125 }  // end namespace chrono
126 
127 #endif
128