Public Member Functions | Protected Member Functions | Protected Attributes | Private Types | Private Member Functions | List of all members
SurgSim::Physics::Fem2DElementTriangle Class Reference

2D FemElement based on a triangle with a constant thickness More...

#include <SurgSim/Physics/Fem2DElementTriangle.h>

Inheritance diagram for SurgSim::Physics::Fem2DElementTriangle:
SurgSim::Physics::FemElement

Public Member Functions

 Fem2DElementTriangle ()
 Constructor. More...
 
 Fem2DElementTriangle (std::array< size_t, 3 > nodeIds)
 Constructor. More...
 
 Fem2DElementTriangle (std::shared_ptr< FemElementStructs::FemElementParameter > elementData)
 Constructor for FemElement object factory. More...
 
void setThickness (double thickness)
 Sets the triangle's thickness. More...
 
double getThickness () const
 Gets the triangle's thickness. More...
 
void initialize (const SurgSim::Math::OdeState &state) override
 Initialize the FemElement once everything has been set. More...
 
double getVolume (const SurgSim::Math::OdeState &state) const override
 Gets the element volume based on the input state (in m-3) More...
 
SurgSim::Math::Vector computeCartesianCoordinate (const SurgSim::Math::OdeState &state, const SurgSim::Math::Vector &naturalCoordinate) const override
 Computes a given natural coordinate in cartesian coordinates. More...
 
SurgSim::Math::Vector computeNaturalCoordinate (const SurgSim::Math::OdeState &state, const SurgSim::Math::Vector &cartesianCoordinate) const override
 Computes a natural coordinate given a global coordinate. More...
 
- Public Member Functions inherited from SurgSim::Physics::FemElement
 FemElement ()
 Constructor. More...
 
virtual ~FemElement ()
 Virtual destructor. More...
 
size_t getNumDofPerNode () const
 Gets the number of degree of freedom per node. More...
 
size_t getNumNodes () const
 Gets the number of nodes connected by this element. More...
 
size_t getNodeId (size_t elementNodeId) const
 Gets the elementNodeId-th node id. More...
 
const std::vector< size_t > & getNodeIds () const
 Gets the node ids for this element. More...
 
void setYoungModulus (double E)
 Sets the Young modulus (in N.m-2) More...
 
double getYoungModulus () const
 Gets the Young modulus (in N.m-2) More...
 
void setPoissonRatio (double nu)
 Sets the Poisson ratio (unitless) More...
 
double getPoissonRatio () const
 Gets the Poisson ratio (unitless) More...
 
void setMassDensity (double rho)
 Sets the mass density (in Kg.m-3) More...
 
double getMassDensity () const
 Gets the mass density (in Kg.m-3) More...
 
double getMass (const SurgSim::Math::OdeState &state) const
 Gets the element mass based on the input state (in Kg) More...
 
virtual void addForce (SurgSim::Math::Vector *F, double scale=1.0) const
 Adds the element force (computed for a given state) to a complete system force vector F (assembly) More...
 
virtual void addMass (SurgSim::Math::SparseMatrix *M, double scale=1.0) const
 Adds the element mass matrix M (computed for a given state) to a complete system mass matrix M (assembly) More...
 
virtual void addDamping (SurgSim::Math::SparseMatrix *D, double scale=1.0) const
 Adds the element damping matrix D (= -df/dv) (comuted for a given state) to a complete system damping matrix D (assembly) More...
 
virtual void addStiffness (SurgSim::Math::SparseMatrix *K, double scale=1.0) const
 Adds the element stiffness matrix K (= -df/dx) (computed for a given state) to a complete system stiffness matrix K (assembly) More...
 
virtual void addFMDK (SurgSim::Math::Vector *F, SurgSim::Math::SparseMatrix *M, SurgSim::Math::SparseMatrix *D, SurgSim::Math::SparseMatrix *K) const
 Adds the element force vector, mass, stiffness and damping matrices (computed for a given state) into a complete system data structure F, M, D, K (assembly) More...
 
virtual void addMatVec (double alphaM, double alphaD, double alphaK, const SurgSim::Math::Vector &x, SurgSim::Math::Vector *F) const
 Adds the element matrix-vector contribution F += (alphaM.M + alphaD.D + alphaK.K).x (computed for a given state) into a complete system data structure F (assembly) More...
 
bool isValidCoordinate (const SurgSim::Math::Vector &naturalCoordinate) const
 Determines whether a given natural coordinate is valid. More...
 
template<typename DerivedSub , typename T , int Opt, typename Index >
void assembleMatrixBlocks (const DerivedSub &subMatrix, const std::vector< size_t > blockIds, size_t blockSize, Eigen::SparseMatrix< T, Opt, Index > *matrix, bool initialize=true) const
 Helper method to add a sub-matrix made of squared-blocks into a matrix, for the sake of clarity. More...
 
void updateFMDK (const Math::OdeState &state, int options)
 Update the FemElement based on the given state. More...
 

Protected Member Functions

void initializeMembers ()
 Initializes variables needed before Initialize() is called. More...
 
SurgSim::Math::Matrix33d computeRotation (const SurgSim::Math::OdeState &state)
 Computes the triangle element's rotation given a state. More...
 
