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:

## 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 FactoryType & | getFactory ()

## 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

  nodeIds | An 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

  elementData | A FemElement struct defining this triangle element with respect to a DeformableRepresentationState which is passed to the initialize method.

Exceptions

  SurgSim::Framework::AssertionFailure | if 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,eta | The 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,eta | The 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,eta | The 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,eta | The 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,eta | The 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

override virtual

Computes a given natural coordinate in cartesian coordinates.

Parameters

  state | The state at which to transform coordinates
  naturalCoordinate | The 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
| ) |  |

protected virtual

Computes the triangle's local mass matrix.

Parameters

  | state | The state to compute the local mass matrix from
  [out] | localMassMatrix | The 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

  | state | The state to compute the local membrane part mass matrix from
  [out] | localMassMatrix | The 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

  | state | The state to compute the local plate part mass matrix from
  [out] | localMassMatrix | The 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
| ) |  |

protected virtual

Computes the triangle's local stiffness matrix.

Parameters

  | state | The state to compute the local stiffness matrix from
  [out] | localStiffnessMatrix | The 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

  | state | The state to compute the mass matrix from
  [out] | massMatrix | The 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

override virtual

Computes a natural coordinate given a global coordinate.

Parameters

  state | The state at which to transform coordinates
  cartesianCoordinate | The 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

  state | The 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

  restState | the 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

  | state | The state to compute the stiffness matrix from
  [out] | stiffnessMatrix | The stiffness matrix to store the result into

## doUpdateFMDK()

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

override protected virtual

Update the FemElement based on the given state.

Parameters

  state | \((x, v)\) the current position and velocity to evaluate the various terms with
  options | Flag 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

override virtual

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

Parameters

  state | The 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 | ) |

override virtual

Initialize the FemElement once everything has been set.

Parameters

  state | The 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

  thickness | The 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:

- SurgSim/Physics/Fem2DElementTriangle.h
- SurgSim/Physics/Fem2DElementTriangle.cpp
- SurgSim/Physics/Fem2DElementTriangle_computeLocalMembraneMass.dox
- SurgSim/Physics/Fem2DElementTriangle_computeLocalPlateMass.dox