Displacement Finite Strain Element

Preliminaries

This element formulation is implemented in the class Marmot::Elements::DisplacementFiniteStrainULElement, which represents a displacement-based finite element considering large deformations. The formulation is valid for general three-dimensional geometries and can be combined with arbitrary material models. It supports both body forces and surface tractions.

This class inherits from the abstract base classes MarmotElement and MarmotGeometryElement, which provide common functionality for all finite elements in Marmot. In addition, it serves as a base class for specializations such as plane strain and axisymmetric elements.

Throughout the following derivations, Einstein summation convention is used, i.e., repeated indices imply summation. Uppercase indices refer to material coordinates, while lowercase indices refer to spatial coordinates. Accordingly, lowercase subscripts preceded by a comma denote spatial derivatives with respect to the current (deformed) configuration, e.g., \((\bullet)_{,i} = \partial (\bullet)/\partial x_i\), while uppercase subscripts preceded by a comma denote material derivatives with respect to the reference (undeformed) configuration, e.g., \((\bullet)_{,I} = \partial (\bullet)/\partial X_I\). Following this convention, the deformation gradient and its determinant is defined as

\[F_{iI} = x_{i,I}\ ,\ J = \text{det}\,F_{iI}.\]

Theory

This element computes the element contribution to the global residual vector \(\mathbf{r}_{Aj}\) and its derivative with respect to the nodal displacement vector \(\mathbf{q}_{Bk}\), i.e., the element contribution to the global stiffness matrix \(\partial \mathbf{r}_{Aj}/\partial \mathbf{q}_{Bk}\) for a displacement-based finite element considering large deformations. The global residual vector is derived by discretizing the weak form for linear momentum, given as

\[\mathbf{r}_{Aj} = \int_{V_0}\,\mathbf{N}_{A,i}\,\tau_{ij}\,dV_0 - \int_{V_0}\,\mathbf{N}_{A}\,f_{j}\,dV_0 - \int_{\bar{A}}\, \mathbf{N}_A\, \bar{t}_j \, d\bar{A},\]

with the Kirchhoff stress \(\tau_{ij} = J\,t_{ij}\), the body force \(f_j\), and the prescribed traction \(\bar{t}_j\) on the Neumann boundary \(\bar{A}\). The global stiffness matrix is computed as

\[\begin{split} \frac{\partial \mathbf{r}_{Aj}}{\partial \mathbf{q}_{Bk}} &= \mathbf{K}_{AjBk} = \int_{V_0} \mathbf{N}_{A,i}\, \frac{\partial \tau_{ij}}{\partial F_{kK}}\,\mathbf{N}_{B,K} - \mathbf{N}_{A,k}\,\mathbf{N}_{B,i}\,\tau_{ij}\,dV_0\, -\\ &- \int_{\bar{A}}\,\mathbf{N}_A\,\bar{t}_i\, \left(\delta_{ij}\delta_{lk} - \delta_{ik}\delta_{lj}\right)\, \mathbf{N}_{B,l}\, d\bar{A},\end{split}\]

where \(\delta_{ij}\) is the Kronecker delta and \(\partial \tau_{ij}/\partial F_{kK}\) the material tangent obtained from the respective material model.

Implementation

template<int nDim, int nNodes>
class DisplacementFiniteStrainULElement : public MarmotElement, public MarmotGeometryElement<nDim, nNodes>

Inheritance diagram for Marmot::Elements::DisplacementFiniteStrainULElement:

Collaboration diagram for Marmot::Elements::DisplacementFiniteStrainULElement:

Implementation of a displacement-based finite element for geometrically nonlinear analysis.

Template Parameters:

• nDim – Number of spatial dimensions (1, 2, 3)

• nNodes – Number of nodes of the element

Public Types

enum SectionType

SectionType is an enum class which involves the following cases:

Values:

enumerator PlaneStress

Plane stress section for 2D elements.

enumerator PlaneStrain

Plane strain section for 2D elements.

enumerator Solid