virtual void computeLocalStiffness (const SurgSim::Math::OdeState &state, Eigen::Matrix< double, 18, 18 > *localStiffnessMatrix)
 Computes the triangle's local stiffness matrix. More...
 
void computeStiffness (const SurgSim::Math::OdeState &state, SurgSim::Math::Matrix *stiffnessMatrix)
 Computes the triangle's stiffness matrix. More...
 
virtual void computeLocalMass (const SurgSim::Math::OdeState &state, Eigen::Matrix< double, 18, 18 > *localMassMatrix)
 Computes the triangle's local mass matrix. More...
 
void computeMass (const SurgSim::Math::OdeState &state, SurgSim::Math::Matrix *massMatrix)
 Computes the triangle's mass matrix. More...
 
void doUpdateFMDK (const Math::OdeState &state, int options) override
 Update the FemElement based on the given state. More...
 
void computeShapeFunctionsParameters (const SurgSim::Math::OdeState &restState)
 Compute the various shape functions (membrane and plate deformations) parameters. More...
 
std::array< double, 9 > batozDhxDxi (double xi, double eta) const
 Batoz derivative dHx/dxi. More...
 
std::array< double, 9 > batozDhxDeta (double xi, double eta) const
 Batoz derivative dHx/deta. More...
 
std::array< double, 9 > batozDhyDxi (double xi, double eta) const
 Batoz derivative dHy/dxi. More...
 
std::array< double, 9 > batozDhyDeta (double xi, double eta) const
 Batoz derivative dHy/deta. More...
 
Matrix39Type batozStrainDisplacement (double xi, double eta) const
 Batoz strain displacement matrix evaluated at a given point. More...
 
- Protected Member Functions inherited from SurgSim::Physics::FemElement
void setNumDofPerNode (size_t numDofPerNode)
 Sets the number of degrees of freedom per node. More...
 
void initializeFMDK ()
 Initialize f, M, D, K variables. More...
 
virtual void doInitializeFMDK ()
 Function to be overridden by the derived classes to initialize the f, M, D, K variables. More...
 

Protected Attributes

Eigen::Matrix< double, 18, 1 > m_x0
 The element's rest state. More...
 
SurgSim::Math::Matrix33d m_initialRotation
 Initial rotation matrix for the element. More...
 
Eigen::Matrix< double, 18, 18 > m_MLocal
 Stiffness matrix (in local coordinate frame) More...
 
Eigen::Matrix< double, 18, 18 > m_KLocal
 Stiffness matrix (in local coordinate frame) More...
 
double m_restArea
 The triangle rest area. More...
 
double m_thickness
 Thickness of the element. More...
 
SurgSim::Math::Matrix33d m_membraneShapeFunctionsParameters
 Membrane (in-plane) deformation. More...
 
SurgSim::Math::Vector3d m_xij
 Batoz variable \(x_{ij} = x_i - x_j\). More...
 
SurgSim::Math::Vector3d m_yij
 Batoz variable \(y_{ij} = y_i - y_j\). More...
 
SurgSim::Math::Vector3d m_lij_sqr
 Batoz variable \(l_{ij}^2 = x_{ij}^2 + y_{ij}^2\). More...
 
SurgSim::Math::Vector3d m_ak
 Batoz variable \(a_{k} = -x_{ij}/l_i^2\). More...
 
SurgSim::Math::Vector3d m_bk
 Batoz variable \(b_{k} = 3/4 x_{ij} y_{ij}/l_{ij}2\). More...
 
SurgSim::Math::Vector3d m_ck
 Batoz variable \(c_{k} = (1/4 x_{ij}^2 - 1/2 y_{ij}^2)/l_{ij}^2\). More...
 
SurgSim::Math::Vector3d m_dk
 Batoz variable \(d_{k} = -y_{ij}/l_{ij}^2\). More...
 
SurgSim::Math::Vector3d m_ek
 Batoz variable \(e_{k} = (1/4 y_{ij}^2 - 1/2 x_{ij}^2)/l_{ij}^2\). More...
 
SurgSim::Math::Vector3d m_Pk
 Batoz variable \(P_{k} = -6x_{ij}/l_{ij}^2 = 6.\textrm{m_a}_k\). More...
 
SurgSim::Math::Vector3d m_qk
 Batoz variable \(q_{k} = 3x_{ij}y_{ij}/l_{ij}^2 = 4.\textrm{m_b}_k\). More...
 
SurgSim::Math::Vector3d m_tk
 Batoz variable \(t_{k} = -6y_{ij}/l_{ij}^2 = 6.\textrm{m_d}_k\). More...
 
SurgSim::Math::Vector3d m_rk
 Batoz variable \(r_{k} = 3y_{ij}^2/l_{ij}^2\). More...
 
SurgSim::Math::Matrix m_integral_dT_d
 Plate mass matrix: integral terms related to the dof \((z)\). More...
 
SurgSim::Math::Matrix m_integralHyiHyj
 Plate mass matrix: integral terms related to the dof \((\theta_x)\). More...
 
SurgSim::Math::Matrix m_integralHxiHxj
 Plate mass matrix: integral terms related to the dof \((\theta_y)\). More...
 
- Protected Attributes inherited from SurgSim::Physics::FemElement
size_t m_numDofPerNode
 Number of degree of freedom per node for this element. More...
 
