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, Alessandro Tasora 13 // ============================================================================= 14 15 #ifndef CH_ELEMENT_TETRA_COROT_4_H 16 #define CH_ELEMENT_TETRA_COROT_4_H 17 18 #include <cmath> 19 20 #include "chrono/fea/ChElementTetrahedron.h" 21 #include "chrono/fea/ChElementGeneric.h" 22 #include "chrono/fea/ChElementCorotational.h" 23 #include "chrono/fea/ChContinuumPoisson3D.h" 24 #include "chrono/fea/ChNodeFEAxyz.h" 25 #include "chrono/fea/ChNodeFEAxyzP.h" 26 #include "chrono/core/ChTensors.h" 27 28 namespace chrono { 29 namespace fea { 30 31 /// @addtogroup fea_elements 32 /// @{ 33 34 /// Tetrahedron FEA element with 4 nodes. 35 /// This is a classical element with linear displacement, hence with constant stress and constant strain. It can be 36 /// easily used for 3D FEA problems. 37 class ChApi ChElementTetraCorot_4 : public ChElementTetrahedron, 38 public ChElementGeneric, 39 public ChElementCorotational, 40 public ChLoadableUVW { 41 public: 42 using ShapeVector = ChMatrixNM<double, 1, 4>; 43 44 ChElementTetraCorot_4(); 45 ~ChElementTetraCorot_4(); 46 GetNnodes()47 virtual int GetNnodes() override { return 4; } GetNdofs()48 virtual int GetNdofs() override { return 4 * 3; } GetNodeNdofs(int n)49 virtual int GetNodeNdofs(int n) override { return 3; } 50 GetVolume()51 double GetVolume() { return Volume; } 52 GetNodeN(int n)53 virtual std::shared_ptr<ChNodeFEAbase> GetNodeN(int n) override { return nodes[n]; } 54 55 /// Return the specified tetrahedron node (0 <= n <= 3). GetTetrahedronNode(int n)56 virtual std::shared_ptr<ChNodeFEAxyz> GetTetrahedronNode(int n) override { return nodes[n]; } 57 58 virtual 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 63 // 64 // FEM functions 65 // 66 67 /// Update element at each time step. 68 virtual void Update() override; 69 70 /// Fills the N shape function matrix with the 71 /// values of shape functions at r,s,t 'volume' coordinates, where 72 /// r=1 at 2nd vertex, s=1 at 3rd, t = 1 at 4th. All ranging in [0...1]. 73 /// The last (u, u=1 at 1st vertex) is computed form the first 3 as 1.0-r-s-t. 74 /// NOTE! actually N should be a 3row, 12 column sparse matrix, 75 /// as N = [n1*eye(3) n2*eye(3) n3*eye(3) n4*eye(3)]; , 76 /// but to avoid wasting zero and repeated elements, here 77 /// it stores only the n1 n2 n3 n4 values in a 1 row, 4 columns matrix. 78 void ShapeFunctions(ShapeVector& N, double r, double s, double t); 79 80 /// Fills the D vector (displacement) with the currentfield values at the nodes of the element, with proper 81 /// ordering. If the D vector has not the size of this->GetNdofs(), it will be resized.For corotational elements, 82 /// field is assumed in local reference! 83 virtual void GetStateBlock(ChVectorDynamic<>& mD) override; 84 85 double ComputeVolume(); 86 87 /// Computes the local STIFFNESS MATRIX of the element: 88 /// K = Volume * [B]' * [D] * [B] 89 virtual void ComputeStiffnessMatrix(); 90 91 /// compute large rotation of element for corotational approach 92 virtual void UpdateRotation() override; 93 94 /// Sets H as the global stiffness matrix K, scaled by Kfactor. Optionally, also 95 /// superimposes global damping matrix R, scaled by Rfactor, and global mass matrix M multiplied by Mfactor. 96 virtual void ComputeKRMmatricesGlobal(ChMatrixRef H, 97 double Kfactor, 98 double Rfactor = 0, 99 double Mfactor = 0) override; 100 101 /// Computes the internal forces (ex. the actual position of nodes is not in relaxed reference position) and set 102 /// values in the Fi vector. 103 virtual void ComputeInternalForces(ChVectorDynamic<>& Fi) override; 104 105 // 106 // Custom properties functions 107 // 108 109 /// Set the material of the element SetMaterial(std::shared_ptr<ChContinuumElastic> my_material)110 void SetMaterial(std::shared_ptr<ChContinuumElastic> my_material) { Material = my_material; } GetMaterial()111 std::shared_ptr<ChContinuumElastic> GetMaterial() { return Material; } 112 113 /// Get the partial derivatives matrix MatrB and the StiffnessMatrix GetMatrB()114 const ChMatrixDynamic<>& GetMatrB() const { return MatrB; } GetStiffnessMatrix()115 const ChMatrixDynamic<>& GetStiffnessMatrix() const { return StiffnessMatrix; } 116 117 /// Returns the strain tensor (note that the tetrahedron 4 nodes is a linear 118 /// element, thus the strain is constant in the entire volume). 119 /// The tensor is in the original undeformed unrotated reference. 120 ChStrainTensor<> GetStrain(); 121 122 /// Returns the stress tensor (note that the tetrahedron 4 nodes is a linear 123 /// element, thus the stress is constant in the entire volume). 124 /// The tensor is in the original undeformed unrotated reference. 125 ChStressTensor<> GetStress(); 126 127 /// This function computes and adds corresponding masses to ElementBase member m_TotalMass 128 void ComputeNodalMass() override; 129 130 // 131 // Functions for interfacing to the solver 132 // (***not needed, thank to bookkeeping in parent class ChElementGeneric) 133 134 // 135 // Functions for ChLoadable interface 136 // 137 138 /// Gets the number of DOFs affected by this element (position part) LoadableGet_ndof_x()139 virtual int LoadableGet_ndof_x() override { return 4 * 3; } 140 141 /// Gets the number of DOFs affected by this element (speed part) LoadableGet_ndof_w()142 virtual int LoadableGet_ndof_w() override { return 4 * 3; } 143 144 /// Gets all the DOFs packed in a single vector (position part) 145 virtual void LoadableGetStateBlock_x(int block_offset, ChState& mD) override; 146 147 /// Gets all the DOFs packed in a single vector (speed part) 148 virtual void LoadableGetStateBlock_w(int block_offset, ChStateDelta& mD) override; 149 150 /// Increment all DOFs using a delta. 151 virtual void LoadableStateIncrement(const unsigned int off_x, 152 ChState& x_new, 153 const ChState& x, 154 const unsigned int off_v, 155 const ChStateDelta& Dv) override; 156 157 /// Number of coordinates in the interpolated field: here the {x,y,z} displacement Get_field_ncoords()158 virtual int Get_field_ncoords() override { return 3; } 159 160 /// Get the number of DOFs sub-blocks. GetSubBlocks()161 virtual int GetSubBlocks() override { return 4; } 162 163 /// Get the offset of the specified sub-block of DOFs in global vector. GetSubBlockOffset(int nblock)164 virtual unsigned int GetSubBlockOffset(int nblock) override { return nodes[nblock]->NodeGetOffset_w(); } 165 166 /// Get the size of the specified sub-block of DOFs in global vector. GetSubBlockSize(int nblock)167 virtual unsigned int GetSubBlockSize(int nblock) override { return 3; } 168 169 /// Check if the specified sub-block of DOFs is active. IsSubBlockActive(int nblock)170 virtual bool IsSubBlockActive(int nblock) const override { return !nodes[nblock]->GetFixed(); } 171 172 /// Get the pointers to the contained ChVariables, appending to the mvars vector. 173 virtual void LoadableGetVariables(std::vector<ChVariables*>& mvars) override; 174 175 /// Evaluate N'*F , where N is some type of shape function 176 /// evaluated at U,V,W coordinates of the volume, each ranging in 0..+1 177 /// F is a load, N'*F is the resulting generalized load 178 /// Returns also det[J] with J=[dx/du,..], that might be useful in gauss quadrature. 179 virtual void ComputeNF(const double U, ///< parametric coordinate in volume 180 const double V, ///< parametric coordinate in volume 181 const double W, ///< parametric coordinate in volume 182 ChVectorDynamic<>& Qi, ///< Return result of N'*F here, maybe with offset block_offset 183 double& detJ, ///< Return det[J] here 184 const ChVectorDynamic<>& F, ///< Input F vector, size is = n.field coords. 185 ChVectorDynamic<>* state_x, ///< if != 0, update state (pos. part) to this, then evaluate Q 186 ChVectorDynamic<>* state_w ///< if != 0, update state (speed part) to this, then evaluate Q 187 ) override; 188 189 /// This is needed so that it can be accessed by ChLoaderVolumeGravity GetDensity()190 virtual double GetDensity() override { return this->Material->Get_density(); } 191 192 /// If true, use quadrature over u,v,w in [0..1] range as tetrahedron volumetric coords, with z=1-u-v-w 193 /// otherwise use quadrature over u,v,w in [-1..+1] as box isoparametric coords. IsTetrahedronIntegrationNeeded()194 virtual bool IsTetrahedronIntegrationNeeded() override { return true; } 195 196 private: 197 /// Initial setup: set up the element's parameters and matrices 198 virtual void SetupInitial(ChSystem* system) override; 199 200 std::vector<std::shared_ptr<ChNodeFEAxyz> > nodes; 201 std::shared_ptr<ChContinuumElastic> Material; 202 ChMatrixDynamic<> MatrB; // matrix of shape function's partial derivatives 203 ChMatrixDynamic<> StiffnessMatrix; // undeformed local stiffness matrix 204 ChMatrixNM<double, 4, 4> mM; // for speeding up corotational approach 205 double Volume; 206 207 public: 208 EIGEN_MAKE_ALIGNED_OPERATOR_NEW 209 }; 210 211 // ----------------------------------------------------------------------------- 212 213 /// Tetrahedron FEM element with 4 nodes for scalar fields (for Poisson-like problems). 214 /// This is a classical element with linear displacement. 215 /// ***EXPERIMENTAL*** 216 class ChApi ChElementTetraCorot_4_P : public ChElementGeneric, public ChElementCorotational, public ChLoadableUVW { 217 public: 218 using ShapeVector = ChMatrixNM<double, 1, 4>; 219 220 ChElementTetraCorot_4_P(); ~ChElementTetraCorot_4_P()221 ~ChElementTetraCorot_4_P() {} 222 GetNnodes()223 virtual int GetNnodes() override { return 4; } GetNdofs()224 virtual int GetNdofs() override { return 4 * 1; } GetNodeNdofs(int n)225 virtual int GetNodeNdofs(int n) override { return 1; } 226 GetVolume()227 double GetVolume() { return Volume; } 228 GetNodeN(int n)229 virtual std::shared_ptr<ChNodeFEAbase> GetNodeN(int n) override { return nodes[n]; } 230 231 virtual void SetNodes(std::shared_ptr<ChNodeFEAxyzP> nodeA, 232 std::shared_ptr<ChNodeFEAxyzP> nodeB, 233 std::shared_ptr<ChNodeFEAxyzP> nodeC, 234 std::shared_ptr<ChNodeFEAxyzP> nodeD); 235 236 // 237 // FEM functions 238 // 239 240 /// Update element at each time step. 241 virtual void Update() override; 242 243 /// Fills the N shape function matrix with the 244 /// values of shape functions at zi parametric coordinates, where 245 /// z0=1 at 1st vertex, z1=1 at second, z2 = 1 at third (volumetric shape functions). 246 /// The 4th is computed form the first 3. All ranging in [0...1]. 247 /// NOTE! actually N should be a 3row, 12 column sparse matrix, 248 /// as N = [n1*eye(3) n2*eye(3) n3*eye(3) n4*eye(3)]; , 249 /// but to avoid wasting zero and repeated elements, here 250 /// it stores only the n1 n2 n3 n4 values in a 1 row, 4 columns matrix! 251 virtual void ShapeFunctions(ShapeVector& N, double z0, double z1, double z2); 252 253 /// Fills the D vector with the current field values at the nodes of the element, with proper ordering. If the D 254 /// vector has not the size of this->GetNdofs(), it will be resized. For corotational elements, field is assumed in 255 /// local reference! 256 virtual void GetStateBlock(ChVectorDynamic<>& mD) override; 257 258 double ComputeVolume(); 259 260 /// Computes the local STIFFNESS MATRIX of the element: 261 /// K = Volume * [B]' * [D] * [B] 262 virtual void ComputeStiffnessMatrix(); 263 264 // compute large rotation of element for corotational approach 265 // Btw: NOT really needed for Poisson problems 266 virtual void UpdateRotation() override; 267 268 /// Sets H as the global stiffness matrix K, scaled by Kfactor. Optionally, also 269 /// superimposes global damping matrix R, scaled by Rfactor, and global mass matrix M multiplied by Mfactor. 270 virtual void ComputeKRMmatricesGlobal(ChMatrixRef H, 271 double Kfactor, 272 double Rfactor = 0, 273 double Mfactor = 0) override; 274 275 /// Computes the internal 'pseudo-forces' and set values in the Fi vector. 276 virtual void ComputeInternalForces(ChVectorDynamic<>& Fi) override; 277 278 // 279 // Custom properties functions 280 // 281 282 /// Set the material of the element SetMaterial(std::shared_ptr<ChContinuumPoisson3D> my_material)283 void SetMaterial(std::shared_ptr<ChContinuumPoisson3D> my_material) { Material = my_material; } GetMaterial()284 std::shared_ptr<ChContinuumPoisson3D> GetMaterial() { return Material; } 285 286 /// Get the partial derivatives matrix MatrB and the StiffnessMatrix GetMatrB()287 const ChMatrixDynamic<>& GetMatrB() const { return MatrB; } GetStiffnessMatrix()288 const ChMatrixDynamic<>& GetStiffnessMatrix() const { return StiffnessMatrix; } 289 290 /// Returns the gradient of P (note that the tetrahedron 4 nodes is a linear 291 /// element, thus the gradient is constant in the entire volume). 292 /// It is in the original undeformed unrotated reference. 293 ChVectorN<double, 3> GetPgradient(); 294 295 // 296 // Functions for interfacing to the solver 297 // (***not needed, thank to bookkeeping in parent class ChElementGeneric) 298 299 // 300 // Functions for ChLoadable interface 301 // 302 303 /// Gets the number of DOFs affected by this element (position part) LoadableGet_ndof_x()304 virtual int LoadableGet_ndof_x() override { return 4 * 3; } 305 306 /// Gets the number of DOFs affected by this element (speed part) LoadableGet_ndof_w()307 virtual int LoadableGet_ndof_w() override { return 4 * 3; } 308 309 /// Gets all the DOFs packed in a single vector (position part) 310 virtual void LoadableGetStateBlock_x(int block_offset, ChState& mD) override; 311 312 /// Gets all the DOFs packed in a single vector (speed part) 313 virtual void LoadableGetStateBlock_w(int block_offset, ChStateDelta& mD) override; 314 315 /// Increment all DOFs using a delta. 316 virtual void LoadableStateIncrement(const unsigned int off_x, 317 ChState& x_new, 318 const ChState& x, 319 const unsigned int off_v, 320 const ChStateDelta& Dv) override; 321 322 /// Number of coordinates in the interpolated field: here the {t} temperature Get_field_ncoords()323 virtual int Get_field_ncoords() override { return 1; } 324 325 /// Get the number of DOFs sub-blocks. GetSubBlocks()326 virtual int GetSubBlocks() override { return 4; } 327 328 /// Get the offset of the specified sub-block of DOFs in global vector. GetSubBlockOffset(int nblock)329 virtual unsigned int GetSubBlockOffset(int nblock) override { return nodes[nblock]->NodeGetOffset_w(); } 330 331 /// Get the size of the specified sub-block of DOFs in global vector. GetSubBlockSize(int nblock)332 virtual unsigned int GetSubBlockSize(int nblock) override { return 1; } 333 334 /// Check if the specified sub-block of DOFs is active. IsSubBlockActive(int nblock)335 virtual bool IsSubBlockActive(int nblock) const override { return !nodes[nblock]->GetFixed(); } 336 337 /// Get the pointers to the contained ChVariables, appending to the mvars vector. 338 virtual void LoadableGetVariables(std::vector<ChVariables*>& mvars) override; 339 340 /// Evaluate N'*F , where N is some type of shape function 341 /// evaluated at U,V,W coordinates of the volume, each ranging in 0..+1 342 /// F is a load, N'*F is the resulting generalized load 343 /// Returns also det[J] with J=[dx/du,..], that might be useful in gauss quadrature. 344 virtual void ComputeNF(const double U, ///< parametric coordinate in volume 345 const double V, ///< parametric coordinate in volume 346 const double W, ///< parametric coordinate in volume 347 ChVectorDynamic<>& Qi, ///< Return result of N'*F here, maybe with offset block_offset 348 double& detJ, ///< Return det[J] here 349 const ChVectorDynamic<>& F, ///< Input F vector, size is = n.field coords. 350 ChVectorDynamic<>* state_x, ///< if != 0, update state (pos. part) to this, then evaluate Q 351 ChVectorDynamic<>* state_w ///< if != 0, update state (speed part) to this, then evaluate Q 352 ) override; 353 354 /// Return 0 if not supportable by ChLoaderVolumeGravity GetDensity()355 virtual double GetDensity() override { return 0; } 356 357 /// If true, use quadrature over u,v,w in [0..1] range as tetrahedron volumetric coords, with z=1-u-v-w 358 /// otherwise use quadrature over u,v,w in [-1..+1] as box isoparametric coords. IsTetrahedronIntegrationNeeded()359 virtual bool IsTetrahedronIntegrationNeeded() override { return true; } 360 361 private: 362 /// Initial setup: set up the element's parameters and matrices 363 virtual void SetupInitial(ChSystem* system) override; 364 365 std::vector<std::shared_ptr<ChNodeFEAxyzP> > nodes; 366 std::shared_ptr<ChContinuumPoisson3D> Material; 367 ChMatrixDynamic<> MatrB; // matrix of shape function's partial derivatives 368 ChMatrixDynamic<> StiffnessMatrix; // local stiffness matrix 369 ChMatrixNM<double, 4, 4> mM; // for speeding up corotational approach 370 double Volume; 371 372 public: 373 EIGEN_MAKE_ALIGNED_OPERATOR_NEW 374 }; 375 376 /// @} fea_elements 377 378 } // end namespace fea 379 } // end namespace chrono 380 381 #endif 382