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 = 5)
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)
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())
pcoarse_solver = SmoothedAggregationCoarseSolver(; B)
builder = PMultigridPreconBuilder(DofHandlerHierarchy([dh_coarse, dh]), 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" _, ch_gal = IterativeSolvers.cg!(ΔΔu, K, g; Pl, maxiter = 100, log=true, verbose=false)
@info "Galerkin CG iterations: $(ch_gal.iters)"
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: 3993
[ Info: Galerkin CG iterations: 13
-------------------------------------------------------------------------------------------
Analysis with 750 elements Time Allocations
----------------------- ------------------------
Tot / % measured: 603ms / 54.0% 52.0MiB / 63.6%
Section ncalls time %tot avg alloc %tot avg
-------------------------------------------------------------------------------------------
Setup preconditioner 1 125ms 38.3% 125ms 28.3MiB 85.6% 28.3MiB
pmultigrid hierarchy 1 52.3ms 16.1% 52.3ms 8.34MiB 25.2% 8.34MiB
RAP numeric 1 35.3ms 10.9% 35.3ms 752B 0.0% 752B
coarse solver setup 1 17.0ms 5.2% 17.0ms 8.29MiB 25.1% 8.29MiB
extend_hierarchy! 2 16.8ms 5.2% 8.41ms 8.24MiB 24.9% 4.12MiB
improve candidates 2 8.19ms 2.5% 4.09ms 0.00B 0.0% 0.00B
RAP 2 4.59ms 1.4% 2.29ms 1.33MiB 4.0% 682KiB
restriction setup 2 2.57ms 0.8% 1.29ms 4.09MiB 12.4% 2.05MiB
strength 2 869μs 0.3% 434μs 2.39MiB 7.2% 1.19MiB
fit candidates 2 567μs 0.2% 284μs 384KiB 1.1% 192KiB
aggregation 2 28.6μs 0.0% 14.3μs 18.2KiB 0.1% 9.11KiB
smoother setup 2 942ns 0.0% 471ns 704B 0.0% 352B
coarse solver setup 1 112μs 0.0% 112μs 13.0KiB 0.0% 13.0KiB
prologue 1 7.74μs 0.0% 7.74μs 35.8KiB 0.1% 35.8KiB
smoother setup 1 831ns 0.0% 831ns 576B 0.0% 576B
build prolongator 1 36.5ms 11.2% 36.5ms 13.5MiB 40.9% 13.5MiB
setup transfer operator 1 33.7ms 10.4% 33.7ms 13.5MiB 40.8% 13.5MiB
assemble transfer operator 1 2.59ms 0.8% 2.59ms 1.44KiB 0.0% 1.44KiB
row normalization 1 144μs 0.0% 144μs 0.00B 0.0% 0.00B
RAP symbolic 1 35.7ms 11.0% 35.7ms 6.43MiB 19.4% 6.43MiB
build restriction 1 611ns 0.0% 611ns 48.0B 0.0% 48.0B
GMRES 1 96.0ms 29.5% 96.0ms 696KiB 2.1% 696KiB
CG 1 36.7ms 11.3% 36.7ms 94.6KiB 0.3% 94.6KiB
Galerkin GMRES 1 33.0ms 10.1% 33.0ms 1.57MiB 4.7% 1.57MiB
Coarse solve 13 10.4ms 3.2% 800μs 95.3KiB 0.3% 7.33KiB
Presmoother 26 4.39ms 1.3% 169μs 0.00B 0.0% 0.00B
Postsmoother 26 4.35ms 1.3% 167μs 0.00B 0.0% 0.00B
Residual eval 26 1.30ms 0.4% 49.8μs 0.00B 0.0% 0.00B
Restriction 26 156μs 0.0% 6.01μs 0.00B 0.0% 0.00B
Prolongation 26 121μs 0.0% 4.65μs 0.00B 0.0% 0.00B
Coarse solve 13 18.1μs 0.0% 1.39μs 14.0KiB 0.0% 1.08KiB
Presmoother 13 7.65ms 2.4% 588μs 624B 0.0% 48.0B
Postsmoother 13 7.59ms 2.3% 584μs 624B 0.0% 48.0B
Residual eval 13 2.49ms 0.8% 192μs 0.00B 0.0% 0.00B
Restriction 13 874μs 0.3% 67.2μs 0.00B 0.0% 0.00B
Prolongation 13 739μs 0.2% 56.9μs 0.00B 0.0% 0.00B
Galerkin CG 1 32.8ms 10.1% 32.8ms 603KiB 1.8% 603KiB
Coarse solve 13 10.3ms 3.2% 792μs 95.3KiB 0.3% 7.33KiB
Presmoother 26 4.32ms 1.3% 166μs 0.00B 0.0% 0.00B
Postsmoother 26 4.26ms 1.3% 164μs 0.00B 0.0% 0.00B
Residual eval 26 1.31ms 0.4% 50.3μs 0.00B 0.0% 0.00B
Restriction 26 184μs 0.1% 7.07μs 0.00B 0.0% 0.00B
Prolongation 26 121μs 0.0% 4.65μs 0.00B 0.0% 0.00B
Coarse solve 13 20.8μs 0.0% 1.60μs 14.0KiB 0.0% 1.08KiB
Postsmoother 13 7.70ms 2.4% 592μs 0.00B 0.0% 0.00B
Presmoother 13 7.59ms 2.3% 584μs 0.00B 0.0% 0.00B
Residual eval 13 2.53ms 0.8% 194μs 0.00B 0.0% 0.00B
Restriction 13 905μs 0.3% 69.6μs 0.00B 0.0% 0.00B
Prolongation 13 744μs 0.2% 57.2μs 0.00B 0.0% 0.00B
export 1 2.36ms 0.7% 2.36ms 1.82MiB 5.5% 1.82MiB
-------------------------------------------------------------------------------------------This page was generated using Literate.jl.