std::vector< size_t > m_nodeIds
 Node ids connected by this element. More...
 
double m_rho
 Mass density (in Kg.m-3) More...
 
double m_E
 Young modulus (in N.m-2) More...
 
double m_nu
 Poisson ratio (unitless) More...
 
SurgSim::Math::Vector m_f
 The force vector. More...
 
SurgSim::Math::Matrix m_M
 The mass matrix. More...
 
SurgSim::Math::Matrix m_D
 The damping matrix. More...
 
bool m_useDamping
 Flag to specify of the damping is used. More...
 
SurgSim::Math::Matrix m_K
 The stiffness matrix. More...
 

Private Types

typedef Eigen::Matrix< double, 3, 3 > Matrix33Type
 
typedef Eigen::Matrix< double, 3, 6 > Matrix36Type
 
typedef Eigen::Matrix< double, 6, 6 > Matrix66Type
 
typedef Eigen::Matrix< double, 3, 9 > Matrix39Type
 
typedef Eigen::Matrix< double, 9, 9 > Matrix99Type
 

Private Member Functions

void computeLocalMembraneMass (const SurgSim::Math::OdeState &state, Eigen::Matrix< double, 18, 18 > *localMassMatrix)
 Computes the triangle's local membrane part mass matrix. More...
 
void computeLocalPlateMass (const SurgSim::Math::OdeState &state, Eigen::Matrix< double, 18, 18 > *localMassMatrix)
 Computes the triangle's local plate part mass matrix. More...
 
void computeIntegral_dTd ()
 Computes the integral terms d^T.d over the parametrized triangle area. More...
 
void computeIntegral_HxHxT ()
 Computes the integral terms Hy.Hy^T over the parametrized triangle area. More...
 
void computeIntegral_HyHyT ()
 Computes the integral terms Hx.Hx^T over the parametrized triangle area. More...
 

Additional Inherited Members

- Public Types inherited from SurgSim::Physics::FemElement
typedef SurgSim::Framework::ObjectFactory1< FemElement, std::shared_ptr< FemElementStructs::FemElementParameter > > FactoryType
 
- Static Public Member Functions inherited from SurgSim::Physics::FemElement
static FactoryTypegetFactory ()
 

Detailed Description

2D FemElement based on a triangle with a constant thickness

The triangle is modeled as a shell (6DOF) which is decomposed into a membrane (in-plane 2DOF \((x, y)\)) and a plate (bending/twisting 3DOF \((z, \theta_x, \theta_y)\)). The thin-plate assumption does not consider the drilling dof \(\theta_z\). The system includes the drilling DOF for completeness but does not assign any mass or stiffness to it.

The membrane (in-plane) equations (mass and stiffness) are following "Theory of Matrix Structural Analysis" from J.S. Przemieniecki.

The thin-plate (bending) equations (mass and stiffness) are following "A Study Of Three-Node Triangular Plate Bending Elements", Jean-Louis Batoz Numerical Methods in Engineering, vol 15, 1771-1812 (1980).
The plate mass matrix is not detailed in the above paper, but the analytical equations have been derived from it. Moreover, to account for contribution of the displacement along z to the plate mass matrix, we used a cubic expression of this displacement given in: "Shell elements: modelizations DKT, DST, DKTG and Q4g", Code_Aster, 2013, Thomas De Soza.

Note
The element is considered to have a constant thickness.
The element uses linear elasticity (not visco-elasticity), so it does not have any damping.

Member Typedef Documentation

◆ Matrix33Type

typedef Eigen::Matrix<double, 3, 3> SurgSim::Physics::Fem2DElementTriangle::Matrix33Type
private

◆ Matrix36Type

typedef Eigen::Matrix<double, 3, 6> SurgSim::Physics::Fem2DElementTriangle::Matrix36Type
private

◆ Matrix39Type

typedef Eigen::Matrix<double, 3, 9> SurgSim::Physics::Fem2DElementTriangle::Matrix39Type
private

◆ Matrix66Type

typedef Eigen::Matrix<double, 6, 6> SurgSim::Physics::Fem2DElementTriangle::Matrix66Type
private

◆ Matrix99Type

typedef Eigen::Matrix<double, 9, 9> SurgSim::Physics::Fem2DElementTriangle::Matrix99Type
private

Constructor & Destructor Documentation

◆ Fem2DElementTriangle() [1/3]

SurgSim::Physics::Fem2DElementTriangle::Fem2DElementTriangle ( )

Constructor.

◆ Fem2DElementTriangle() [2/3]

SurgSim::Physics::Fem2DElementTriangle::Fem2DElementTriangle ( std::array< size_t, 3 >  nodeIds)
explicit

Constructor.

Parameters
nodeIdsAn array of 3 node ids defining this triangle element with respect to a DeformableRepresentationState which is passed to the initialize method.

◆ Fem2DElementTriangle() [3/3]

SurgSim::Physics::Fem2DElementTriangle::Fem2DElementTriangle ( std::shared_ptr< FemElementStructs::FemElementParameter elementData)
explicit

Constructor for FemElement object factory.

Parameters
elementDataA FemElement struct defining this triangle element with respect to a DeformableRepresentationState which is passed to the initialize method.
Exceptions
SurgSim::Framework::AssertionFailureif nodeIds has a size different than 3