Solid (homogeneous) section for 1D, 2D and 3D elements.

using ParentGeometryElement = MarmotGeometryElement<nDim, nNodes>

Parent element class for geometry related operations as, e.g., shape functions.

using Material = MarmotMaterialFiniteStrain

Material class for finite strain material formulations.

using JacobianSized = typename ParentGeometryElement::JacobianSized

Sized matrix type used for Jacobian matrix (inherited from ParentGeometryElement)

using NSized = typename ParentGeometryElement::NSized

Sized matrix type used for shape function matrix (inherited from ParentGeometryElement)

using dNdXiSized = typename ParentGeometryElement::dNdXiSized

Sized matrix typed used for shape function derivatives (inherited from ParentGeometryElement)

using XiSized = typename ParentGeometryElement::XiSized

Sized vector type used for local coordinates (inherited from ParentGeometryElement)

using RhsSized = Eigen::Matrix<double, sizeLoadVector, 1>

Sized vector type used for the right hand side of the global equation system (negative element residual vector)

using KSizedMatrix = Eigen::Matrix<double, sizeLoadVector, sizeLoadVector>

Sized matrix type used for the element stiffness matrix.

using USizedVector = Eigen::Matrix<double, bsU, 1>

Sized vector type used for the element displacement vector.

using ForceSized = Eigen::Matrix<double, nDim, 1>

Sized vector type used for force vectors.

Public Functions

DisplacementFiniteStrainULElement(int elementID, Marmot::FiniteElement::Quadrature::IntegrationTypes integrationType, SectionType sectionType)

Constructor of the displacement-based finite strain element.

Parameters:

• elementID[in] – Element ID (label) of the element

• integrationType[in] – Integration type of the element

• sectionType[in] – Section type of the element

virtual int getNumberOfRequiredStateVars()

Get the total number of required state variables of the element.

virtual std::vector<std::vector<std::string>> getNodeFields()

Get the nodal fields of the element.

virtual std::vector<int> getDofIndicesPermutationPattern()

Get the permutation pattern of the element degrees of freedom.

inline virtual int getNNodes()

Get the number of nodes of the element.

Returns:

Number of nodes of the element

inline virtual int getNSpatialDimensions()

Get the number of spatial dimensions of the element.

Returns:

Number of spatial dimensions of the element

inline virtual int getNDofPerElement()

Get the number of degrees of freedom per node of the element.

Returns:

Number of degrees of freedom per node of the element

inline virtual std::string getElementShape()

Get the shape of the element.

Returns:

Shape of the element

virtual void assignStateVars(double *managedStateVars, int nStateVars)

Assign the state variable vector of the element.

Parameters:

• managedStateVars[in] – Pointer to the state variable vector of the element

• nStateVars[in] – Number of state variables of the element

virtual void assignProperty(const ElementProperties &MarmotElementProperty)

Assign the element properties of the element.

Parameters:

MarmotElementProperty[in] – Element properties

virtual void assignProperty(const MarmotMaterialSection &MarmotElementProperty)

Assign the material section of the element.

Parameters:

MarmotElementProperty[in] – Material section

virtual void assignNodeCoordinates(const double *coordinates)

Assign the nodal coordinates of the element.

Parameters:

coordinates[in] – Pointer to the nodal coordinates of the element

virtual void initializeYourself()

Initialize the element.

virtual void setInitialConditions(StateTypes state, const double *values)

Set the initial conditions of the element.

Parameters:

• state[in] – Type of the initial state

• values[in] – Pointer to the values defining the initial state

virtual void computeDistributedLoad(MarmotElement::DistributedLoadTypes loadType, double *P, double *K, const int elementFace, const double *load, const double *QTotal, const double *time, double dT)

Compute the contributions of distributed loads to the element residual vector and stiffness matrix.

