artisynth.core.femmodels

Class FemUtilities

java.lang.Object

artisynth.core.femmodels.FemUtilities

public class FemUtilities
extends java.lang.Object

Provides some general utilities for FEM computations. Some of these involve the B matrix formed from the shape function gradient dNdx. If this gradient is represented as g, and if the stress tensor is represented with the vector

   sigma = [ sxx syy szz sxy syz sxz ]
 

then the B matrix takes the form

      [ gx  0  0 ]           00 11 22 01 12 02
      [  0 gy  0 ]
      [  0  0 gz ]     T   [ gx  0  0 gy  0 gz ]
  B = [ gy gx  0 ]    B  = [  0 gy  0 gx gz  0 ]
      [  0 gz gy ]         [  0  0 gz  0 gy gx ]
      [ gz  0 gx ]
 

Constructor Summary

Constructors

Constructor and Description
FemUtilities()

Method Summary

Modifier and Type	Method and Description
static void	addDilationalStiffness(Matrix3d K, double kp, Matrix3x1 GT_i, Matrix3x1 GT_j) Adds dilational stiffness to the node-to-node stiffness matrix Kij.
static void	addDilationalStiffness(Matrix3d K, double kp, Vector3d intGi, Vector3d intGj) Adds dilational stiffness to the node-to-node stiffness matrix Kij.
static void	addDilationalStiffness(Matrix3d K, MatrixNd Rinv, MatrixBlock GT_i, MatrixBlock GT_j) Adds dilational stiffness to the node-to-node stiffness matrix Kij.
static void	addGeometricStiffness(Matrix3d K, Vector3d gi, SymmetricMatrix3d sig, Vector3d gj, double dv) Adds the geometric stiffness defined by
static void	addIncompressibilityStiffness(Matrix3d K, double s, Vector3d intGi, Vector3d intGj) Adds stiffness for the incompressibility constraint to the node-to-node stiffness matrix Kij.
static void	addMaterialStiffness(Matrix3d K, Vector3d gi, double E, double nu, Vector3d gj, double dv) Adds a weighted node-to-node stiffness to the matrix Kij via the formula
static void	addMaterialStiffness(Matrix3d K, Vector3d gi, Matrix6d D, SymmetricMatrix3d sig, Vector3d gj, double dv) Adds a weighted node-to-node stiffness to the matrix Kij via the formula
static void	addMaterialStiffness(Matrix3d K, Vector3d gi, Matrix6d D, Vector3d gj, double dv) Adds a weighted node-to-node stiffness to the matrix Kij via the formula
static void	addMembraneMaterialStiffness(Matrix3d K00, Vector3d idN, Vector3d jdN, double dv, Matrix3d invJ, SymmetricMatrix3d matStress, Matrix6d matTangent)
static void	addMembraneStressForce(Vector3d f, SymmetricMatrix3d sig, double dv, double dNdr, double dNds, Matrix3d invJ)
static void	addPressureStiffness(Matrix3d K, Vector3d gi, double p, Vector3d gj, double dv) Adds pressure stiffness to the matrix Kij via the formula
static void	addShellMaterialStiffness(Matrix3d K00, Matrix3d K01, Matrix3d K10, Matrix3d K11, double iN, double jN, Vector3d idN, Vector3d jdN, double dv, double t, Matrix3d invJ, SymmetricMatrix3d matStress, Matrix6d matTangent) Add weighted material stiffness for this (i,j) node neighbor pair (represented by 6x6 stiffness block), relative to a particular integration point of the shell element.
static void	addShellStressForce(Vector3d f, Vector3d df, SymmetricMatrix3d sig, double t, double dv, double N, double dNdr, double dNds, Matrix3d invJ) Adds the material+geometric forces on a node resulting from a given stress at a given shell integration point.
static void	addStressForce(Vector3d f, Vector3d g, SymmetricMatrix3d sig, double dv) Adds the force on a node resulting from a given stress and shape function gradient g.
static void	addToIncompressConstraints(MatrixBlock blk, double[] H, Vector3d GNx, double dv) Adds H^T GNx dv to a matrix block, where H is a row vector of weight values, and dv is a volume differential.
static double	computeJ2(SymmetricMatrix3d M) Compute the second deviatoric invariant J2 https://en.wikipedia.org/wiki/Cauchy_stress_tensor#Stress_deviator_tensor
static double	computeVonMisesStrain(SymmetricMatrix3d strain) Compute the Von Mises strain equivalent according to http://www.continuummechanics.org/vonmisesstress.html which is equivalent to https://dianafea.com/manuals/d944/Analys/node405.html The Von Mises Strain Equivalent is equal to sqrt(4/3 J2), where J2 is the second invariant of the average deviatoric strain for the node.
static double	computeVonMisesStress(SymmetricMatrix3d sig) Compute the Von Mises Stress criterion https://en.wikipedia.org/wiki/Von_Mises_yield_criterion The Von Mises Stress is equal to sqrt(3 J2), where J2 is the second invariant of the average deviatoric strain for the node.
static FemNode3d[][]	triangulate6NodeFace(FaceNodes3d face)
static FemNode3d[][]	triangulate8NodeFace(FaceNodes3d face)
static int[]	triangulateFaceIndices(int[] faces) Helper method to triangulate a set of faces.

