Linear elasticity
This is an example of a p-multigrid preconditioner with Chebyshev smoothing for the 3D linear elasticity problem.
Problem formulation
The strong form of the static balance of linear-momentum at finite strain is given by
\[\nabla \cdot \boldsymbol{\sigma} + \bold{g} = \bold{0}\]
where $\boldsymbol{\sigma}$ and $\bold{g}$ are the stress and forcing functions, respectively.
We have the weak form
\[\int_{\Omega} \nabla \bold{v} : \boldsymbol{\sigma} - \int_{\partial \Omega} v \cdot \left( \boldsymbol{\sigma} \cdot \hat{\bold{n}} \right) - \int_{\Omega} \bold{v} \cdot \bold{g} = 0.\]
The constitutive law (stress-strain relationship) can be written as
\[\boldsymbol{\sigma} = \bold{C} : \boldsymbol{\epsilon}\]
where $\boldsymbol{\epsilon}$ is the infinitesimal strain tensor
\[\boldsymbol{\epsilon} = \frac{1}{2} \left( \nabla \bold{u} + \nabla \bold{u}^T \right)\]
and the elasticity tensor $\bold{C}$ is given by
\[\bold{C} = \begin{pmatrix} \lambda + 2\mu & \lambda & \lambda & & & \\ \lambda & \lambda + 2\mu & \lambda & & & \\ \lambda & \lambda & \lambda + 2\mu & & & \\ & & & \mu & & \\ & & & & \mu & \\ & & & & & \mu \end{pmatrix}.\]
$\lambda$ and $\mu$ are the Lamé parameters, given by
\[\lambda = \frac{E \nu}{\left( 1 + \nu \right) \left( 1 - 2 \nu \right)} \\ \mu = \frac{E}{2 \left( 1 + \nu \right)}\]
where $E$ is the Young's modulus and $\nu$ is the Poisson's ratio for the materiel.
For a full discussion of the formulation of the linear elasticity problem, see the libCEED documentation.
LFAToolkit code
# ------------------------------------------------------------------------------
# linear elasticity 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
λ = E * ν / ((1 + ν) * (1 - 2 * ν)) # Lamé parameters
μ = E / (2 * (1 + ν))
function linearelasticityweakform(deltadu::Array{Float64}, w::Array{Float64})
# strain
dϵ = (deltadu + deltadu') / 2
# strain energy
dσ_11 = (λ + 2μ) * dϵ[1, 1] + λ * dϵ[2, 2] + λ * dϵ[3, 3]
dσ_22 = λ * dϵ[1, 1] + (λ + 2μ) * dϵ[2, 2] + λ * dϵ[3, 3]
dσ_33 = λ * dϵ[1, 1] + λ * dϵ[2, 2] + (λ + 2μ) * dϵ[3, 3]
dσ_12 = μ * dϵ[1, 2]
dσ_13 = μ * dϵ[1, 3]
dσ_23 = μ * dϵ[2, 3]
dσ = [dσ_11 dσ_12 dσ_13; dσ_12 dσ_22 dσ_23; dσ_13 dσ_23 dσ_33]
# delta dv
deltadv = dσ * w[1]
return [deltadv']
end
# linear elasticity 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(linearelasticityweakform, 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))
# ------------------------------------------------------------------------------