Member Function Documentation

◆ batozDhxDeta()

std::array< double, 9 > SurgSim::Physics::Fem2DElementTriangle::batozDhxDeta ( double  xi,
double  eta 
) const
protected

Batoz derivative dHx/deta.

Parameters
xi,etaThe parametric coordinate (in [0 1] and xi+eta<1.0)
Returns
The vector dHx/deta evaluated at (xi, eta)

◆ batozDhxDxi()

std::array< double, 9 > SurgSim::Physics::Fem2DElementTriangle::batozDhxDxi ( double  xi,
double  eta 
) const
protected

Batoz derivative dHx/dxi.

Parameters
xi,etaThe parametric coordinate (in [0 1] and xi+eta<1.0)
Returns
The vector dHx/dxi evaluated at (xi, eta)

◆ batozDhyDeta()

std::array< double, 9 > SurgSim::Physics::Fem2DElementTriangle::batozDhyDeta ( double  xi,
double  eta 
) const
protected

Batoz derivative dHy/deta.

Parameters
xi,etaThe parametric coordinate (in [0 1] and xi+eta<1.0)
Returns
The vector dHy/deta evaluated at (xi, eta)

◆ batozDhyDxi()

std::array< double, 9 > SurgSim::Physics::Fem2DElementTriangle::batozDhyDxi ( double  xi,
double  eta 
) const
protected

Batoz derivative dHy/dxi.

Parameters
xi,etaThe parametric coordinate (in [0 1] and xi+eta<1.0)
Returns
The vector dHy/dxi evaluated at (xi, eta)

◆ batozStrainDisplacement()

Fem2DElementTriangle::Matrix39Type SurgSim::Physics::Fem2DElementTriangle::batozStrainDisplacement ( double  xi,
double  eta 
) const
protected

Batoz strain displacement matrix evaluated at a given point.

Parameters
xi,etaThe parametric coordinate (in [0 1] and xi+eta<1.0)
Returns
The 3x9 strain displacement matrix evaluated at (xi, eta)

◆ computeCartesianCoordinate()

SurgSim::Math::Vector SurgSim::Physics::Fem2DElementTriangle::computeCartesianCoordinate ( const SurgSim::Math::OdeState state,
const SurgSim::Math::Vector naturalCoordinate 
) const
overridevirtual

Computes a given natural coordinate in cartesian coordinates.

Parameters
stateThe state at which to transform coordinates
naturalCoordinateThe coordinates to transform
Returns
The resultant cartesian coordinates

Implements SurgSim::Physics::FemElement.

◆ computeIntegral_dTd()

void SurgSim::Physics::Fem2DElementTriangle::computeIntegral_dTd ( )
private

Computes the integral terms d^T.d over the parametrized triangle area.

This integral is required in the plate mass matrix computation. The displacement along z is w(x, y) = [d1 d2 d3 d4 d5 d6 d7 d8 d9].U = d.U with di cubic shape functions and U nodal plate displacements.

◆ computeIntegral_HxHxT()

void SurgSim::Physics::Fem2DElementTriangle::computeIntegral_HxHxT ( )
private

Computes the integral terms Hy.Hy^T over the parametrized triangle area.

This integral is required in the plate mass matrix computation. The displacement along thetay is Thetay(x, y) = -dw/dx = betax = Hx^T.U with Hxi quadratic shape functions and U nodal plate displacements.

◆ computeIntegral_HyHyT()

void SurgSim::Physics::Fem2DElementTriangle::computeIntegral_HyHyT ( )
private

Computes the integral terms Hx.Hx^T over the parametrized triangle area.

This integral is required in the plate mass matrix computation. The displacement along thetax is Thetax(x, y) = dw/dy = -betay = -Hy^T.U with Hyi quadratic shape functions and U nodal plate displacements.

◆ computeLocalMass()

void SurgSim::Physics::Fem2DElementTriangle::computeLocalMass ( const SurgSim::Math::OdeState state,
Eigen::Matrix< double, 18, 18 > *  localMassMatrix 
)
protectedvirtual

Computes the triangle's local mass matrix.

Parameters
stateThe state to compute the local mass matrix from
[out]localMassMatrixThe local mass matrix to store the result into

◆ computeLocalMembraneMass()

void SurgSim::Physics::Fem2DElementTriangle::computeLocalMembraneMass ( const SurgSim::Math::OdeState state,
Eigen::Matrix< double, 18, 18 > *  localMassMatrix 
)
private

Computes the triangle's local membrane part mass matrix.

Parameters
stateThe state to compute the local membrane part mass matrix from
[out]localMassMatrixThe local mass matrix to store the result into


The mass matrix is derived from the kinetic energy, which depends on the displacements \(\mathbf{u}\).
In the case of triangle membrane, the displacements are a subset of the full shell displacements:

\[ \displaystyle \mathbf{u}(x, y) = \left( \begin{array}{c} u_x(x, y) \\ u_y(x, y) \end{array} \right) = \left( \begin{array}{c} \sum_{i=0}^2 N_i(x, y).u_x^i \\ \sum_{i=0}^2 N_i(x, y).u_y^i \end{array} \right) = \left( \begin{array}{cccccc} N_0(x, y) & 0 & N_1(x, y) & 0 & N_2(x, y) & 0 \\ 0 & N_0(x, y) & 0 & N_1(x, y) & 0 & N_2(x, y) \end{array} \right) \left( \begin{array}{c} u_x^0 \\ u_y^0 \\ u_x^1 \\ u_y^1 \\ u_x^2 \\ u_y^2 \end{array} \right) = a.\mathbf{U} \]