Methods inherited from class java.lang.Object

equals, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Constructor Detail

FemUtilities

public FemUtilities()

Method Detail

addStressForce

public static void addStressForce(Vector3d f,
                                  Vector3d g,
                                  SymmetricMatrix3d sig,
                                  double dv)

Adds the force on a node resulting from a given stress and shape function gradient g. This is done by computing f = B^T sig, where B is the matrix formed from the shape function gradient.

addMaterialStiffness

public static void addMaterialStiffness(Matrix3d K,
                                        Vector3d gi,
                                        Matrix6d D,
                                        SymmetricMatrix3d sig,
                                        Vector3d gj,
                                        double dv)

Adds a weighted node-to-node stiffness to the matrix Kij via the formula

  Kij = (Bi^T D Bj + (gi^T sig gj) I)*dv
 

where gi and gj are the gradients for shape functions i and j, Bi and Bj are the B matrices formed from gi and gj, D is the derivative of the constitutive relationship, and dv is the weighting term.

addMaterialStiffness

public static void addMaterialStiffness(Matrix3d K,
                                        Vector3d gi,
                                        Matrix6d D,
                                        Vector3d gj,
                                        double dv)

Adds a weighted node-to-node stiffness to the matrix Kij via the formula

  Kij = Bi^T D Bj dv
 

where gi and gj are the gradients for shape functions i and j, Bi and Bj are the B matrices formed from gi and gj, D is the derivative of the constitutive relationship, and dv is the weighting term.

addGeometricStiffness

public static void addGeometricStiffness(Matrix3d K,
                                         Vector3d gi,
                                         SymmetricMatrix3d sig,
                                         Vector3d gj,
                                         double dv)

Adds the geometric stiffness defined by

 
 Kij += (gi^T sig gj) I*dv
 

addMaterialStiffness

public static void addMaterialStiffness(Matrix3d K,
                                        Vector3d gi,
                                        double E,
                                        double nu,
                                        Vector3d gj,
                                        double dv)

Adds a weighted node-to-node stiffness to the matrix Kij via the formula

  Kij = Bi^T D Bj dv
 

where gi and gj are the gradients for shape functions i and j, Bi and Bj are the B matrices formed from gi and gj, D is the linear stiffness relationship associated with Youngs modulus E and Poissons ratio nu, and dv is the weighting term.

addPressureStiffness

public static void addPressureStiffness(Matrix3d K,
                                        Vector3d gi,
                                        double p,
                                        Vector3d gj,
                                        double dv)

Adds pressure stiffness to the matrix Kij via the formula

  Kij = Bi^T D Bj
 

where gi and gj are the gradients for shape functions i and j, Bi and Bj are the B matrices formed from gi and gj, and D is created from the pressure p according to D_ii = -p and D_ij = p for i != j and i, j < 3.

addDilationalStiffness

public static void addDilationalStiffness(Matrix3d K,
                                          double kp,
                                          Vector3d intGi,
                                          Vector3d intGj)

Adds dilational stiffness to the node-to-node stiffness matrix Kij. This is calculated via the formula

  Kij += kp (intGi intGj^T)
 

where intGi and intGj are the integrals (over the element) of the gradients for shape functions i and j, and kp is (typically) the bulkModulus divided by the restVolume.

addDilationalStiffness

public static void addDilationalStiffness(Matrix3d K,
                                          double kp,
                                          Matrix3x1 GT_i,
                                          Matrix3x1 GT_j)

Adds dilational stiffness to the node-to-node stiffness matrix Kij. This is calculated via the formula

  Kij += kp * (GT_i GT_j^T)
 

where GT_i and GT_j are the incompressibility constraints for nodes i and j and kp is (typically) the bulkModulus divided by the restVolume.

addDilationalStiffness

public static void addDilationalStiffness(Matrix3d K,
                                          MatrixNd Rinv,
                                          MatrixBlock GT_i,
                                          MatrixBlock GT_j)

Adds dilational stiffness to the node-to-node stiffness matrix Kij. This is calculated via the formula

  Kij += GT_i Rinv GT_j^T
 

where GT_i and GT_j are the incompressibility constraints for nodes i and j and Kp is a diagonal matrix of stiffness pressures.

addIncompressibilityStiffness

public static void addIncompressibilityStiffness(Matrix3d K,
                                                 double s,
                                                 Vector3d intGi,
                                                 Vector3d intGj)

Adds stiffness for the incompressibility constraint to the node-to-node stiffness matrix Kij. This is calculated via the formula

  Kij += kp (intGi intGj^T - intGj intGi^T)
 

where intGi and intGj are the integrals (over the element) of the gradients for shape functions i and j, and s is (typically) the determinant of F times the current pressure.

addToIncompressConstraints

