Nonlinear Elasticity
The full explanation for the underlying FEM theory in this example can be found in the Hyperelasticity tutorial of the Ferrite.jl documentation.
Implementation
The following code is based on the Hyperelasticity tutorial from the Ferrite.jl documentation, with some comments removed for brevity. There are two main modifications:
- Second-order
Lagrangeshape functions are used for field approximation:ip = Lagrange{RefTriangle,2}()^2. - Four quadrature points are used to accommodate the second-order shape functions:
qr = QuadratureRule{RefTriangle}(4).
using Ferrite, Tensors, TimerOutputs, IterativeSolvers
using FerriteMultigrid
TimerOutputs.enable_debug_timings(AlgebraicMultigrid)
TimerOutputs.enable_debug_timings(FerriteMultigrid)
struct NeoHooke
μ::Float64
λ::Float64
end
function Ψ(C, mp::NeoHooke)
μ = mp.μ
λ = mp.λ
Ic = tr(C)
J = sqrt(det(C))
return μ / 2 * (Ic - 3 - 2 * log(J)) + λ / 2 * (J - 1)^2
end
function constitutive_driver(C, mp::NeoHooke)
# Compute all derivatives in one function call
∂²Ψ∂C², ∂Ψ∂C = Tensors.hessian(y -> Ψ(y, mp), C, :all)
S = 2.0 * ∂Ψ∂C
∂S∂C = 2.0 * ∂²Ψ∂C²
return S, ∂S∂C
end;
function assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN)
# Reinitialize cell values, and reset output arrays
reinit!(cv, cell)
fill!(ke, 0.0)
fill!(ge, 0.0)
b = Vec{3}((0.0, -0.5, 0.0)) # Body force
tn = 0.1 # Traction (to be scaled with surface normal)
ndofs = getnbasefunctions(cv)
for qp in 1:getnquadpoints(cv)
dΩ = getdetJdV(cv, qp)
# Compute deformation gradient F and right Cauchy-Green tensor C
∇u = function_gradient(cv, qp, ue)
F = one(∇u) + ∇u
C = tdot(F) # F' ⋅ F
# Compute stress and tangent
S, ∂S∂C = constitutive_driver(C, mp)
P = F ⋅ S
I = one(S)
∂P∂F = otimesu(I, S) + 2 * F ⋅ ∂S∂C ⊡ otimesu(F', I)
# Loop over test functions
for i in 1:ndofs
# Test function and gradient
δui = shape_value(cv, qp, i)
∇δui = shape_gradient(cv, qp, i)
# Add contribution to the residual from this test function
ge[i] += (∇δui ⊡ P - δui ⋅ b) * dΩ
∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation
for j in 1:ndofs
∇δuj = shape_gradient(cv, qp, j)
# Add contribution to the tangent
ke[i, j] += (∇δui∂P∂F ⊡ ∇δuj) * dΩ
end
end
end
# Surface integral for the traction
for facet in 1:nfacets(cell)
if (cellid(cell), facet) in ΓN
reinit!(fv, cell, facet)
for q_point in 1:getnquadpoints(fv)
t = tn * getnormal(fv, q_point)
dΓ = getdetJdV(fv, q_point)
for i in 1:ndofs
δui = shape_value(fv, q_point, i)
ge[i] -= (δui ⋅ t) * dΓ
end
end
end
end
return
end;
function assemble_global!(K, g, dh, cv, fv, mp, u, ΓN)
n = ndofs_per_cell(dh)
ke = zeros(n, n)
ge = zeros(n)
# start_assemble resets K and g
assembler = start_assemble(K, g)
# Loop over all cells in the grid
for cell in CellIterator(dh)
global_dofs = celldofs(cell)
ue = u[global_dofs] # element dofs
assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN)
assemble!(assembler, global_dofs, ke, ge)
end
return
end;Near Null Space (NNS)
In multigrid methods for problems with vector-valued unknowns, such as elasticity, the near null space represents the low energy mode or the smooth error that needs to be captured in the coarser grid when using SA-AMG (Smoothed Aggregation Algebraic Multigrid), more on the topic can be found in Schroder [1].
For 3D linear elasticity problems, the rigid body modes are:
- Translation in the x-direction,
- Translation in the y-direction,
- Translation in the z-direction,
- Rotation about the x-axis (i.e., $x_1$): each point (x, y, z) is mapped to (0, -z, y).
- Rotation about the y-axis (i.e., $x_2$): each point (x, y, z) is mapped to (z, 0, -x).
- Rotation about the z-axis (i.e., $x_3$): each point (x, y, z) is mapped to (-y, x, 0).
The function create_nns constructs the NNS matrix B ∈ ℝ^{n × 6}, where n is the number of degrees of freedom (DOFs) for the case of p = 1 (i.e., linear interpolation), because B is only relevant for AMG.
function create_nns(dh, fieldname = first(dh.field_names))
@assert length(dh.field_names) == 1 "Only a single field is supported for now."
coords_flat = zeros(ndofs(dh))
apply_analytical!(coords_flat, dh, fieldname, x -> x)
coords = reshape(coords_flat, (length(coords_flat) ÷ 3, 3))
grid = dh.grid
B = zeros(Float64, ndofs(dh), 6)
B[1:3:end, 1] .= 1 # x - translation
B[2:3:end, 2] .= 1 # y - translation
B[3:3:end, 3] .= 1 # z - translation
# rotations
x = coords[:, 1]
y = coords[:, 2]
z = coords[:, 3]
# Around x
B[2:3:end, 4] .= -z
B[3:3:end, 4] .= y
# Around y
B[1:3:end, 5] .= z
B[3:3:end, 5] .= -x
# Around z
B[1:3:end, 6] .= -y
B[2:3:end, 6] .= x
return B
end
function _solve(N = 10)
reset_timer!()
# Generate a grid
L = 1.0
left = zero(Vec{3})
right = L * ones(Vec{3})
grid = generate_grid(Tetrahedron, (N, N, N), left, right)
# Material parameters
E = 10.0
ν = 0.3
μ = E / (2(1 + ν))
λ = (E * ν) / ((1 + ν) * (1 - 2ν))
mp = NeoHooke(μ, λ)
# Finite element base
ip = Lagrange{RefTetrahedron, 2}()^3
qr = QuadratureRule{RefTetrahedron}(4)
qr_facet = FacetQuadratureRule{RefTetrahedron}(3)
cv = CellValues(qr, ip)
fv = FacetValues(qr_facet, ip)
# DofHandler
dh = DofHandler(grid)
add!(dh, :u, ip) # Add a displacement field
close!(dh)
function rotation(X, t)
θ = pi / 3 # 60°
x, y, z = X
return t * Vec{3}(
(
0.0,
L / 2 - y + (y - L / 2) * cos(θ) - (z - L / 2) * sin(θ),
L / 2 - z + (y - L / 2) * sin(θ) + (z - L / 2) * cos(θ),
)
)
end
ch = ConstraintHandler(dh)
# Add a homogeneous boundary condition on the "clamped" edge
dbc = Dirichlet(:u, getfacetset(grid, "right"), (x, t) -> [0.0, 0.0, 0.0], [1, 2, 3])
add!(ch, dbc)
dbc = Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> rotation(x, t), [1, 2, 3])
add!(ch, dbc)
close!(ch)
t = 0.5
Ferrite.update!(ch, t)
# Neumann part of the boundary
ΓN = union(
getfacetset(grid, "top"),
getfacetset(grid, "bottom"),
getfacetset(grid, "front"),
getfacetset(grid, "back"),
)
# Pre-allocation of vectors for the solution and Newton increments
_ndofs = ndofs(dh)
un = zeros(_ndofs) # previous solution vector
u = zeros(_ndofs)
Δu = zeros(_ndofs)
ΔΔu = zeros(_ndofs)
apply!(un, ch)
# Create sparse matrix and residual vector
K = allocate_matrix(dh)
g = zeros(_ndofs)
# FIXME this needs better integration
dh_coarse = DofHandler(grid)
add!(dh_coarse, :u, Lagrange{RefTetrahedron, 1}()^3) # Add a displacement field
close!(dh_coarse)
B = create_nns(dh_coarse)
config_gal = pmultigrid_config(coarse_strategy = Galerkin())
fe_space = FESpace(dh, cv, ch)
pcoarse_solver = SmoothedAggregationCoarseSolver(; B)
builder = PMultigridPreconBuilder(fe_space, config_gal; pcoarse_solver)
# Perform Newton iterations
newton_itr = -1
NEWTON_TOL = 1.0e-8
NEWTON_MAXITER = 30
@info ndofs(dh)
while true
newton_itr += 1
# Construct the current guess
u .= un .+ Δu
# Compute residual and tangent for current guess
assemble_global!(K, g, dh, cv, fv, mp, u, ΓN)
# Apply boundary conditions
apply_zero!(K, g, ch)
# Compute the residual norm and compare with tolerance
normg = norm(g)
if normg < NEWTON_TOL
break
elseif newton_itr > NEWTON_MAXITER
error("Reached maximum Newton iterations, aborting")
end
# Compute increment using conjugate gradients
fill!(ΔΔu, 0.0)
@timeit "Setup preconditioner" Pl = builder(K)[1]
@timeit "Galerkin CG" IterativeSolvers.cg!(ΔΔu, K, g; Pl, maxiter = 100, verbose=false)
fill!(ΔΔu, 0.0)
@timeit "Galerkin GMRES" IterativeSolvers.gmres!(ΔΔu, K, g; Pl, maxiter = 100, verbose=false)
fill!(ΔΔu, 0.0)
@timeit "CG" IterativeSolvers.cg!(ΔΔu, K, g; maxiter = 1000, verbose=false)
fill!(ΔΔu, 0.0)
@timeit "GMRES" IterativeSolvers.gmres!(ΔΔu, K, g; maxiter = 1000, verbose=false)
apply_zero!(ΔΔu, ch)
Δu .-= ΔΔu
break
end
# Save the solution
@timeit "export" begin
VTKGridFile("hyperelasticity", dh) do vtk
write_solution(vtk, dh, u)
end
end
print_timer(title = "Analysis with $(getncells(grid)) elements", linechars = :ascii)
return u
end
u = _solve();[ Info: 27783
---------------------------------------------------------------------------------------
Analysis with 6000 elements Time Allocations
----------------------- ------------------------
Tot / % measured: 4.97s / 86.4% 441MiB / 78.8%
Section ncalls time %tot avg alloc %tot avg
---------------------------------------------------------------------------------------
GMRES 1 1.78s 41.5% 1.78s 4.68MiB 1.3% 4.68MiB
CG 1 893ms 20.8% 893ms 652KiB 0.2% 652KiB
Setup preconditioner 1 858ms 20.0% 858ms 319MiB 91.7% 319MiB
pmultigrid hierarchy 1 858ms 20.0% 858ms 319MiB 91.7% 319MiB
extend_hierarchy! 1 817ms 19.0% 817ms 238MiB 68.4% 238MiB
build prolongator 1 796ms 18.5% 796ms 150MiB 43.2% 150MiB
assembly 1 796ms 18.5% 796ms 150MiB 43.1% 150MiB
row normalization 1 82.7μs 0.0% 82.7μs 0.00B 0.0% 0.00B
RAP 1 21.5ms 0.5% 21.5ms 87.7MiB 25.2% 87.7MiB
build restriction 1 410ns 0.0% 410ns 48.0B 0.0% 48.0B
coarse solver setup 1 39.4ms 0.9% 39.4ms 78.1MiB 22.5% 78.1MiB
extend_hierarchy! 3 39.4ms 0.9% 13.1ms 77.9MiB 22.4% 26.0MiB
RAP 3 12.1ms 0.3% 4.04ms 33.7MiB 9.7% 11.2MiB
restriction setup 3 11.8ms 0.3% 3.93ms 38.1MiB 10.9% 12.7MiB
improve candidates 3 10.3ms 0.2% 3.43ms 0.00B 0.0% 0.00B
fit candidates 3 4.12ms 0.1% 1.37ms 2.79MiB 0.8% 952KiB
strength 3 950μs 0.0% 317μs 2.99MiB 0.9% 1.00MiB
aggregation 3 80.4μs 0.0% 26.8μs 144KiB 0.0% 47.9KiB
coarse solver setup 1 48.2μs 0.0% 48.2μs 6.57KiB 0.0% 6.57KiB
prologue 1 13.2μs 0.0% 13.2μs 219KiB 0.1% 219KiB
ml setup 1 4.68μs 0.0% 4.68μs 224B 0.0% 224B
coarsen order 1 1.11ms 0.0% 1.11ms 2.49MiB 0.7% 2.49MiB
Galerkin GMRES 1 384ms 8.9% 384ms 13.4MiB 3.9% 13.4MiB
Postsmoother 19 146ms 3.4% 7.70ms 0.00B 0.0% 0.00B
Presmoother 19 146ms 3.4% 7.68ms 0.00B 0.0% 0.00B
Residual eval 19 28.4ms 0.7% 1.49ms 0.00B 0.0% 0.00B
Coarse solve 19 23.5ms 0.5% 1.24ms 684KiB 0.2% 36.0KiB
Presmoother 57 8.94ms 0.2% 157μs 0.00B 0.0% 0.00B
Postsmoother 57 8.20ms 0.2% 144μs 0.00B 0.0% 0.00B
Residual eval 57 2.92ms 0.1% 51.3μs 0.00B 0.0% 0.00B
Restriction 57 1.57ms 0.0% 27.6μs 0.00B 0.0% 0.00B
Prolongation 57 1.31ms 0.0% 23.1μs 0.00B 0.0% 0.00B
Coarse solve 19 29.7μs 0.0% 1.56μs 0.00B 0.0% 0.00B
Prolongation 19 1.07ms 0.0% 56.4μs 0.00B 0.0% 0.00B
Restriction 19 894μs 0.0% 47.1μs 0.00B 0.0% 0.00B
Galerkin CG 1 362ms 8.4% 362ms 5.34MiB 1.5% 5.34MiB
Postsmoother 19 141ms 3.3% 7.40ms 0.00B 0.0% 0.00B
Presmoother 19 136ms 3.2% 7.15ms 0.00B 0.0% 0.00B
Residual eval 19 27.8ms 0.6% 1.46ms 0.00B 0.0% 0.00B
Coarse solve 19 23.2ms 0.5% 1.22ms 684KiB 0.2% 36.0KiB
Presmoother 57 8.78ms 0.2% 154μs 0.00B 0.0% 0.00B
Postsmoother 57 8.16ms 0.2% 143μs 0.00B 0.0% 0.00B
Residual eval 57 2.85ms 0.1% 50.0μs 0.00B 0.0% 0.00B
Restriction 57 1.56ms 0.0% 27.4μs 0.00B 0.0% 0.00B
Prolongation 57 1.27ms 0.0% 22.2μs 0.00B 0.0% 0.00B
Coarse solve 19 32.5μs 0.0% 1.71μs 0.00B 0.0% 0.00B
Prolongation 19 993μs 0.0% 52.2μs 0.00B 0.0% 0.00B
Restriction 19 853μs 0.0% 44.9μs 0.00B 0.0% 0.00B
export 1 14.3ms 0.3% 14.3ms 4.88MiB 1.4% 4.88MiB
---------------------------------------------------------------------------------------This page was generated using Literate.jl.