with \(\mathbf{U}\) the triangle membrane nodal displacements vector (subset of the full shell nodal displacements) and \(N_i(x, y) = a_i + b_i.x + c_i.y\) the triangle linear shape functions.

Considering a constant mass density \(\rho\) and a constant thickness \(t\), the kinetic energy of the triangle membrane is \(T = \int_V \frac{1}{2} \rho(x, y, z) \mathbf{\dot{u}}^T.\mathbf{\dot{u}} \, dV = \frac{1}{2}.\rho.t.\int_A \left(\dot{\mathbf{U}}^T.a^T\right).\left(a.\dot{\mathbf{U}}\right) \, dA\).

The triangle membrane mass matrix can be derived from the kinetic energy using the Lagrange equation of motion:

\[ \displaystyle \frac{d \frac{\partial T}{\partial \dot{\mathbf{U}}}}{d t} = \frac{1}{2}.\rho.t.\frac{d \frac{\partial \displaystyle \int_A \left(\dot{\mathbf{U}}^T.a^T\right).\left(a.\dot{\mathbf{U}}\right) \, dA}{\partial \dot{\mathbf{U}}}}{d t} = \frac{1}{2}.\rho.t.\frac{d \displaystyle \int_A 2.a^T.a.\dot{\mathbf{U}} \, dA}{d t} = \left( \rho.t.\int_A a^T.a \, dA \right).\ddot{\mathbf{U}} \]

Therefore, the mass matrix is:

\[ M = \rho.t.\int_A a^T.a \, dA \]

To integrate over the triangle area, we operate the following variables change:

\[ \left\{ \begin{array}{ccc} x(\xi, \eta) & = & x_0 + \xi.(x_1 - x_0) + \eta.(x_2 - x_0) \\ y(\xi, \eta) & = & y_0 + \xi.(y_1 - y_0) + \eta.(y_2 - y_0) \end{array} \right. \]

Therefore, we have \( J = \left( \begin{array}{cc} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{array} \right) = \left( \begin{array}{cc} (x_1 - x_0) & (x_2 - x_0) \\ (y_1 - y_0) & (y_2 - y_0) \\ \end{array} \right) \) and \(|J| = 2A\). Therefore the mass matrix becomes

\[ M = \rho.t.\int_0^1 \int_0^{1-\eta} a^T.a |J| \, d\xi \, d\eta = 2.A.\rho.t.\int_0^1 \int_0^{1-\eta} \left( \begin{array}{cccccc} N0.N0 & 0 & N0.N1 & 0 & N0.N2 & 0 \\ 0 & N0.N0 & 0 & N0.N1 & 0 & N0.N2 \\ N1.N0 & 0 & N1.N1 & 0 & N1.N2 & 0 \\ 0 & N1.N0 & 0 & N1.N1 & 0 & N1.N2 \\ N2.N0 & 0 & N2.N1 & 0 & N2.N2 & 0 \\ 0 & N2.N0 & 0 & N2.N1 & 0 & N2.N2 \end{array} \right) \, d\xi \, d\eta \]

Considering the triangle in its local space (i.e. \(x_0=y_0=y_1=0\), \(x_1>0\), \(y_2>0\)), we get:

\[ M = 2.A.\rho.t. \left( \begin{array}{cccccc} \frac{1.0}{12.0} & 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{24.0} & 0 \\ 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{24.0} \\ \frac{1.0}{24.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{24.0} & 0 \\ 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{24.0} \\ \frac{1.0}{24.0} & 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{12.0} & 0 \\ 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{24.0} & 0 & \frac{1.0}{12.0} \end{array} \right) = \rho.A.t. \left( \begin{array}{cccccc} \frac{1.0}{6.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{12.0} & 0 \\ 0 & \frac{1.0}{6.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{12.0} \\ \frac{1.0}{12.0} & 0 & \frac{1.0}{6.0} & 0 & \frac{1.0}{12.0} & 0 \\ 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{6.0} & 0 & \frac{1.0}{12.0} \\ \frac{1.0}{12.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{6.0} & 0 \\ 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{12.0} & 0 & \frac{1.0}{6.0} \end{array} \right) \]

which needs to be assembled properly in the full shell element mass matrix.

◆ computeLocalPlateMass()

void SurgSim::Physics::Fem2DElementTriangle::computeLocalPlateMass ( const SurgSim::Math::OdeState state,
Eigen::Matrix< double, 18, 18 > *  localMassMatrix 
)
private

Computes the triangle's local plate part mass matrix.

Parameters
stateThe state to compute the local plate part mass matrix from
[out]localMassMatrixThe local mass matrix to store the result into
Note
The plate mass matrix is composed of 3 matrices associated respectively to
displacements in direction (z, thetax, thetay).


The mass matrix is derived from the kinetic energy, which depends on the displacements \(\mathbf{u}\).
In the case of triangle plate, the displacements are a subset of the full shell displacements:

\[ \displaystyle \mathbf{u}(x, y, z) = \left( \begin{array}{c} u(x, y, z) = z.\beta_x(x, y) \\ v(x, y, z) = z.\beta_y(x, y) \\ w(x, y) \end{array} \right) \]

Note that we use the plate formulation DKT (Discrete Kirchhoff Theory) developped by Batoz in 1980. In this formulation, the displacement \(w\) is never formaly expressed, we only know that it is a cubic polynomial function of \((x, y)\). To complete the plate mass matrix calculation, we tried to use the expression of \(w\) from Przemieniecki book "Theory of Matrix Structural Analysis" Chapter 11.10, equation 11.57 (i.e. \(w(x, y) = \mathbf{d}.C^{-1}.\mathbf{U}\)). Unfortunately, it can occur that the matrix \(C\) becomes singular (example of a right isoceles triangle, i.e. \(x_0=y_0=y_1=x_2=0\) and \(x_1=y_2=1\)). To overcome this issue, we prefered the Code_Aster approach developped in "Shell elements: modelizations DKT, DST, DKTG and Q4g", Code_Aster, 2013, Thomas De Soza, where a cubic well defined expression of \(w(x, y)\) is given as: \(w(x, y) = \mathbf{d}.\mathbf{U}\) where \(d_i\) are cubic shape functions.

Therefore, it leads to the following plate displacements:

\[ \mathbf{u}(x, y, z) = \left( \begin{array}{c} u(x, y, z) = z.\beta_x(x, y) = z.\mathbf{H_x}^T(\xi, \eta).\mathbf{U} \\ v(x, y, z) = z.\beta_y(x, y) = z.\mathbf{H_y}^T(\xi, \eta).\mathbf{U} \\ w(x, y) = \mathbf{d}.\mathbf{U} \end{array} \right) = \left( \begin{array}{c} z.\mathbf{H_x}^T(\xi, \eta) \\ z.\mathbf{H_y}^T(\xi, \eta) \\ \mathbf{d} \end{array} \right) . \mathbf{U} = a.\mathbf{U} \]

with \(\mathbf{U}\) the triangle plate nodal displacements vector (subset of the full shell nodal displacements) and \(\mathbf{H_x}\) and \(\mathbf{H_y}\) are the Batoz shape functions.

Considering a constant mass density \(\rho\) and a constant thickness \(t\), the kinetic energy of the triangle plate is \(T = \int_V \frac{1}{2} \rho(x, y, z) \mathbf{\dot{u}}^T.\mathbf{\dot{u}} \, dV = \frac{1}{2}.\rho.\int_{-t/2}^{+t/2} \int_A \left(\dot{\mathbf{U}}^T.a^T\right).\left(a.\dot{\mathbf{U}}\right) \, dA \, dz \).

The triangle plate mass matrix can be derived from the kinetic energy using the Lagrange equation of motion:

\[ \displaystyle \frac{d \frac{\partial T}{\partial \dot{\mathbf{U}}}}{d t} = \frac{1}{2}.\rho.\frac{d \frac{\partial \displaystyle \int_{-t/2}^{+t/2} \int_A \left(\dot{\mathbf{U}}^T.a^T\right).\left(a.\dot{\mathbf{U}}\right) \, dA \, dz}{\partial \dot{\mathbf{U}}}}{d t} = \frac{1}{2}.\rho.\frac{d \displaystyle \int_{-t/2}^{+t/2} \int_A 2.a^T.a.\dot{\mathbf{U}} \, dA \, dz}{d t} \\ = \left( \rho.\int_{-t/2}^{+t/2} \int_A a^T.a \, dA \, dz \right).\ddot{\mathbf{U}} = \left( \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_x}).(z.\mathbf{H_x}^T) + (z.\mathbf{H_y}).(z.\mathbf{H_y}^T) + \mathbf{d}^T).\mathbf{d} \, dA \, dz \right).\ddot{\mathbf{U}} \]

Therefore, the mass matrix is:

\[ M = \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_x}).(z.\mathbf{H_x}^T) \, dA \, dz + \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_y}).(z.\mathbf{H_y}^T) \, dA \, dz + \rho.\int_{-t/2}^{+t/2} \int_A \mathbf{d}^T.\mathbf{d} \, dA \, dz \]

It clearly appears that the triangle plate mass matrix is composed of 3 matrices, each providing the contribution of a displacement (i.e. \((u, v, w)\)) associated to a given degree of freedom (i.e. respectively \((\theta_y, \theta_x, z)\)). Let's note these 3 matrices \(M_u\), \(M_v\) and \(M_w\).

To integrate over the triangle area, we operate the following variables change:

\[ \left\{ \begin{array}{ccc} x(\xi, \eta) & = & x_0 + \xi.(x_1 - x_0) + \eta.(x_2 - x_0) \\ y(\xi, \eta) & = & y_0 + \xi.(y_1 - y_0) + \eta.(y_2 - y_0) \end{array} \right. \]

Therefore, we have \( J = \left( \begin{array}{cc} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{array} \right) = \left( \begin{array}{cc} (x_1 - x_0) & (x_2 - x_0) \\ (y_1 - y_0) & (y_2 - y_0) \\ \end{array} \right) \) and \(|J| = 2A\). Therefore the mass matrices become (noting that the shape functions \(d_i\) do not depend on \(z\)):

\[ \begin{array}{ccc} M_u & = & \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_x}).(z.\mathbf{H_x}^T) \, dA \, dz = \rho.\int_{-t/2}^{+t/2} z^2 \, dz.\int_A \mathbf{H_x}.\mathbf{H_x}^T \, dA \\ & = & \rho.\frac{t^3}{12}.\int_0^1 \int_0^{1-\eta} \mathbf{H_x}.\mathbf{H_x}^T |J| \, d\xi d\eta = 2.A.\rho.\frac{t^3}{12}.\int_0^1 \int_0^{1-\eta} \mathbf{H_x}.\mathbf{H_x}^T \, d\xi d\eta \\ M_v & = & \rho.\int_{-t/2}^{+t/2} \int_A (z.\mathbf{H_y}).(z.\mathbf{H_y}^T) \, dA \, dz = \rho.\int_{-t/2}^{+t/2} z^2 \, dz.\int_A \mathbf{H_y}.\mathbf{H_y}^T \, dA \\ & = & \rho.\frac{t^3}{12}.\int_0^1 \int_0^{1-\eta} \mathbf{H_y}.\mathbf{H_y}^T |J| \, d\xi d\eta = 2.A.\rho.\frac{t^3}{12}.\int_0^1 \int_0^{1-\eta} \mathbf{H_y}.\mathbf{H_y}^T \, d\xi d\eta \\ M_w & = & \rho.\int_{-t/2}^{+t/2} \int_A \mathbf{d}^T.\mathbf{d} \, dA \, dz = 2.A.\rho.t.\int_0^1 \int_0^{1-\eta} \mathbf{d}^T.\mathbf{d} \, d\xi d\eta \end{array} \]

To evaluate the 3 integral terms

\[ \begin{array}{ccccc} \int_0^1 \int_0^{1-\eta} \mathbf{H_x}.\mathbf{H_x}^T \, d\xi d\eta & , & \int_0^1 \int_0^{1-\eta} \mathbf{H_y}.\mathbf{H_y}^T \, d\xi d\eta & , & \int_0^1 \int_0^{1-\eta} \mathbf{d}^T.\mathbf{d} \, d\xi d\eta \end{array} \]

we used a formal mathematical package and tested the corresponding code against a numerically exact evaluation using the appropriate Gauss quadrature on the triangle.

The resulting triangle plate mass matrix needs to be assembled properly in the full shell element mass matrix.

◆ computeLocalStiffness()

void SurgSim::Physics::Fem2DElementTriangle::computeLocalStiffness ( const SurgSim::Math::OdeState state,
Eigen::Matrix< double, 18, 18 > *  localStiffnessMatrix 
)
protectedvirtual

Computes the triangle's local stiffness matrix.

Parameters
stateThe state to compute the local stiffness matrix from
[out]localStiffnessMatrixThe local stiffness matrix to store the result into

◆ computeMass()

void SurgSim::Physics::Fem2DElementTriangle::computeMass ( const SurgSim::Math::OdeState state,
SurgSim::Math::Matrix massMatrix 
)
protected

Computes the triangle's mass matrix.

Parameters
stateThe state to compute the mass matrix from
[out]massMatrixThe mass matrix to store the result into

◆ computeNaturalCoordinate()

SurgSim::Math::Vector SurgSim::Physics::Fem2DElementTriangle::computeNaturalCoordinate ( const SurgSim::Math::OdeState state,
const SurgSim::Math::Vector cartesianCoordinate 
) const
overridevirtual

Computes a natural coordinate given a global coordinate.

Parameters
stateThe state at which to transform coordinates
cartesianCoordinateThe coordinates to transform
Returns
The resultant natural coordinates

Implements SurgSim::Physics::FemElement.

◆ computeRotation()

SurgSim::Math::Matrix33d SurgSim::Physics::Fem2DElementTriangle::computeRotation ( const SurgSim::Math::OdeState state)
protected

Computes the triangle element's rotation given a state.

Parameters
stateThe state to compute the rotation from
Returns
The rotation matrix of the element in the given state

◆ computeShapeFunctionsParameters()

void SurgSim::Physics::Fem2DElementTriangle::computeShapeFunctionsParameters ( const SurgSim::Math::OdeState restState)
protected

Compute the various shape functions (membrane and plate deformations) parameters.

Parameters
restStatethe rest state to compute the shape functions parameters from

◆ computeStiffness()

void SurgSim::Physics::Fem2DElementTriangle::computeStiffness ( const SurgSim::Math::OdeState state,
SurgSim::Math::Matrix stiffnessMatrix 
)
protected

Computes the triangle's stiffness matrix.

Parameters
stateThe state to compute the stiffness matrix from
[out]stiffnessMatrixThe stiffness matrix to store the result into

◆ doUpdateFMDK()

void SurgSim::Physics::Fem2DElementTriangle::doUpdateFMDK ( const Math::OdeState state,
int  options 
)
overrideprotectedvirtual

Update the FemElement based on the given state.

Parameters
state\((x, v)\) the current position and velocity to evaluate the various terms with
optionsFlag to specify which of the F, M, D, K needs to be updated

