General Gradient-Enhanced Displacement Finite Element

Theory

This element implements a small strain gradient enhanced displacement based Finite Element formulation. The formulation couples the standard displacement field with one or multiple nonlocal scalar fields.

In addition to the standard static equilibrium equation, one or multiple nonlocal scalar fields are solved for.

\[\nabla \sig + \boldsymbol f = \boldsymbol 0\]

Here \(\sig\) is the Cauchy stress tensor and \(\boldsymbol f\) the body force vector per unit volume. Each nonlocal variable introduces an additional balance equation, which is solved simultaneously with the static equilibrium equation. The general balance equation for the i-th nonlocal variable \(\knl_i\) is defined as

\[\knl_i - \nabla ( c( \knl_i )\, \nabla \knl_i ) = \kl_i (\eps,\, \knl_i )\, ,\]

where \(\knl\) is the nonlocal variable, \(\kl\) the local equivalent, \(c( \knl_i )\) the gradient interaction parameter, and \(\kl_i\) is the local driving variable for the nonlocal field. The element formulation supports a variable length scale parameter which can be used to model damage dependent interactions and may depend on the state of the nonlocal field: \(c = c(\knl)\). This implementation is versatile and is suited for both implicit gradient enhanced models and phase field formulations to simulate failure and strain localization. Depending on the specific material model, this field is governed either by the screened Poisson equation for implicit gradient regularisation,

\[\knl - c \, \nabla^2 \knl = \kl(\eps)\, ,\]

or by the phase field crack evolution equation:

\[d - l^2 \, \nabla^2 d = s(d, \mathcal{Y})\, .\]

For the screened Poisson equation \(\knl\) is the nonlocal equivalent strain-like variable, \(\kl\) the local equivalent, and \(c=l^2\) the gradient length scale parameter. In the context of phase field models, \(d\) represents the phase field variable, \(l\) is the internal length scale, and \(s(d, \mathcal{Y})\) is a local source function driven by the thermodynamic driving force \(\mathcal{Y}\). Homogeneous Neumann boundary conditions are assumed for the nonlocal scalar fields on the element boundary, i.e., \((\nabla \knl \boldsymbol)^T\, n = 0\) or \((\nabla d)^T \, \boldsymbol n = 0\), where \(\boldsymbol n\) is the outward unit normal.

The nodal unknown vector uses a field wise structure:

\[\mathbf{q} = [\mathbf{\qu}, \mathbf{\qk}]^\mathsf{T}\]

The incremental strain is computed from the displacement increment via the kinematic matrix \(\mathbf{B}\). The nonlocal field is interpolated using the shape function for the nonlocal field \(\Nk\) which may be configured to be an order lower than the displacement shape functions.

\[\Delta \eps = \mathbf{B}\, \Delta \mathbf{\qu}, \qquad \knl = \Nk\, \mathbf{\qk} % \knl = \Nk_A\, \qk_A,\, ,\]

Here \(\mathbf{B}\) is the strain-displacement matrix, \(\Delta \mathbf{\qu}\) the nodal displacement increment, and \(\Delta \eps\) is the linearized strain tensor in Voigt notation.

The discretized weak form of the governing equations is given by

\[ \begin{align}\begin{aligned}\delta \qu_{Aj} \left( \int_{V} \mathbf{N}_{A,i}\, \sig_{ij} \, dV -\int_{V} \mathbf{N}_{A}\, f_j \, dV -\int_{S_\mathrm{t}} \mathbf{N}_A\, \bar{t}_j \, d S_\mathrm{t} \right) = \delta \qu_{Aj}\, \ru_{Aj} = 0\, ,\\\delta \qk_A \int_V \mathbf{N}_A \, \knl + c \, \mathbf{N}_{A,i}\, \knl_{,i} - \mathbf{N}_A\, \kl \, dV = \delta \qk_A \, \rk_A = 0\, ,\end{aligned}\end{align} \]