For a given distributed load vector \(\bar{\boldsymbol{t}}^{(n+1)}\) at the current time step \(t^{(n+1)} = t^{(n)} + \Delta\,t\), compute the distributed load contribution to the negative element residual vector (right hand side of global newton) \(\int_\bar{A}\,\mathbf{N}_A\,\bar{t}_j\,d\bar{A}\) and the element stiffness matrix \(-\int_\bar{A}\,\mathbf{N}_{A}\,\bar{t}_i\left(\delta_{ij}\delta_{lk} - \delta_{ik}\delta_{lj}\right)\,\mathbf{N}_{B,l}\,d\bar{A}\).

Parameters:

• loadType[in] – Type of the distributed load, e.g., pressure or surface traction

• P[in, out] – Pointer to the element residual vector (right hand side of the global equation system)

• K[in, out] – Pointer to the element stiffness matrix

• elementFace[in] – Local face number of the element where the distributed load is applied

• load[in] – Pointer to the distributed load vector

• QTotal[in] – Pointer to the total element displacement vector at the current time step

• time[in] – Pointer to the time at the beginning of the current time step

• dT[in] – Length of the current time step

virtual void computeBodyForce(double *P, double *K, const double *load, const double *QTotal, const double *time, double dT)

Compute the contributions of body forces to the element residual vector and stiffness matrix.

For a given body force vector \(\mathbf{f}^{(n+1)}\) at the current time step \(t^{(n+1)} = t^{(n)} + \Delta\,t\), compute the body force contribution for the negative element residual vector (right hand side of global newton) \(\int_{V_0}\,\mathbf{N}_A\,f_j\,dV_0\). The stiffness matrix contribution is zero and thus not computed.

Parameters:

• P[in, out] – Pointer to the element residual vector (right hand side of the global equation system)

• K[in, out] – Pointer to the element stiffness matrix

• load[in] – Pointer to the body force vector

• QTotal[in] – Pointer to the total element displacement vector at the current time step

• time[in] – Pointer to the time at the beginning of the current time step

• dT[in] – Length of the current time step

virtual void computeYourself(const double *QTotal, const double *dQ, double *Pe, double *Ke, const double *time, double dT, double &pNewdT)

Compute the negative element residual vector (right hand side of global newton) and stiffness matrix.

For a given displacement \(\mathbf{q}^{(n+1)}\) at the current time step \(t^{(n+1)} = t^{(n)} + \Delta\,t\), compute the internal work contribution for the negative element residual vector (right hand side of global newton) \(-\int_{V_0}\,\mathbf{N}_{A,i}\,\tau_{ij}\,dV_0\) and the element stiffness matrix \(\int_{V_0}\,\mathbf{N}_{A,i}\,\frac{\partial \tau_{ij}}{\partial F_{kK}}\,\mathbf{N}_{B,K}\,-\,\mathbf{N}_{A,k\,}\mathbf{N}_{B,i}\,\tau_{ij}\,dV_0\).

Parameters:

• QTotal[in] – Pointer to the total element displacement vector at the current time step

• dQ[in] – Pointer to the increment of the element displacement vector at the current time step

• Pe[in, out] – Pointer to the negative element residual vector (right hand side of global newton)

• Ke[in, out] – Pointer to the element stiffness matrix

• time[in] – Pointer to the time at the beginning of the current time step

• dT[in] – Length of the current time step

• pNewdT[in, out] – Suggested length of the next time step

virtual StateView getStateView(const std::string &stateName, int qpNumber)

Get a view to a state variable at a specific quadrature point of the element.

Parameters:

• stateName[in] – Name of the state variable

• qpNumber[in] – Number of the quadrature point where the state variable is stored

Returns:

View to the requested state variable at the specified quadrature point

virtual std::vector<double> getCoordinatesAtCenter()

Get the coordinates of the element center.

Returns:

Coordinates of the element center

virtual std::vector<std::vector<double>> getCoordinatesAtQuadraturePoints()

Get the coordinates of all quadrature points of the element.

Returns:

Coordinates of all quadrature points of the element

virtual int getNumberOfQuadraturePoints()

Get the number of quadrature points of the element.

Returns:

Number of quadrature points of the element

