Neo-Hookean hyperelasticity
This is an example of a p-multigrid preconditioner with Chebyshev smoothing for the 3D Neo-Hookean hyperelasticity problem.
Problem formulation
The strong form of the static balance of linear-momentum at finite strain is given by
\[-\nabla_x \cdot \bold{P} - \rho_0 \bold{g} = \bold{0}\]
where $-\nabla_k$ is the gradient with respect to the reference configuration, $\bold{P}$ is the first Piola-Kirchhoff stress tensor, $\rho_0$ is the reference mass density, and $\bold{g}$ is the forcing function.
The first Piola-Kirchhoff stress tensor is given by
\[\bold{P} = \bold{F} \bold{S}\]
where $\bold{F}$ is the deformation gradient and $\bold{S}$ is the second Piola-Kirchhoff stress tensor. In this example, the second Piola-Kirchhoff stress tensor is given by the Neo-Hookean model.
For a full discussion of the formulation of the 3D Neo-Hookean hyperelasticity problem, see the libCEED documentation.
LFAToolkit code
# ------------------------------------------------------------------------------
# hyperelasticity example
# ------------------------------------------------------------------------------
using LFAToolkit
using LinearAlgebra
# setup
mesh = Mesh3D(1.0, 1.0, 1.0)
finep = 2
coarsep = 1
numbercomponents = 3
dimension = 3
finebasis = TensorH1LagrangeBasis(finep + 1, finep + 1, numbercomponents, dimension)
coarsebasis = TensorH1LagrangeBasis(coarsep + 1, finep + 1, numbercomponents, dimension)
ctofbasis = TensorH1LagrangeBasis(
coarsep + 1,
finep + 1,
numbercomponents,
dimension,
collocatedquadrature = true,
)
# constants
E = 1E6 # Young's modulus
ν = 0.3 # Poisson's ratio
K = E / (3 * (1 - 2 * ν)) # bulk modulus
λ = E * ν / ((1 + ν) * (1 - 2 * ν)) # Lamé parameters
μ = E / (2 * (1 + ν))
# state
gradu = [1; 2; 3] * ones(1, 3);
function neohookeanweakform(deltadu::Array{Float64}, w::Array{Float64})
# dP = dF S + F dS
# deformation gradient
F = gradu + I
J = det(F)
# Green-Lagrange strain tensor
E = (gradu * gradu' + gradu' * gradu) / 2
# right Cauchy-Green tensor
C = 2 * E + I
C_inv = C^-1
# second Piola-Kirchhoff
S = λ * log(J) * C_inv + 2 * μ * C_inv * E
# delta du
deltadu = deltadu'
# dF
dF = deltadu + I
# deltaE
deltaE = (deltadu * deltadu' + deltadu' * deltadu) / 2
# dS
dS = λ * sum(C_inv .* deltaE) * C_inv + 2 * (μ - λ * log(J)) * C_inv * deltaE * C_inv
# dP
dP = (dF * S + F * dS) * w[1]
return [dP']
end
# linearized Neo-Hookean operators
function makeoperator(basis::TensorBasis)
inputs = [
OperatorField(basis, [EvaluationMode.gradient], "gradent of deformation"),
OperatorField(basis, [EvaluationMode.quadratureweights], "quadrature weights"),
]
outputs = [
OperatorField(
basis,
[EvaluationMode.gradient],
"test function gradient of deformation",
),
]
return Operator(neohookeanweakform, mesh, inputs, outputs)
end
fineoperator = makeoperator(finebasis)
coarseoperator = makeoperator(coarsebasis)
# Chebyshev smoother
chebyshev = Chebyshev(fineoperator)
# p-multigrid preconditioner
multigrid = PMultigrid(fineoperator, coarseoperator, chebyshev, [ctofbasis])
# compute operator symbols
A = computesymbols(multigrid, [3], [1, 1], [π, π, π])
eigenvalues = real(eigvals(A))
# ------------------------------------------------------------------------------