Implements SurgSim::Physics::FemElement.

◆ getThickness()

double SurgSim::Physics::Fem2DElementTriangle::getThickness ( ) const

Gets the triangle's thickness.

Returns
The thickness of the triangle

◆ getVolume()

double SurgSim::Physics::Fem2DElementTriangle::getVolume ( const SurgSim::Math::OdeState state) const
overridevirtual

Gets the element volume based on the input state (in m-3)

Parameters
stateThe state to compute the volume with
Returns
The volume of this element (in m-3)

Implements SurgSim::Physics::FemElement.

◆ initialize()

void SurgSim::Physics::Fem2DElementTriangle::initialize ( const SurgSim::Math::OdeState state)
overridevirtual

Initialize the FemElement once everything has been set.

Parameters
stateThe state to initialize the FemElement with

Reimplemented from SurgSim::Physics::FemElement.

◆ initializeMembers()

void SurgSim::Physics::Fem2DElementTriangle::initializeMembers ( )
protected

Initializes variables needed before Initialize() is called.

◆ setThickness()

void SurgSim::Physics::Fem2DElementTriangle::setThickness ( double  thickness)

Sets the triangle's thickness.

Parameters
thicknessThe thickness of the triangle

Member Data Documentation

◆ m_ak

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_ak
protected

Batoz variable \(a_{k} = -x_{ij}/l_i^2\).

◆ m_bk

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_bk
protected

Batoz variable \(b_{k} = 3/4 x_{ij} y_{ij}/l_{ij}2\).

◆ m_ck

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_ck
protected

Batoz variable \(c_{k} = (1/4 x_{ij}^2 - 1/2 y_{ij}^2)/l_{ij}^2\).

◆ m_dk

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_dk
protected

Batoz variable \(d_{k} = -y_{ij}/l_{ij}^2\).

◆ m_ek

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_ek
protected

Batoz variable \(e_{k} = (1/4 y_{ij}^2 - 1/2 x_{ij}^2)/l_{ij}^2\).

◆ m_initialRotation

SurgSim::Math::Matrix33d SurgSim::Physics::Fem2DElementTriangle::m_initialRotation
protected

Initial rotation matrix for the element.

◆ m_integral_dT_d

SurgSim::Math::Matrix SurgSim::Physics::Fem2DElementTriangle::m_integral_dT_d
protected

Plate mass matrix: integral terms related to the dof \((z)\).

◆ m_integralHxiHxj

SurgSim::Math::Matrix SurgSim::Physics::Fem2DElementTriangle::m_integralHxiHxj
protected

Plate mass matrix: integral terms related to the dof \((\theta_y)\).

◆ m_integralHyiHyj

SurgSim::Math::Matrix SurgSim::Physics::Fem2DElementTriangle::m_integralHyiHyj
protected

Plate mass matrix: integral terms related to the dof \((\theta_x)\).

◆ m_KLocal

Eigen::Matrix<double, 18, 18> SurgSim::Physics::Fem2DElementTriangle::m_KLocal
protected

Stiffness matrix (in local coordinate frame)

◆ m_lij_sqr

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_lij_sqr
protected

Batoz variable \(l_{ij}^2 = x_{ij}^2 + y_{ij}^2\).

◆ m_membraneShapeFunctionsParameters

SurgSim::Math::Matrix33d SurgSim::Physics::Fem2DElementTriangle::m_membraneShapeFunctionsParameters
protected

Membrane (in-plane) deformation.

DOF simulated: (x, y). "Theory of Matrix Structural Analysis" from J.S. Przemieniecki. Shape functions are \(f_i(x, y) = a_i + b_i.x + c_i.y\), here we store \((a_i, b_i, c_i)\) on each row.

◆ m_MLocal

Eigen::Matrix<double, 18, 18> SurgSim::Physics::Fem2DElementTriangle::m_MLocal
protected

Stiffness matrix (in local coordinate frame)

◆ m_Pk

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_Pk
protected

Batoz variable \(P_{k} = -6x_{ij}/l_{ij}^2 = 6.\textrm{m_a}_k\).

◆ m_qk

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_qk
protected

Batoz variable \(q_{k} = 3x_{ij}y_{ij}/l_{ij}^2 = 4.\textrm{m_b}_k\).

◆ m_restArea

double SurgSim::Physics::Fem2DElementTriangle::m_restArea
protected

The triangle rest area.

◆ m_rk

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_rk
protected

Batoz variable \(r_{k} = 3y_{ij}^2/l_{ij}^2\).

◆ m_thickness

double SurgSim::Physics::Fem2DElementTriangle::m_thickness
protected

Thickness of the element.

◆ m_tk

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_tk
protected

Batoz variable \(t_{k} = -6y_{ij}/l_{ij}^2 = 6.\textrm{m_d}_k\).

◆ m_x0

Eigen::Matrix<double, 18, 1> SurgSim::Physics::Fem2DElementTriangle::m_x0
protected

The element's rest state.

◆ m_xij

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_xij
protected

Batoz variable \(x_{ij} = x_i - x_j\).

◆ m_yij

SurgSim::Math::Vector3d SurgSim::Physics::Fem2DElementTriangle::m_yij
protected

Batoz variable \(y_{ij} = y_i - y_j\).


The documentation for this class was generated from the following files: