Nonlinear Mechanics and Fields

div(r::Union{VectorField,TensorField})

Solves the divergence of the vector field or tensor field r. An alternative way to solve div is to use ∇ as a differencial operator.

Return: ScalarField or VectorField

Types:

• r: VectorField or TensorField

3D Examples (assumes problem is set as in the ∇ doc setup)

# Assumes a 3D mesh with physical group "body".

# 1) Divergence of a 3D vector field → ScalarField
v1(X,Y,Z) = X
v2(X,Y,Z) = Y
v3(X,Y,Z) = Z
v = vectorField(problem, [field("body", fx=v1, fy=v2, fz=v3)])
D1 = div(v)
D2 = ∇ ⋅ v

# 2) Divergence of a 3D tensor field → VectorField
fsz(X, Y, Z) = 10 - Z
S = tensorField(problem, [field("body", fz=fsz)])
b1 = -div(S)
b2 = -S ⋅ ∇

IIPiolaKirchhoff(energy, C::TensorField, params)

Computes the nodal Second Piola–Kirchhoff (II PK) stress tensor field from a given nodal Right Cauchy–Green deformation tensor field.

For each time step and each node, the II PK stress is evaluated as

S = ∂ψ(C, params) / ∂C

using the provided free energy function energy and material parameters params. The stress evaluation is performed directly at the nodal values of C, without Gauss-point integration or spatial averaging.

This function is intended for post-processing, visualization, and diagnostic purposes. The resulting stress field is energy-consistent with the constitutive model but is not suitable for assembling internal forces or tangents, which require Gauss-point evaluation.

Arguments:

• energy::Function: Free energy density ψ(C, params).
• C::TensorField: Nodal Right Cauchy–Green deformation tensor field.
• params::NamedTuple: Material parameters passed to energy.

Requirements:

• The problem must be three-dimensional (dim = pdim = 3).
• The input field C must be defined nodally.

Returns:

• TensorField: Nodal Second Piola–Kirchhoff stress field with the same time steps and model as C.

Notes:

• The stress is computed independently at each node and time step.
• No interpolation, integration, or smoothing is performed.
• For use in weak forms or Newton iterations, Gauss-point evaluation should be used instead.

curl(r::VectorField)

Solves the rotation of the vector field r. An alternative way to solve curl is to use ∇ as a differencial operator.

Return: VectorField

Types:

• r: VectorField

3D Example (assumes problem is set as in the ∇ doc setup)

# Assumes a 3D mesh with physical group "body".
vx(X, Y, Z) = 0
vy(X, Y, Z) = X
vz(X, Y, Z) = 0
v = vectorField(problem, [field("body", fx=vx, fy=vy, fz=vz)])
D1 = curl(v)
D2 = ∇ × v

externalTangentFollower(
    problem::Problem,
    loads::Vector{BoundaryCondition};
    F::TensorField
)

Assembles the external tangent stiffness matrix associated with follower (configuration-dependent) surface loads for a finite deformation analysis.

The function computes the linearization of the external force vector with respect to the displacement field when the applied tractions rotate and/or stretch together with the deforming body.

Arguments

• problem::Problem Finite element problem definition.

• loads::Vector{BoundaryCondition} Boundary conditions defining surface tractions or pressures acting on the current configuration. Supported load types:

  • surface traction components (fx, fy, fz),
  • surface pressure (p, only for 3D solids).

  Each boundary condition is interpreted as a follower load, i.e. its direction and/or magnitude may depend on the deformation.

Keyword arguments

• F::TensorField Nodal deformation gradient field. Each node stores the full deformation gradient F, which is interpolated to Gauss points and used to:
  • rotate the applied tractions,
  • compute the Jacobian and inverse deformation gradient,
  • form the consistent external tangent matrix.

Mathematical formulation

The external force vector for follower tractions may be written as

f_ext = ∫_Γ Nᵀ · t(F) dΓ

where the traction vector depends on the deformation gradient, typically through a relation of the form

t = J F⁻ᵀ t₀    or    t = F · t_loc

The present implementation linearizes this expression consistently, leading to an external tangent contribution

K_ext = ∂f_ext / ∂u

which includes:

• the variation of the Jacobian and inverse deformation gradient,
• the variation of the traction due to rotation/stretching with F.

Both pressure loads and vector-valued surface tractions are supported.

Notes

• This function assembles only the tangent contribution of follower loads. The corresponding external force vector must be assembled separately, e.g. using loadVector with deformation-dependent tractions.

• Dead loads (loads fixed in the reference configuration) must not be passed to this function. For dead loads, the external tangent is identically zero.

• The formulation assumes a Total Lagrange–type description and is intended for large-rotation and large-deformation problems.

• For pressure loads, the reference normal vector is obtained from the undeformed geometry and appropriately transformed.

Returns

• SystemMatrix Sparse global external tangent stiffness matrix of size (ndof × ndof).

Typical usage

Kext = externalTangentFollower(
    problem,
    loads;
    F = Ffield
)

K = Kmat + Kgeo - Kext

See also

• loadVector
• internalForceVector
• materialTangentMatrix
• initialStressVector

grad(r::Union{ScalarField,VectorField})

Solves the gradient of the scalar field or vector field r. An alternative way to solve grad is to use ∇ as a differencial operator.

Return: VectorField or TensorField

Types:

• r: ScalarField or VectorField

3D Examples

# Assumes a 3D mesh with physical group "body".

# 1) Gradient of a 3D scalar field → VectorField
f(X,Y,Z) = X^2 + Y*Z
S = scalarField(problem, [field("body", f=f)])
G1 = grad(S)
G2 = ∇(S)

# 2) Gradient of a 3D vector field → TensorField
vx(X,Y,Z) = X
vy(X,Y,Z) = Y
vz(X,Y,Z) = Z
V = vectorField(problem, [field("body", fx=vx, fy=vy, fz=vz)])
T1 = grad(V)
T2 = V ∘ ∇

initialStressMatrix(
    problem::Problem;
    stress::TensorField,
    energy::Function,
    C::TensorField,
    params
)

Assembles the geometric (initial stress) stiffness matrix associated with a given stress field, for finite deformation analysis in a Total Lagrange–type formulation.

The function computes the matrix

K_geo = ∫_Ω G(S) dΩ

where G(S) represents the geometric stiffness contribution induced by the supplied stress tensor field stress. The stress measure itself (P, S, or σ) is not interpreted by the routine; it is treated purely as a second-order tensor field whose contraction with displacement gradients yields the geometric stiffness.

Arguments

• problem::Problem Finite element problem definition. Must satisfy:
  • problem.dim == problem.pdim

Keyword arguments

• stress::TensorField Nodal stress tensor field.

  Each node stores a full dim×dim tensor, which is interpolated to Gauss points using standard Lagrange shape functions. The physical meaning of the tensor is left to the caller; typical choices include:

  • First Piola–Kirchhoff stress P,
  • Second Piola–Kirchhoff stress S,
  • Cauchy stress σ,

  provided that the chosen stress measure is consistent with the kinematic description used elsewhere in the formulation.

• energy::Function Free-energy density ψ(C, params) of a hyperelastic material. If provided, the stress field is computed at Gauss points via Tensors.jl from the supplied C and stress must be omitted.

• C::TensorField Nodal field of the right Cauchy–Green tensor C = FᵀF. Required when using the energy-based path.

• params Optional parameter container passed through to energy.

Mathematical formulation

At each Gauss point, the stress tensor S_gp is obtained by interpolation from nodal values. The geometric stiffness contribution is computed as

g_ab = (∇N_a)ᵀ · S_gp · (∇N_b)

and distributed to the displacement components as

(K_geo)_ai,bi = g_ab

with no coupling between different displacement directions (i = i only).

The element contribution reads

K_e = ∫ (∇N)ᵀ · S · (∇N) dΩ

and is assembled into the global matrix.

Notes

• This function assembles only the geometric (stress-dependent) part of the tangent stiffness matrix. The material (constitutive) tangent must be assembled separately, e.g. via materialTangentMatrix.

• The routine is stress-measure agnostic: it does not enforce symmetry or specific push-forward/pull-back operations. Ensuring energetic and kinematic consistency between the stress field and the rest of the formulation is the responsibility of the caller.

• If energy is supplied, C is interpolated to Gauss points and the stress is obtained from the free-energy function using Tensors.jl. This enables arbitrary hyperelastic laws without explicitly providing a stress field.

• The formulation corresponds to the classical initial stress stiffness encountered in large-deformation and stability (buckling) analyses.

Returns

• SystemMatrix Sparse global geometric stiffness matrix of size (ndof × ndof).

Typical usage

Kgeo = initialStressMatrix(
    problem;
    stress = Sfield   # TensorField (P, S, or σ)
)

K = Kmat + Kgeo

See also

• materialTangentMatrix
• internalForceVector
• TensorField

internalForceVector(
    problem::Problem;
    P::TensorField,
    F::TensorField,
    energy::Function,
    params
)

Assembles the internal force vector associated with a given first-order stress tensor field P, using a finite deformation formulation. Alternatively, P can be computed from a free-energy density energy (hyperelastic formulation) via Tensors.jl.

The function computes

f_int = ∫_Ω B_Pᵀ · P dΩ

where P is interpolated to Gauss points and contracted with the gradients of the shape functions. The routine is stress-measure agnostic: it treats P purely as a second-order tensor field and does not enforce any specific constitutive interpretation.

Arguments

• problem::Problem Finite element problem definition. The routine currently supports only:

problem.dim  == 3
problem.pdim == 3

• P::TensorField

Nodal stress tensor field. Each node stores a full dim×dim tensor, which is interpolated to Gauss points using Lagrange shape functions.

Typical choices for P include:

• First Piola–Kirchhoff stress,
• Second Piola–Kirchhoff stress (with compatible kinematics),
• Cauchy stress (if used consistently).

The function does not distinguish between these cases; ensuring consistency with the rest of the formulation is the responsibility of the caller.

• F::TensorField Nodal deformation gradient field. Required when using the energy-based path to compute P internally.

• energy::Function Free-energy density ψ(C, params) of a hyperelastic material. If provided, the routine constructs C = FᵀF, computes S via Tensors.jl, and forms P = F * S at Gauss points.

• params Optional parameter container passed through to energy.

Mathematical formulation

At each Gauss point, the stress tensor is obtained as

P_gp = Σ_a N_a P_a

The element internal force contribution is computed as

(f_e)*{ai} = ∫ (P_gp)*{iJ} (∇N_a)_J dΩ

where:

• a denotes the node index,
• i denotes the displacement component,
• ∇N_a is the gradient of the shape function in the reference configuration.

The resulting element vector is assembled into the global internal force vector.

Notes

• This function assembles only the internal force vector. The associated tangent contributions must be assembled separately via:

  • materialTangentMatrix (constitutive tangent),
  • initialStressMatrix (geometric tangent).

• The stress field P is interpolated using standard Lagrange shape functions. No symmetry or push-forward/pull-back operations are applied internally.

• If energy is supplied, P is computed at Gauss points from the free-energy function using Tensors.jl, enabling arbitrary hyperelastic laws without explicitly providing a stress field.

• The formulation corresponds to a Total Lagrange–type internal force expression when P is interpreted as the first Piola–Kirchhoff stress.

Returns

• VectorField Global internal force vector of size (ndof), returned as a VectorField with one time step.

Typical usage

f_int = internalForceVector(problem, Pfield)

R = f_int - f_ext

See also

• materialTangentMatrix
• initialStressMatrix
• externalTangentFollower
• TensorField

materialTangentMatrix(
    problem::Problem;
    F::TensorField,
    C::AbstractMatrix,
    energy::Function,
    params
)

Assembles the material (constitutive) tangent stiffness matrix for a 3D solid under finite deformation (Total Lagrange formulation), using a user-provided deformation gradient field F and a material tangent matrix C given in 6×6 Mandel/Voigt form. Alternatively, it can compute the material tangent from a free-energy density energy (hyperelastic formulation) via Tensors.jl.

This function builds the matrix

K_mat = ∫_Ω B(F)ᵀ · C · B(F) dΩ

where B(F) is the nonlinear strain–displacement matrix associated with the Green–Lagrange strain, and C is the material tangent ∂S/∂E expressed in Mandel notation.

The routine is material-agnostic: it does not assume any specific constitutive law. The material behavior is entirely defined by the supplied C, which may be spatially varying.

Arguments

• problem::Problem The finite element problem definition. Must satisfy:
  • problem.dim == 3
  • problem.pdim == 3

Keyword arguments

• F::TensorField Nodal deformation gradient field. Each node stores the full 3×3 deformation gradient F, which is interpolated to Gauss points during integration.

• C::AbstractMatrix (size 6×6) Material tangent matrix in Mandel/Voigt notation.

  Each entry C[i,j] may be either:

  • a Number (constant material tangent), or
  • a ScalarField (spatially varying material tangent component).

  All entries must satisfy:

C[i,j] isa Number || C[i,j] isa ScalarField

• energy::Function Free-energy density ψ(C, params) of a hyperelastic material, evaluated at the right Cauchy–Green tensor C = FᵀF. If provided, the tangent is computed at Gauss points via Tensors.jl and C must be omitted.