where \(\delta \qu_{Aj}\) and \(\delta \qk_A\) are the virtual nodal variations of the displacement and nonlocal field, respectively. \(\ru_{Aj}\) and \(\rk_A\) denote the nodal residuals associated with the respective degrees of freedom. Uppercase subscripts, for example \((\bullet)_A\), represent nodal indices that denote specific points in the discretised domain, while lowercase subscripts, for example \((\bullet)_i\), denote the spatial dimensions. Internal forces and consistent tangents are evaluated by Gauss quadrature:

\[ \begin{align}\begin{aligned}\begin{split}\mathbf{K}_e = \begin{bmatrix} \mathbf{K}_{uu} & \mathbf{K}_{uk} \\ \mathbf{K}_{ku} & \mathbf{K}_{kk} \end{bmatrix},\qquad \mathbf{\fuint} = \sum_{qp} \mathbf{B}^\mathsf{T}\, \sig\, J_0\, w_{qp}\, ,\end{split}\\\mathbf{\fk} = \sum_{qp} \left (\mathbf{\Nk}^\mathsf{T}\, \knl\, + c\, \partial_\mathbf{x} \mathbf{\Nk}^\mathsf{T}\, \partial_\mathbf{x} \mathbf{\Nk}\, \mathbf{\qk} - \mathbf{\Nk}^\mathsf{T}\, \kl \right ) J_0\, w_{qp} \, .\end{aligned}\end{align} \]

Here \(\mathbf{K}_e\) the element tangent stiffness in block structure, \(J_0 = \det \mathbf{J}\) the Jacobian determinant and \(w_{qp}\) the quadrature weight, \(\mathbf{\fuint}\) the element internal force vector for the displacement degrees of freedom, \(\mathbf{\fk}\) the element force vector for the nonlocal degrees of freedom, \(\sig\) the Cauchy stress tensor in Voigt notation, and \(\sum_{qp}\) denotes summation over quadrature points.

The stiffness submatrices \(\mathbf{K}_{uu}\), \(\mathbf{K}_{u\knl}\), \(\mathbf{K}_{\knl u}\) and \(\mathbf{K}_{\knl\knl}\) are evaluated using the following expressions:

\[\begin{split}\begin{aligned} & \mathbf{K}_{uu} = \sum_{qp} \mathbf{B}^\mathsf{T} \frac{\partial \mathbf{\sig}}{\partial \mathbf{\eps}} \mathbf{B}\, J_0\, w_{qp}, \\ & \mathbf{K}_{u\knl} = \sum_{qp} \mathbf{B}^\mathsf{T} \frac{\partial \mathbf{\sig}}{\partial \knl} \mathbf{\Nk}\, J_0\, w_{qp}, \\ & \mathbf{K}_{\knl u} = -\sum_{qp} \mathbf{\Nk}^\mathsf{T} \frac{\partial \kl}{\partial \mathbf{\eps}} \mathbf{B}\, J_0\, w_{qp}, \\ & \mathbf{K}_{\knl\knl} = \sum_{qp} \left( \mathbf{\Nk}^\mathsf{T}\, \mathbf{\Nk} + c \, \partial_\mathbf{x} \mathbf{\Nk}^\mathsf{T} \, \partial_\mathbf{x} \mathbf{\Nk} + \frac{\partial c}{\partial \knl} \, \partial_\mathbf{x} \mathbf{\Nk}^\mathsf{T} \, \partial_\mathbf{x} \mathbf{\Nk}\, \qk \, \mathbf{\Nk} - \mathbf{\Nk}^\mathsf{T}\, \frac{\partial \kl}{\partial \knl} \right) J_0\, w_{qp}\, . \end{aligned}\end{split}\]

The derivative \(\partial c / \partial \knl\) is only non-zero if a variable length scale is used, and the derivative \(\partial \kl / \partial \knl\) is only required in phase field formulations where it represents the derivative of the source function with respect to the crack phase field \(\frac{\partial s(d, \mathcal{Y})}{\partial d}\).

The body-force vector is obtained using

\[\mathbf{\fextb} = \sum_{qp} \mathbf{N}^\mathsf{T} \mathbf{f}\, J_0 w\, ,\]

and surface pressures on a boundary face \(\Gamma_e\) are integrated as

\[\mathbf{\fextp} = - \int_{\Gamma_e} p \, \mathbf{N}^\mathsf{T} \, \mathbf{n} \, \mathrm{d}\Gamma.\]

where \(\Gamma_e\) is the element boundary, \(\mathbf{n}\) its outward unit normal, \(\mathbf{N}\) the shape-function matrix of the displacement field, \(p\) the pressure magnitude, and \(\mathbf{f}\) the body-force vector per unit volume. The minus sign for pressure follows an outward-normal convention.

For 2D and 1D problems, the integration measure \(J_0 w\) is scaled by thickness \(t\) or cross-sectional area \(A\), respectively: \(J_0 w \leftarrow J_0 w\, t\) (2D), \(J_0 w \leftarrow J_0 w\, A\) (1D). Here, \(t\) denotes thickness and \(A\) denotes cross-sectional area.

Constitutive updates are carried out per quadrature point using a Marmot material model of the following types:

• 3D solid

• Plane stress

• Plane strain

• 1D uniaxial stress

Each quadrature point stores stress, strain and a material state vector; the element accumulates the history incrementally and may suggest a reduced time step if required by the material.

Implementation

template<int nDim, int nNodes, int nNonlocalVariables = 1, int nNonLocalNodes = nNodes>
class GeneralGradientEnhancedDisplacementFiniteElement : public MarmotElement

Inheritance diagram for Marmot::Elements::GeneralGradientEnhancedDisplacementFiniteElement:

Collaboration diagram for Marmot::Elements::GeneralGradientEnhancedDisplacementFiniteElement:

General Gradient Enhanced Displacement-based finite element template.

Uses linearized kinematics (small strains) and supports 1D, 2D, and 3D formulations via a section assumption. Holds quadrature-point state and delegates constitutive updates to Marmot materials while assembling element residuals and tangents, taking non-local effects into account.

Template Parameters:

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

• nNodes – Number of element nodes for the displacement field.

• nNonlocalVariables – Number of nonlocal variables.

• nNonLocalNodes – Number of nonlocal element nodes. Usually the order of the nonlocal interpolation is similar or one order lower than the local interpolation, so nNonLocalNodes <= nNodes.

Public Types

enum SectionType

Kinematic section assumption used by the element.

Controls which constitutive call is performed at quadrature points.

• PlaneStress

• PlaneStrain

• Solid

Values:

enumerator PlaneStress

enumerator PlaneStrain

enumerator Solid

using ParentGeometryElement = MarmotGeometryElement<nDim, nNodes>

using JacobianSized = typename ParentGeometryElement::JacobianSized

using NSized = typename ParentGeometryElement::NSized

using dNdXiSized = typename ParentGeometryElement::dNdXiSized

using BSized = typename ParentGeometryElement::BSized

using XiSized = typename ParentGeometryElement::XiSized

using RhsSized = Matrix<double, sizeLoadVector, 1>

using KeSizedMatrix = Matrix<double, sizeLoadVector, sizeLoadVector>

using USizedVector = Matrix<double, sizeDoFU, 1>

using KSizedVector = Matrix<double, sizeDoFK, 1>

using CSized = Matrix<double, ParentGeometryElement::voigtSize, ParentGeometryElement::voigtSize>

using Voigt = Matrix<double, ParentGeometryElement::voigtSize, 1>

using ParentGeometryElementK = MarmotGeometryElement<nDim, nNonLocalNodes>

using JacobianSizedK = typename ParentGeometryElementK::JacobianSized

using NSizedK = typename ParentGeometryElementK::NSized

using dNdXiSizedK = typename ParentGeometryElementK::dNdXiSized

using BSizedK = typename ParentGeometryElementK::BSized

Public Functions

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

Construct element with ID, quadrature rule and section assumption.

Parameters:

• elementID – Unique element label.

• integrationType – Integration (quadrature) rule.

• sectionType – Section assumption (1D/2D/3D).

virtual int getNumberOfRequiredStateVars()

Total number of required state variables for this element (sum over all quadrature points).

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

Node-level fields exposed by the element. Returns [“displacement”, “nonlocal field”, “strain symmetric”] for the respective nodes.

virtual std::vector<int> getDofIndicesPermutationPattern()

The permutation pattern for the residual vector and the stiffness matrix to aggregate all entries in order to resemble the defined fields nodewise.

inline virtual int getNNodes()

Number of nodes of this element type.

inline virtual int getNSpatialDimensions()

Number of spatial dimensions.

inline virtual int getNDofPerElement()

Number of degrees of freedom per element.

inline virtual std::string getElementShape()

Geometric shape of the element (as reported by the local geometry element).

virtual void assignStateVars(double *stateVars, int nStateVars)

Map the provided element state vector to all quadrature points.

virtual void assignProperty(const ElementProperties &marmotElementProperty)

Assign element properties (e.g., thickness in 2D, area in 1D).

virtual void assignProperty(const MarmotMaterialSection &marmotElementProperty)

Assign material section and instantiate per-quadrature-point materials.

virtual void assignNodeCoordinates(const double *coordinates)

Provide nodal coordinates to local and non-local geometry elements.

virtual void initializeYourself()

Precompute geometry-related quantities at quadrature points (B, detJ, J0xW, N_K, dNdX_K).

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

Initialize state or materials.

Parameters:

• state – MarmotMaterialInitialization, GeostaticStress or MarmotMaterialStateVars.

• values – For GeostaticStress: [sigmaY(z1), y1, sigmaY(z2), y2, kx, kz].

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

Assemble distributed surface loads on a boundary face.

Pressure contributions are integrated on the boundary \(\Gamma_e\):

\[\mathbf{\fextp} = - \int_{\Gamma_e} p\, \mathbf{N}^\mathsf{T} \mathbf{n}\, \mathrm{d}\Gamma\, . \]

Parameters:

• loadType – Pressure or SurfaceTraction.

• P – Element RHS contribution (accumulated).

• K – Optional stiffness contribution (unused).

• elementFace – Boundary face index.

• load – Pressure magnitude or traction vector (size nDim).

• QTotal – Total DOF vector (unused).

• time – Current time data forwarded to materials.

• dT – Time increment.

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

Assemble body force contribution.

Integrates \(\mathbf{\fextb} = \int_{\Omega_e} \mathbf{N}^\mathsf{T} \mathbf{f}\, \mathrm{d}\Omega\).

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

Compute internal force and consistent tangent stiffness.

Uses the small-strain relation \(\Delta \eps = \mathbf{B}\, \Delta \mathbf{\qu}\), interpolates the non-local variables as \( \knl = \Nk\, \mathbf{\qk}\). Internal forces and consistent tangents are evaluated by Gauss quadrature:

\[\begin{split}\mathbf{K}_e = \begin{bmatrix} \mathbf{K}_{uu} & \mathbf{K}_{uk} \\ \mathbf{K}_{ku} & \mathbf{K}_{kk} \end{bmatrix},\qquad \mathbf{\fuint} = \sum_{qp} \mathbf{B}^\mathsf{T}\, \sig\, J_0\, w_{qp}\, . \end{split}\]

\[\mathbf{\fk} = \sum_{qp} \left (\mathbf{\Nk}^\mathsf{T}\, \knl\, + c\, \partial_\mathbf{x} \mathbf{\Nk}^\mathsf{T}\,\partial_\mathbf{x} \mathbf{\Nk}\, \mathbf{\qk} - \mathbf{\Nk}^\mathsf{T}\, \kl \right ) J_0\, w_{qp} \, . \]

The stiffness submatrices are evaluated using the following expressions:

\[\mathbf{K}_{uu} = \sum_{qp} \mathbf{B}^\mathsf{T} \frac{\partial \mathbf{\sig}}{\partial \mathbf{\eps}} \mathbf{B}\, J_0\, w_{qp},\quad \mathbf{K}_{u\knl} = \sum_{qp} \mathbf{B}^\mathsf{T} \frac{\partial \mathbf{\sig}}{\partial \knl} \mathbf{\Nk}\, J_0\, w_{qp},\quad \mathbf{K}_{\knl u} = -\sum_{qp} \mathbf{\Nk}^\mathsf{T} \frac{\partial \kl}{\partial \mathbf{\eps}} \mathbf{B}\, J_0\, w_{qp}, \]

\[\mathbf{K}_{\knl\knl} = \sum_{qp} \left( \mathbf{\Nk}^\mathsf{T}\, \mathbf{\Nk} + c \, \partial_\mathbf{x} \mathbf{\Nk}^\mathsf{T} \, \partial_\mathbf{x} \mathbf{\Nk} + \frac{\partial c}{\partial \knl} \, \partial_\mathbf{x} \mathbf{\Nk}^\mathsf{T} \, \partial_\mathbf{x} \mathbf{\Nk}\, \qk \, \mathbf{\Nk} - \mathbf{\Nk}^\mathsf{T}\, \frac{\partial \kl}{\partial \knl} \right) J_0\, w_{qp}\, . \]

If pNewdT<1, the routine returns early to signal time step reduction.

Parameters:

• QTotal – Total displacement vector in field-wise format: \(\mathbf{q} = [\mathbf{\qu}, \mathbf{\qk}]^\mathsf{T}\).

• dQ – Incremental displacement.

• Pe – Internal force vector (accumulated).

• Ke – Tangent stiffness matrix (accumulated).

• time – Time data forwarded to materials.

• dT – Time increment.

• pNewdT – Suggested scaling of dT by the material; if reduced (<1), the routine returns early.

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

Compute internal force without tangents.

Uses the small-strain relation \(\Delta \eps = \mathbf{B}\, \Delta \mathbf{\qu}\), interpolates the non-local variables as \( \knl = \Nk\, \mathbf{\qk}\). Internal forces are evaluated by Gauss quadrature: \mathbf{\fk} = \sum_{qp} \left (\mathbf{\Nk}^\mathsf{T}\, \knl\, + c\, \partial_\mathbf{x} \mathbf{\Nk}^\mathsf{T}\,\partial_\mathbf{x} \mathbf{\Nk}\, \mathbf{\qk} - \mathbf{\Nk}^\mathsf{T}\, \kl \right ) J_0\, w_{qp} \, . \f] If pNewdT<1, the routine returns early to signal time step reduction.

Parameters:

• QTotal – Total displacement vector in field-wise format: \(\mathbf{q} = [\mathbf{\qu}, \mathbf{\qk}]^\mathsf{T}\).

• dQ – Incremental displacement.

• Pe – Internal force vector (accumulated).

• Ke – Tangent stiffness matrix (accumulated).

• time – Time data forwarded to materials.

• dT – Time increment.

• pNewdT – Suggested scaling of dT by the material; if reduced (<1), the routine returns early.

virtual void computeConsistentInertia(double *M)

Compute consistent mass matrix using material density.

\(\mathbf{M}_e = \sum_{qp} \rho\, \mathbf{N}^\mathsf{T}\mathbf{N}\, J_0 w\).

virtual void computeLumpedInertia(double *M)

Compute the lumped (diagonal) mass matrix.

Uses the manifold-based lumping scheme according to Yang et al. (2017) “A rigorous and unified mass lumping scheme for higher-order elements”, CMAME. The lumped mass entries are computed using a weighted shape function \(\hat{N} = \tfrac{1}{2}N + \tfrac{1}{2}N_\mathrm{lin}\), where \(N\) is the high-order shape function and \(N_\mathrm{lin}\) is the corresponding linear (corner-node) shape function on the same element.

