Solid Mechanics Example¶

# dependencies
using LFAToolkit
using LinearAlgebra
using Pkg
Pkg.activate("./")
Pkg.instantiate()
using Plots

# setup
mesh = Mesh3D(1.0, 1.0, 1.0)
finep = 3
coarsep = 2
numbercomponents = 1
dimension = 3
finebasis = TensorH1LagrangeBasis(finep, finep, numbercomponents, dimension)
coarsebasis = TensorH1LagrangeBasis(coarsep, finep, numbercomponents, dimension)
ctofbasis = TensorH1LagrangeBasis(coarsep, finep, numbercomponents, dimension, lagrangequadrature = 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

    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, ctofbasis, ctofbasis])

# full operator symbols
numberruns = 250
maxeigenvalue = 0
θ_min = -π/2
θ_max = 3π/2

# compute and plot smoothing factor
# -- 1D --
if dimension == 1
    # setup
    ω = [0.636]
    v = [1, 1]
    maxeigenvalues = zeros(numberruns)

    # compute
    for i in 1:numberruns
        θ = [θ_min + (θ_max - θ_min)*i/numberruns]
        if abs(θ[1]  % 2π) >  π/128
            A = computesymbols(multigrid, ω, v, θ)
            eigenvalues = [abs(val) for val in eigvals(A)]
            maxeigenvalues[i] = max(eigenvalues...)
            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])
        end
    end

    # plot
    println("max eigenvalue: ", maxeigenvalue)
    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π
    plot(
        xrange,
        xlabel="θ/π",
        xtickfont=font(12, "Courier"),
        maxeigenvalues,
        ytickfont=font(12, "Courier"),
        ylabel="spectral radius",
        linewidth=3,
        legend=:none,
        title="P-Multigrid Error Symbol"
    )
    ylims!(0.0, max(maxeigenvalues...) * 1.1)
# -- 2D --
elseif dimension == 2
    # setup
    ω = [0.636]
    v = [1, 1]
    maxeigenvalues = zeros(numberruns, numberruns)

    # compute
    for i in 1:numberruns, j in 1:numberruns
        θ = [
            θ_min + (θ_max - θ_min)*i/numberruns,
            θ_min + (θ_max - θ_min)*j/numberruns
        ]
        if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π/128
            A = computesymbols(multigrid, ω, v, θ)
            eigenvalues = [abs(val) for val in eigvals(A)]
            maxeigenvalues[i, j] = max(eigenvalues...)
            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])
        end
    end

    # plot
    println("max eigenvalue: ", maxeigenvalue)
    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π
    heatmap(
        xrange,
        xlabel="θ/π",
        xtickfont=font(12, "Courier"),
        xrange,
        ylabel="θ/π",
        ytickfont=font(12, "Courier"),
        maxeigenvalues,
        title="P-Multigrid Error Symbol",
        transpose=true,
        aspect_ratio=:equal
    )
    xlims!(θ_min/π, θ_max/π)
    ylims!(θ_min/π, θ_max/π)
end