• params Optional parameter container passed through to energy.

Mathematical formulation

• Strain measure: Green–Lagrange strain

E = 1/2 (FᵀF − I)

• Variation:

δE = sym(Fᵀ δF)

• Tangent contribution:

δS = C : δE

• Element stiffness contribution:

K_e = ∫ B(F)ᵀ · C · B(F) dΩ

The strain–displacement matrix B(F) is constructed consistently with Mandel notation (shear components scaled by 2).

Notes

• This function assembles only the material part of the tangent

stiffness matrix. The geometric stiffness (stress-dependent part) must be assembled separately.

• The formulation is suitable for:
• hyperelastic materials,
• finite rotations,
• finite strains,

provided that C is consistent with the stress measure used elsewhere.

• If energy is supplied, the routine constructs C = FᵀF at Gauss points and obtains the constitutive tangent from the free-energy function using Tensors.jl. This enables arbitrary hyperelastic laws without providing a closed-form C.

• The material tangent C is evaluated at Gauss points by interpolating

nodal ScalarField values when necessary.

Returns

• SystemMatrix

Sparse global material tangent stiffness matrix of size (ndof × ndof).

Typical usage

Kmat = materialTangentMatrix(
  problem;
  F = Ffield,
  C = Cvoigt   # 6×6 matrix of Number / ScalarField
)

K = Kint + Kmat

See also

• initialStressMatrix

• internalForceVector

• externalTangentFollower

• TensorField

• ScalarField

nodePositionVector(problem)

Returns the position vectors of all mesh nodes as a VectorField (initial configuration). Returning vector is always a 3D vector.

Returns: R

Types:

• problem: Problem
• R: VectorField

rot(r::VectorField)

Solves the rotation of the vector field r. In some countries "rot" denotes the English "curl". (See the curl function.)

Return: VectorField

Types:

• r: VectorField

showDeformationResults(r::VectorField, comp; name=String, visible=Boolean)

Shows deformation result, where r contains the position vectors of nodes in the current configuration.

∇(r::Union{VectorField, ScalarField, TensorField}; nabla=:grad)

Computes derivatives of r.

• If r is a ScalarField and nabla == :grad, returns the gradient (a VectorField).
• If r is a VectorField and nabla == :grad, returns the gradient (a TensorField).
• If r is a VectorField and nabla == :curl, returns the curl (a VectorField).
• If r is a VectorField and nabla == :div, returns the divergence (a ScalarField).
• If r is a TensorField and nabla == :div, returns the divergence (a VectorField).

Returns: ScalarField, VectorField, or TensorField

Types:

• r: ScalarField, VectorField, or TensorField
• nabla: Symbol

3D Examples (assumes problem is set as in the ∇ doc setup)

# One-time 3D setup (assumes examples/Fields/cube.geo exists with physical group "body")
using LowLevelFEM
gmsh.initialize()
gmsh.open("examples/Fields/cube.geo")
mat = material("body", E=210e3, ν=0.3, ρ=7.85e-9)
problem = Problem([mat], type=:Solid)

# 1) Gradient of a 3D scalar field: ∇f → VectorField
f(X,Y,Z) = X^2 + Y*Z
S = scalarField(problem, [field("body", f=f)])
G = ∇(S)  # VectorField with 3 components

# 2) Curl of a 3D vector field: ∇ × v → VectorField
vx(X,Y,Z) = 0
vy(X,Y,Z) = X
vz(X,Y,Z) = 0
V = vectorField(problem, [field("body", fx=vx, fy=vy, fz=vz)])
C = ∇(V, nabla=:curl)  # approx (0, 0, 1) everywhere

# 3) Divergence of a 3D vector field: ∇ ⋅ v → ScalarField
v1(X,Y,Z) = X
v2(X,Y,Z) = Y
v3(X,Y,Z) = Z
V2 = vectorField(problem, [field("body", fx=v1, fy=v2, fz=v3)])
D = ∇(V2, nabla=:div)  # ≈ 3

# 4) Divergence of a 3D tensor field: ∇ · T → VectorField (if T is TensorField)
# For example, a diagonal tensor T with only Tzz = g(Z): div(T) = (0, 0, ∂g/∂Z)
g(Z) = 10 - Z
T = tensorField(problem, [field("body", fz=g)])
DV = ∇(T, nabla=:div)  # VectorField

# Symmetric displacement gradient via operators
# A = (u ∘ ∇ + ∇ ∘ u) / 2
gmsh.finalize()