public static void addToIncompressConstraints(MatrixBlock blk,
                                              double[] H,
                                              Vector3d GNx,
                                              double dv)

Adds H^T GNx dv to a matrix block, where H is a row vector of weight values, and dv is a volume differential. It is assumed that the matrix block has a row size of 3.

triangulate8NodeFace

public static FemNode3d[][] triangulate8NodeFace(FaceNodes3d face)

triangulate6NodeFace

public static FemNode3d[][] triangulate6NodeFace(FaceNodes3d face)

triangulateFaceIndices

public static int[] triangulateFaceIndices(int[] faces)

Helper method to triangulate a set of faces. Each face is given by a set of nodes, arranged counter-clockwise about what would be its outer facing normal if all the nodes were coplanar. Information for all faces is contained in faces, with each entry consisting of the number of nodes, followed by the indices for each node.

addShellMaterialStiffness

public static void addShellMaterialStiffness(Matrix3d K00,
                                             Matrix3d K01,
                                             Matrix3d K10,
                                             Matrix3d K11,
                                             double iN,
                                             double jN,
                                             Vector3d idN,
                                             Vector3d jdN,
                                             double dv,
                                             double t,
                                             Matrix3d invJ,
                                             SymmetricMatrix3d matStress,
                                             Matrix6d matTangent)

Add weighted material stiffness for this (i,j) node neighbor pair (represented by 6x6 stiffness block), relative to a particular integration point of the shell element. This material stiffness also accounts for geometrical stiffness. FEBio: FEElasticShellDomain::ElementStiffness

Parameters:

K00 - front-front stiffness block to accumulate stiffness.

K01 - front-back stiffness block to accumulate stiffness.

K10 - back-front stiffness block to accumulate stiffness.

K11 - back-back stiffness block to accumulate stiffness.

iN - Shape function of i-th node and integration point.

jN - Shape function of j-th node and integration point.

idN - Derivative of shape function (wrt natural coords) of i-th node and integration point.

jdN - Derivative of shape function (wrt natural coords) of j-th node and integration point.

dv - integrationPt.detJ * integrationPt.weight

t - t-component of (r,s,t) integration point coordinates (i.e. x-component of gauss point)

invJ - Current inverse Jacobian at the integration point

matStress - Material stress of integration point.

matTangent - Material tangent of integration point. Postcond: Stiffness is applied to K

addMembraneMaterialStiffness

public static void addMembraneMaterialStiffness(Matrix3d K00,
                                                Vector3d idN,
                                                Vector3d jdN,
                                                double dv,
                                                Matrix3d invJ,
                                                SymmetricMatrix3d matStress,
                                                Matrix6d matTangent)

addShellStressForce

public static void addShellStressForce(Vector3d f,
                                       Vector3d df,
                                       SymmetricMatrix3d sig,
                                       double t,
                                       double dv,
                                       double N,
                                       double dNdr,
                                       double dNds,
                                       Matrix3d invJ)

Adds the material+geometric forces on a node resulting from a given stress at a given shell integration point. FEBio: FEElasticShellDomain::ElementInternalForce

Parameters:

f - Displacement force vector (x,y,z) of node to append

df - Direction force vector (u,w,v) of node to append

sig - Computed material stress

t - t-component of (r,s,t) integration point coordinates (i.e. x-component of gauss point)

dv - integrationPt.detJ * integrationPt.weight

N - Shape function of node and integration point.

dNdr - x (i.e. r) component of derivative of shape function of node and integration point.

dNds - y (i.e. s) component of derivative of shape function of node and integration point.

invJ - Current inverse Jacobian at the integration point Postcond: Stress forces are applied to f and df

addMembraneStressForce

public static void addMembraneStressForce(Vector3d f,
                                          SymmetricMatrix3d sig,
                                          double dv,
                                          double dNdr,
                                          double dNds,
                                          Matrix3d invJ)

computeJ2

public static double computeJ2(SymmetricMatrix3d M)

Compute the second deviatoric invariant J2 https://en.wikipedia.org/wiki/Cauchy_stress_tensor#Stress_deviator_tensor

Parameters:

M - Tensor

Returns:

J2

computeVonMisesStress

public static double computeVonMisesStress(SymmetricMatrix3d sig)

Compute the Von Mises Stress criterion https://en.wikipedia.org/wiki/Von_Mises_yield_criterion The Von Mises Stress is equal to sqrt(3 J2), where J2 is the second invariant of the average deviatoric strain for the node.

Parameters:

sig - stress tensor

Returns:

Von Mises stress

computeVonMisesStrain

public static double computeVonMisesStrain(SymmetricMatrix3d strain)

Compute the Von Mises strain equivalent according to http://www.continuummechanics.org/vonmisesstress.html which is equivalent to https://dianafea.com/manuals/d944/Analys/node405.html The Von Mises Strain Equivalent is equal to sqrt(4/3 J2), where J2 is the second invariant of the average deviatoric strain for the node.

Parameters:

strain - strain tensor

Returns:

von Mises strain Equivalent