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 // Authors: Bryan Peterson, Antonio Recuero, Radu Serban 12 // ============================================================================= 13 // Hexahedronal element with 8 nodes (with EAS) 14 // ============================================================================= 15 16 #ifndef CH_ELEMENT_HEXA_ANCF_3813_H 17 #define CH_ELEMENT_HEXA_ANCF_3813_H 18 19 #include "chrono/fea/ChElementHexahedron.h" 20 #include "chrono/fea/ChElementGeneric.h" 21 #include "chrono/fea/ChNodeFEAxyz.h" 22 #include "chrono/fea/ChContinuumMaterial.h" 23 #include "chrono/core/ChQuadrature.h" 24 #include "chrono/physics/ChLoadable.h" 25 26 namespace chrono { 27 namespace fea { 28 29 /// @addtogroup fea_elements 30 /// @{ 31 32 /// Hexahedronal solid element with 8 nodes (with EAS). 33 /// While technically not an ANCF element, the name is justified because the implementation can use the same ANCF 34 /// machinery. 35 class ChApi ChElementHexaANCF_3813 : public ChElementHexahedron, public ChElementGeneric, public ChLoadableUVW { 36 public: 37 using ShapeVector = ChMatrixNM<double, 1, 8>; 38 39 ChElementHexaANCF_3813(); ~ChElementHexaANCF_3813()40 ~ChElementHexaANCF_3813() {} 41 42 /// Get number of nodes of this element GetNnodes()43 virtual int GetNnodes() override { return 8; } 44 45 /// Get number of degrees of freedom of this element GetNdofs()46 virtual int GetNdofs() override { return 8 * 3; } 47 48 /// Get the number of coordinates from the n-th node used by this element. GetNodeNdofs(int n)49 virtual int GetNodeNdofs(int n) override { return 3; } 50 51 /// Access the n-th node of this element. GetNodeN(int n)52 virtual std::shared_ptr<ChNodeFEAbase> GetNodeN(int n) override { return m_nodes[n]; } 53 54 /// Return the specified hexahedron node (0 <= n <= 7). GetHexahedronNode(int n)55 virtual std::shared_ptr<ChNodeFEAxyz> GetHexahedronNode(int n) override { return m_nodes[n]; } 56 57 /// Specify the nodes of this element. 58 void SetNodes(std::shared_ptr<ChNodeFEAxyz> nodeA, 59 std::shared_ptr<ChNodeFEAxyz> nodeB, 60 std::shared_ptr<ChNodeFEAxyz> nodeC, 61 std::shared_ptr<ChNodeFEAxyz> nodeD, 62 std::shared_ptr<ChNodeFEAxyz> nodeE, 63 std::shared_ptr<ChNodeFEAxyz> nodeF, 64 std::shared_ptr<ChNodeFEAxyz> nodeG, 65 std::shared_ptr<ChNodeFEAxyz> nodeH); 66 67 /// Set the element number. SetElemNum(int kb)68 void SetElemNum(int kb) { m_elementnumber = kb; } 69 70 /// Set EAS internal parameters (stored values). 71 void 72 SetStockAlpha(double a1, double a2, double a3, double a4, double a5, double a6, double a7, double a8, double a9); 73 74 /// Set some element parameters (dimensions). SetInertFlexVec(const ChVector<double> & a)75 void SetInertFlexVec(const ChVector<double>& a) { m_InertFlexVec = a; } 76 GetElemNum()77 int GetElemNum() const { return m_elementnumber; } 78 79 /// Get initial position of the element in matrix form GetInitialPos()80 const ChMatrixNM<double, 8, 3>& GetInitialPos() const { return m_d0; } 81 GetNodeA()82 std::shared_ptr<ChNodeFEAxyz> GetNodeA() const { return m_nodes[0]; } GetNodeB()83 std::shared_ptr<ChNodeFEAxyz> GetNodeB() const { return m_nodes[1]; } GetNodeC()84 std::shared_ptr<ChNodeFEAxyz> GetNodeC() const { return m_nodes[2]; } GetNodeD()85 std::shared_ptr<ChNodeFEAxyz> GetNodeD() const { return m_nodes[3]; } GetNodeE()86 std::shared_ptr<ChNodeFEAxyz> GetNodeE() const { return m_nodes[4]; } GetNodeF()87 std::shared_ptr<ChNodeFEAxyz> GetNodeF() const { return m_nodes[5]; } GetNodeG()88 std::shared_ptr<ChNodeFEAxyz> GetNodeG() const { return m_nodes[6]; } GetNodeH()89 std::shared_ptr<ChNodeFEAxyz> GetNodeH() const { return m_nodes[7]; } 90 GetLengthX()91 double GetLengthX() const { return m_InertFlexVec.x(); } GetLengthY()92 double GetLengthY() const { return m_InertFlexVec.y(); } GetLengthZ()93 double GetLengthZ() const { return m_InertFlexVec.z(); } 94 SetMaterial(std::shared_ptr<ChContinuumElastic> my_material)95 void SetMaterial(std::shared_ptr<ChContinuumElastic> my_material) { m_Material = my_material; } GetMaterial()96 std::shared_ptr<ChContinuumElastic> GetMaterial() const { return m_Material; } 97 /// Set whether material is Mooney-Rivlin (Otherwise linear elastic isotropic) SetMooneyRivlin(bool val)98 void SetMooneyRivlin(bool val) { m_isMooney = val; } 99 /// Set Mooney-Rivlin coefficients SetMRCoefficients(double C1,double C2)100 void SetMRCoefficients(double C1, double C2) { 101 CCOM1 = C1; 102 CCOM2 = C2; 103 } 104 /// Fills the N shape function matrix 105 /// as N = [s1*eye(3) s2*eye(3) s3*eye(3) s4*eye(3)...]; , 106 void ShapeFunctions(ShapeVector& N, double x, double y, double z); 107 void ShapeFunctionsDerivativeX(ShapeVector& Nx, double x, double y, double z); 108 void ShapeFunctionsDerivativeY(ShapeVector& Ny, double x, double y, double z); 109 void ShapeFunctionsDerivativeZ(ShapeVector& Nz, double x, double y, double z); 110 // Functions for ChLoadable interface 111 /// Gets the number of DOFs affected by this element (position part) LoadableGet_ndof_x()112 virtual int LoadableGet_ndof_x() override { return 8 * 3; } 113 114 /// Gets the number of DOFs affected by this element (speed part) LoadableGet_ndof_w()115 virtual int LoadableGet_ndof_w() override { return 8 * 3; } 116 117 /// Gets all the DOFs packed in a single vector (position part) 118 virtual void LoadableGetStateBlock_x(int block_offset, ChState& mD) override; 119 120 /// Gets all the DOFs packed in a single vector (speed part) 121 virtual void LoadableGetStateBlock_w(int block_offset, ChStateDelta& mD) override; 122 123 /// Increment all DOFs using a delta. 124 virtual void LoadableStateIncrement(const unsigned int off_x, 125 ChState& x_new, 126 const ChState& x, 127 const unsigned int off_v, 128 const ChStateDelta& Dv) override; 129 130 /// Number of coordinates in the interpolated field: here the {x,y,z} displacement Get_field_ncoords()131 virtual int Get_field_ncoords() override { return 3; } 132 133 /// Get the number of DOFs sub-blocks. GetSubBlocks()134 virtual int GetSubBlocks() override { return 8; } 135 136 /// Get the offset of the specified sub-block of DOFs in global vector. GetSubBlockOffset(int nblock)137 virtual unsigned int GetSubBlockOffset(int nblock) override { return m_nodes[nblock]->NodeGetOffset_w(); } 138 139 /// Get the size of the specified sub-block of DOFs in global vector. GetSubBlockSize(int nblock)140 virtual unsigned int GetSubBlockSize(int nblock) override { return 3; } 141 142 /// Check if the specified sub-block of DOFs is active. IsSubBlockActive(int nblock)143 virtual bool IsSubBlockActive(int nblock) const override { return !m_nodes[nblock]->GetFixed(); } 144 145 /// Get the pointers to the contained ChVariables, appending to the mvars vector. 146 virtual void LoadableGetVariables(std::vector<ChVariables*>& mvars) override; 147 148 /// Evaluate N'*F , where N is some type of shape function 149 /// evaluated at U,V,W coordinates of the volume, each ranging in -1..+1 150 /// F is a load, N'*F is the resulting generalized load 151 /// Returns also det[J] with J=[dx/du,..], that might be useful in gauss quadrature. 152 virtual void ComputeNF(const double U, ///< parametric coordinate in volume 153 const double V, ///< parametric coordinate in volume 154 const double W, ///< parametric coordinate in volume 155 ChVectorDynamic<>& Qi, ///< Return result of N'*F here, maybe with offset block_offset 156 double& detJ, ///< Return det[J] here 157 const ChVectorDynamic<>& F, ///< Input F vector, size is = n.field coords. 158 ChVectorDynamic<>* state_x, ///< if != 0, update state (pos. part) to this, then evaluate Q 159 ChVectorDynamic<>* state_w ///< if != 0, update state (speed part) to this, then evaluate Q 160 ) override; 161 162 /// This is needed so that it can be accessed by ChLoaderVolumeGravity GetDensity()163 virtual double GetDensity() override { return this->m_Material->Get_density(); } 164 165 private: 166 enum JacobianType { ANALYTICAL, NUMERICAL }; 167 168 // Private Data 169 std::vector<std::shared_ptr<ChNodeFEAxyz> > m_nodes; ///< Element nodes 170 171 std::shared_ptr<ChContinuumElastic> m_Material; ///< Elastic Material 172 173 ChMatrixNM<double, 24, 24> m_StiffnessMatrix; ///< Stiffness matrix 174 ChMatrixNM<double, 24, 24> m_MassMatrix; ///< Mass matrix 175 ChVectorN<double, 8> m_GravForceScale; ///< Gravity scaling matrix used to get the generalized force due to gravity 176 ChVector<double> m_InertFlexVec; ///< for element size (EL,EW,EH) 177 // EAS 178 int m_elementnumber; ///< Element number, for EAS 179 ChMatrixNM<double, 24, 24> m_stock_jac_EAS; ///< EAS Jacobian matrix 180 ChVectorN<double, 9> m_stock_alpha_EAS; ///< EAS previous step internal parameters 181 ChMatrixNM<double, 24, 24> m_stock_KTE; ///< Analytical Jacobian 182 ChMatrixNM<double, 8, 3> m_d0; ///< Initial Coordinate per element 183 JacobianType m_flag_HE; 184 bool m_isMooney; ///< Flag indicating whether the material is Mooney Rivlin 185 double CCOM1; ///< First coefficient for Mooney-Rivlin 186 double CCOM2; ///< Second coefficient for Mooney-Rivlin 187 // Private Methods 188 189 virtual void Update() override; 190 191 /// Fills the D vector with the current field values at the nodes of the element, with proper ordering. 192 /// If the D vector has not the size of this->GetNdofs(), it will be resized. 193 /// {x_a y_a z_a Dx_a Dx_a Dx_a x_b y_b z_b Dx_b Dy_b Dz_b} 194 virtual void GetStateBlock(ChVectorDynamic<>& mD) override; 195 196 /// Computes the STIFFNESS MATRIX of the element: 197 /// K = integral( .... ), 198 /// Note: in this 'basic' implementation, constant section and 199 /// constant material are assumed 200 void ComputeStiffnessMatrix(); 201 /// Computes the MASS MATRIX of the element 202 /// Note: in this 'basic' implementation, constant section and 203 /// constant material are assumed 204 void ComputeMassMatrix(); 205 /// Compute the matrix to scale gravity by to get the generalized gravitational force. 206 void ComputeGravityForceScale(); 207 /// Initial setup. Precompute mass and matrices that do not change during the simulation. 208 virtual void SetupInitial(ChSystem* system) override; 209 /// Sets M as the global mass matrix. ComputeMmatrixGlobal(ChMatrixRef M)210 virtual void ComputeMmatrixGlobal(ChMatrixRef M) override { M = m_MassMatrix; } 211 212 /// Compute the generalized force vector due to gravity using the efficient element specific method 213 virtual void ComputeGravityForces(ChVectorDynamic<>& Fg, const ChVector<>& G_acc) override; 214 215 /// Sets H as the global stiffness matrix K, scaled by Kfactor. Optionally, also 216 /// superimposes global damping matrix R, scaled by Rfactor, and global mass matrix M multiplied by Mfactor. 217 virtual void ComputeKRMmatricesGlobal(ChMatrixRef H, 218 double Kfactor, 219 double Rfactor = 0, 220 double Mfactor = 0) override; 221 222 /// Computes the internal forces (ex. the actual position of nodes is not in relaxed reference position) and set 223 /// values in the Fi vector. 224 virtual void ComputeInternalForces(ChVectorDynamic<>& Fi) override; 225 226 // [EAS] matrix T0 (inverse and transposed) and detJ0 at center are used for Enhanced Assumed Strains alpha 227 void T0DetJElementCenterForEAS(ChMatrixNM<double, 8, 3>& d0, ChMatrixNM<double, 6, 6>& T0, double& detJ0C); 228 // [EAS] Basis function of M for Enhanced Assumed Strain 229 void Basis_M(ChMatrixNM<double, 6, 9>& M, double x, double y, double z); 230 231 friend class Brick_ForceAnalytical; 232 friend class Brick_ForceNumerical; 233 234 public: 235 EIGEN_MAKE_ALIGNED_OPERATOR_NEW 236 }; 237 238 /// @} fea_elements 239 240 } // end namespace fea 241 } // end namespace chrono 242 243 #endif 244