virtual void computeCriticalTimeStepForExplicitDynamics(double &criticalTimeStep, const double *QTotal)

Compute the critical time step for explicit dynamics based on the dilatational wave speed and the element size.

Uses the formula \(\Delta t_\mathrm{crit} = \frac{l_\mathrm{min}}{c_\mathrm{dil}}\), where \(l_\mathrm{min}\) is the minimum characteristic element length and \(c_\mathrm{dil}\) is the dilatational wave speed computed from the material properties at the quadrature points. The minimum time step across all quadrature points is returned.

Parameters:

criticalTimeStep – Output parameter for the computed critical time step.

virtual void computeInternalEnergy(double &internalEnergy)

Compute the internal energy of the element by summing the strain energy contributions from all quadrature points.

Parameters:

internalEnergy – Output parameter for the computed internal energy.

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

Access a named state view at a quadrature point.

Note

Using “sdv” returns the raw material state vector and is deprecated.

virtual std::vector<double> getCoordinatesAtCenter()

Get physical coordinates at the element center.

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

Get physical coordinates at each quadrature point.

virtual int getNumberOfQuadraturePoints()

Number of quadrature points of this element.

Public Members

MarmotGeometryElement<nDim, nNodes> localGeometryElement

MarmotGeometryElement<nDim, nNonLocalNodes> nonLocalGeometryElement

Map<const VectorXd> elementProperties

Element-level properties (e.g., thickness for 2D, area for 1D).

const int elLabel

Element label (ID) used for logging and material creation.

const SectionType sectionType

Section assumption applied by this element instance.

std::vector<QuadraturePoint> qps

Quadrature points owned by the element (one per integration point).

Public Static Attributes

static int nDofPerNodeU = nDim

static int nDofPerNodeK = nNonlocalVariables

static int sizeLoadVector = nNodes * nDofPerNodeU + nNonLocalNodes * nDofPerNodeK

static int nCoordinates = nNodes * nDim

static int sizeDoFU = nNodes * nDofPerNodeU

static int sizeDoFK = nNonLocalNodes * nDofPerNodeK

struct QuadraturePoint

Data and state associated with a quadrature point.

Holds parent coordinates, integration weight, Jacobian determinant, shape functions N and their gradients dNdX, kinematic (strain-displacement) B-matrices for both local and non-local evaluations, and a material instance with managed state variables.

Public Functions

inline int getNumberOfRequiredStateVarsQuadraturePointOnly()

inline int getNumberOfRequiredStateVars()

inline void assignStateVars(double *stateVars, int nStateVars)

inline QuadraturePoint(XiSized xi, double weight)

Public Members

const XiSized xi

const double weight

double detJ

double J0xW

NSized N

dNdXiSized dNdX

BSized B

NSizedK N_K

dNdXiSizedK dNdX_K

BSizedK B_K

std::unique_ptr<QPStateVarManager> managedStateVars

std::unique_ptr<MarmotMaterialGeneralGradientEnhancedHypoElastic<nNonlocalVariables>> material

class QPStateVarManager : public MarmotStateVarVectorManager

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

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

Manager for per-quadrature-point state variables.

Provides named accessors to stress \(\sigma\), strain \(\varepsilon\) and the material state vector. The layout is [stress(6), strain(6), begin of material state(…)] in 3D Voigt notation.

Public Functions

inline QPStateVarManager(double *theStateVarVector, int nStateVars)

Public Members

mVector6d stress

mVector6d strain

double &totalStrainEnergy

double &elasticStrainEnergy

double &dissipation

Eigen::Map<Eigen::VectorXd> materialStateVars

Public Static Functions

static inline int getNumberOfRequiredStateVarsQuadraturePointOnly()

Private Static Attributes

static const auto layout