Public Members

Eigen::Map<const Eigen::VectorXd> elementProperties

Element properties as provided in the input file.

const int elLabel

Element label (ID)

const SectionType sectionType

Section type of the element.

bool hasEigenDeformation

Boolean for indicating whether initial deformation is considered or not.

std::vector<QuadraturePoint> qps

List of quadrature points of the element.

Public Static Attributes

static int nDofPerNodeU = nDim

Displacement degrees of freedom per node.

static int nCoordinates = nNodes * nDim

Total number of coordinates of the element.

static int bsU = nNodes * nDofPerNodeU

Block size of element stiffness matrix and load vector for displacement field U.

static int sizeLoadVector = bsU

Size of element stiffness matrix and load vector.

static int idxU = 0

Starting index of displacement field U in element stiffness matrix and load vector.

struct QuadraturePoint

Structure for storing quadrature point related information.

Public Functions

inline int getNumberOfRequiredStateVarsQuadraturePointOnly()

Get number of required state variables at the quadrature point only.

Returns:

Number of required state variables at the quadrature point only (without material state variables)

inline int getNumberOfRequiredStateVars()

Get total number of required state variables at the quadrature point.

Returns:

Total number of required state variables at the quadrature point (including material state variables)

inline void assignStateVars(double *stateVars, int nStateVars)

Assign the state variable vector at the quadrature point.

Parameters:

• stateVars[in] – Pointer to the state variable vector at the quadrature point

• nStateVars[in] – Number of state variables at the quadrature point

inline QuadraturePoint(XiSized xi, double weight)

Constructor of the quadrature point.

Note

The shape function derivatives w.r.t. material (undeformed) coordinates and the determinant of the undeformed Jacobian times quadrature weight are initialized with zero values.

Parameters:

• xi[in] – Local coordinates of the quadrature point

• weight[in] – Weight of the quadrature point

Public Members

const XiSized xi

Local coordinates of the quadrature point

const double weight

Weight of the quadrature point

dNdXiSized dNdX

Shape function derivatives w.r.t. material (undeformed) coordinates evaluated at the quadrature point

double J0xW

Determinant of the undeformed Jacobian times quadrature weight

std::unique_ptr<QPStateVarManager> managedStateVars

Managed state variables at the quadrature point.

std::unique_ptr<Material> material

Material at the quadrature point.

class QPStateVarManager : public MarmotStateVarVectorManager

Inheritance diagram for Marmot::Elements::DisplacementFiniteStrainULElement::QuadraturePoint::QPStateVarManager:

Collaboration diagram for Marmot::Elements::DisplacementFiniteStrainULElement::QuadraturePoint::QPStateVarManager:

Manager class for handling state variables at the quadrature point.

Public Functions

inline QPStateVarManager(double *theStateVarVector, int nStateVars)

Constructor of the state variable manager at the quadrature point.

Parameters:

• theStateVarVector[in] – Pointer to the state variable vector at the quadrature point

• nStateVars[in] – Number of state variables at the quadrature point

Public Members

Eigen::Map<Marmot::Vector9d> stress

Stress tensor at the quadrature point

double &F0_XX

Deformation gradient component XX for prescribing an initial deformation state

double &F0_YY

Deformation gradient component YY for prescribing an initial deformation state

double &F0_ZZ

Deformation gradient component ZZ for prescribing an initial deformation state

Eigen::Map<Eigen::VectorXd> materialStateVars

Material state variables at the quadrature point

Public Static Functions

static inline int getNumberOfRequiredStateVarsQuadraturePointOnly()

Get number of required state variables at the quadrature point only (without material state variables)

Private Static Attributes

static const auto layout = makeLayout( {{ .name = "stress", .length = 9 },{ .name = "F0 XX", .length = 1 },{ .name = "F0 YY", .length = 1 },{ .name = "F0 ZZ", .length = 1 },{ .name = "begin ofmaterialstate", .length = 0 },} )

Layout of the state variable vector at the quadrature point.