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: Andrea Favali, Radu Serban 13 // ============================================================================= 14 15 #ifndef CH_ELEMENT_HEXA_COROT_8_H 16 #define CH_ELEMENT_HEXA_COROT_8_H 17 18 #include "chrono/fea/ChElementHexahedron.h" 19 #include "chrono/fea/ChElementGeneric.h" 20 #include "chrono/fea/ChElementCorotational.h" 21 #include "chrono/fea/ChGaussIntegrationRule.h" 22 #include "chrono/fea/ChNodeFEAxyz.h" 23 24 namespace chrono { 25 namespace fea { 26 27 /// @addtogroup fea_elements 28 /// @{ 29 30 /// Class for FEA elements of hexahedron type (isoparametric 3D bricks) with 8 nodes. 31 /// This element has a linear displacement field. 32 class ChApi ChElementHexaCorot_8 : public ChElementHexahedron, 33 public ChElementGeneric, 34 public ChElementCorotational, 35 public ChLoadableUVW { 36 public: 37 using ShapeVector = ChMatrixNM<double, 1, 8>; 38 39 ChElementHexaCorot_8(); 40 ~ChElementHexaCorot_8(); 41 GetNnodes()42 virtual int GetNnodes() override { return 8; } GetNdofs()43 virtual int GetNdofs() override { return 8 * 3; } GetNodeNdofs(int n)44 virtual int GetNodeNdofs(int n) override { return 3; } 45 GetVolume()46 double GetVolume() { return Volume; } 47 GetNodeN(int n)48 virtual std::shared_ptr<ChNodeFEAbase> GetNodeN(int n) override { return nodes[n]; } 49 50 /// Return the specified hexahedron node (0 <= n <= 7). GetHexahedronNode(int n)51 virtual std::shared_ptr<ChNodeFEAxyz> GetHexahedronNode(int n) override { return nodes[n]; } 52 53 virtual void SetNodes(std::shared_ptr<ChNodeFEAxyz> nodeA, 54 std::shared_ptr<ChNodeFEAxyz> nodeB, 55 std::shared_ptr<ChNodeFEAxyz> nodeC, 56 std::shared_ptr<ChNodeFEAxyz> nodeD, 57 std::shared_ptr<ChNodeFEAxyz> nodeE, 58 std::shared_ptr<ChNodeFEAxyz> nodeF, 59 std::shared_ptr<ChNodeFEAxyz> nodeG, 60 std::shared_ptr<ChNodeFEAxyz> nodeH); 61 62 // 63 // QUADRATURE functions 64 // 65 SetDefaultIntegrationRule()66 virtual void SetDefaultIntegrationRule() { this->ir->SetIntOnCube(8, &this->GpVector); } 67 SetReducedIntegrationRule()68 virtual void SetReducedIntegrationRule() { this->ir->SetIntOnCube(1, &this->GpVector); } 69 SetIntegrationRule(int nPoints)70 virtual void SetIntegrationRule(int nPoints) { this->ir->SetIntOnCube(nPoints, &this->GpVector); } 71 72 // 73 // FEA functions 74 // 75 76 /// Fills the N shape function matrix with the 77 /// values of shape functions at z0,z1,z2 parametric coordinates, where 78 /// each zi is in [-1...+1] range. 79 /// It stores the Ni(z0,z1,z2) values in a 1 row, 8 columns matrix. 80 void ShapeFunctions(ShapeVector& N, double z0, double z1, double z2); 81 82 /// Fills the D vector (displacement) with the current field values at the nodes of the element, with proper 83 /// ordering. If the D vector has not the size of this->GetNdofs(), it will be resized. For corotational elements, 84 /// field is assumed in local reference! 85 virtual void GetStateBlock(ChVectorDynamic<>& mD) override; 86 87 /// Puts inside 'Jacobian' and 'J1' the Jacobian matrix and the shape functions derivatives matrix of the element. 88 /// The vector "coord" contains the natural coordinates of the integration point. 89 /// in case of hexahedral elements natural coords vary in the classical range -1 ... +1. 90 virtual void ComputeJacobian(ChMatrixDynamic<>& Jacobian, ChMatrixDynamic<>& J1, ChVector<> coord); 91 92 /// Computes the matrix of partial derivatives and puts data in "MatrB" 93 /// evaluated at natural coordinates zeta1,...,zeta4 . Also computes determinant of jacobian. 94 /// note: in case of hexahedral elements natural coord. vary in the range -1 ... +1 95 virtual void ComputeMatrB(ChMatrixDynamic<>& MatrB, double zeta1, double zeta2, double zeta3, double& JacobianDet); 96 97 /// Computes the matrix of partial derivatives and puts data in "GaussPt". 98 /// Stores the determinant of the jacobian in "JacobianDet". 99 virtual void ComputeMatrB(ChGaussPoint* GaussPt, double& JacobianDet); 100 101 /// Computes the global STIFFNESS MATRIX of the element: 102 /// K = sum (w_i * [B]' * [D] * [B]) 103 /// The number of Gauss Point is defined by SetIntegrationRule function (default: 8 Gp). 104 virtual void ComputeStiffnessMatrix(); 105 106 /// Update element at each time step. 107 virtual void Update() override; 108 109 // Compute large rotation of element for corotational approach 110 virtual void UpdateRotation() override; 111 112 /// Returns the strain tensor at given parameters. 113 /// The tensor is in the original undeformed unrotated reference. 114 ChStrainTensor<> GetStrain(double z1, double z2, double z3); 115 116 /// Returns the stress tensor at given parameters. 117 /// The tensor is in the original undeformed unrotated reference. 118 ChStressTensor<> GetStress(double z1, double z2, double z3); 119 120 /// Sets H as the global stiffness matrix K, scaled by Kfactor. Optionally, also 121 /// superimposes global damping matrix R, scaled by Rfactor, and global mass matrix M multiplied by Mfactor. 122 virtual void ComputeKRMmatricesGlobal(ChMatrixRef H, 123 double Kfactor, 124 double Rfactor = 0, 125 double Mfactor = 0) override; 126 127 /// Computes the internal forces (ex. the actual position of nodes is not in relaxed reference position) and set 128 /// values in the Fi vector. 129 virtual void ComputeInternalForces(ChVectorDynamic<>& Fi) override; 130 131 // 132 // Custom properties functions 133 // 134 135 /// Set the material of the element SetMaterial(std::shared_ptr<ChContinuumElastic> my_material)136 void SetMaterial(std::shared_ptr<ChContinuumElastic> my_material) { Material = my_material; } GetMaterial()137 std::shared_ptr<ChContinuumElastic> GetMaterial() { return Material; } 138 139 /// Get the StiffnessMatrix GetStiffnessMatrix()140 const ChMatrixDynamic<>& GetStiffnessMatrix() const { return StiffnessMatrix; } 141 142 /// Get the Nth gauss point GetGaussPoint(int N)143 ChGaussPoint* GetGaussPoint(int N) { return GpVector[N]; } 144 145 // 146 // Functions for interfacing to the solver 147 // (***not needed, thank to bookkeeping in parent class ChElementGeneric) 148 149 // 150 // Functions for ChLoadable interface 151 // 152 153 /// Gets the number of DOFs affected by this element (position part) LoadableGet_ndof_x()154 virtual int LoadableGet_ndof_x() override { return 8 * 3; } 155 156 /// Gets the number of DOFs affected by this element (speed part) LoadableGet_ndof_w()157 virtual int LoadableGet_ndof_w() override { return 8 * 3; } 158 159 /// Gets all the DOFs packed in a single vector (position part) 160 virtual void LoadableGetStateBlock_x(int block_offset, ChState& mD) override; 161 162 /// Gets all the DOFs packed in a single vector (speed part) 163 virtual void LoadableGetStateBlock_w(int block_offset, ChStateDelta& mD) override; 164 165 /// Increment all DOFs using a delta. 166 virtual void LoadableStateIncrement(const unsigned int off_x, 167 ChState& x_new, 168 const ChState& x, 169 const unsigned int off_v, 170 const ChStateDelta& Dv) override; 171 172 /// Number of coordinates in the interpolated field: here the {x,y,z} displacement Get_field_ncoords()173 virtual int Get_field_ncoords() override { return 3; } 174 175 /// Get the number of DOFs sub-blocks. GetSubBlocks()176 virtual int GetSubBlocks() override { return 8; } 177 178 /// Get the offset of the specified sub-block of DOFs in global vector. GetSubBlockOffset(int nblock)179 virtual unsigned int GetSubBlockOffset(int nblock) override { return nodes[nblock]->NodeGetOffset_w(); } 180 181 /// Get the size of the specified sub-block of DOFs in global vector. GetSubBlockSize(int nblock)182 virtual unsigned int GetSubBlockSize(int nblock) override { return 3; } 183 184 /// Check if the specified sub-block of DOFs is active. IsSubBlockActive(int nblock)185 virtual bool IsSubBlockActive(int nblock) const override { return !nodes[nblock]->GetFixed(); } 186 187 /// Get the pointers to the contained ChVariables, appending to the mvars vector. 188 virtual void LoadableGetVariables(std::vector<ChVariables*>& mvars) override; 189 190 /// Evaluate N'*F , where N is some type of shape function 191 /// evaluated at U,V,W coordinates of the volume, each ranging in -1..+1 192 /// F is a load, N'*F is the resulting generalized load 193 /// Returns also det[J] with J=[dx/du,..], that might be useful in gauss quadrature. 194 virtual void ComputeNF(const double U, ///< parametric coordinate in volume 195 const double V, ///< parametric coordinate in volume 196 const double W, ///< parametric coordinate in volume 197 ChVectorDynamic<>& Qi, ///< Return result of N'*F here, maybe with offset block_offset 198 double& detJ, ///< Return det[J] here 199 const ChVectorDynamic<>& F, ///< Input F vector, size is = n.field coords. 200 ChVectorDynamic<>* state_x, ///< if != 0, update state (pos. part) to this, then evaluate Q 201 ChVectorDynamic<>* state_w ///< if != 0, update state (speed part) to this, then evaluate Q 202 ) override; 203 204 /// This is needed so that it can be accessed by ChLoaderVolumeGravity GetDensity()205 virtual double GetDensity() override { return this->Material->Get_density(); } 206 207 private: SetupInitial(ChSystem * system)208 virtual void SetupInitial(ChSystem* system) override { ComputeStiffnessMatrix(); } 209 210 std::vector<std::shared_ptr<ChNodeFEAxyz> > nodes; 211 std::shared_ptr<ChContinuumElastic> Material; 212 ChMatrixDynamic<> StiffnessMatrix; 213 214 ChGaussIntegrationRule* ir; 215 std::vector<ChGaussPoint*> GpVector; 216 double Volume; 217 }; 218 219 /// @} fea_elements 220 221 } // end namespace fea 222 } // end namespace chrono 223 224 #endif 225