Nearly Incompressible Hyperelasticity

Tip

This example is also available as a Jupyter notebook: quasi_incompressible_hyperelasticity.ipynb

Introduction

In this example we study quasi- or nearly-incompressible hyperelasticity using the stable Taylor-Hood approximation. In spirit, this example is the nonlinear analogue of incompressible_elasticity and the incompressible analogue of hyperelasticity. Much of the code therefore follows from the above two examples. The problem is formulated in the undeformed or reference configuration with the displacement $\mathbf{u}$ and pressure $p$ being the unknown fields. We now briefly outline the formulation. Consider the standard hyperelasticity problem

\[ \mathbf{u} = \argmin_{\mathbf{v}\in\mathcal{K}(\Omega)}\Pi(\mathbf{v}),\quad \text{where}\quad \Pi(\mathbf{v}) = \int_\Omega \Psi(\mathbf{v}) \ \mathrm{d}\Omega\ .\]

where $\mathcal{K}(\Omega)$ is a suitable function space.

For clarity of presentation we ignore any non-zero surface tractions and body forces and instead consider only applied displacements (i.e. non-homogeneous dirichlet boundary conditions). Moreover we stick our attention to the standard Neo-Hookean stored energy density

\[ \Psi(\mathbf{u}) = \frac{\mu}{2}\left(I_1 - 3 \right) - \mu \log(J) + \frac{\lambda}{2}\left( J - 1\right){}^2,\]

where $I_1 = \mathrm{tr}(\mathbf{C}) = \mathrm{tr}(\mathbf{F}^\mathrm{T} \mathbf{F}) = F_{ij}F_{ij}$ and $J = \det(\mathbf{F})$ denote the standard invariants of the deformation gradient tensor $\mathbf{F} = \mathbf{I}+\nabla_{\mathbf{X}} \mathbf{u}$. The above problem is ill-posed in the limit of incompressibility (or near-incompressibility), namely when

\[ \lambda/\mu \rightarrow +\infty.\]

In order to alleviate the problem, we consider the partial legendre transform of the strain energy density $\Psi$ with respect to $J = \det(\mathbf{F})$, namely

\[ \widehat{\Psi}(\mathbf{u}, p) = \sup_{J} \left[ p(J - 1) - \frac{\mu}{2}\left(I_1 - 3 \right) + \mu \log(J) - \frac{\lambda}{2}\left( J - 1\right){}^2 \right].\]

The supremum, say $J^\star$, can be calculated in closed form by the first order optimailty condition $\partial\widehat{\Psi}/\partial J = 0$. This gives

\[ J^\star(p) = \frac{\lambda + p + \sqrt{(\lambda + p){}^2 + 4 \lambda \mu }}{(2 \lambda)}.\]

Furthermore, taking the partial legendre transform of $\widehat{\Psi}$ once again, gives us back the original problem, i.e.

\[ \Psi(\mathbf{u}) = \Psi^\star(\mathbf{u}, p) = \sup_{p} \left[ p(J - 1) - p(J^\star - 1) + \frac{\mu}{2}\left(I_1 - 3 \right) - \mu \log(J^\star) + \frac{\lambda}{2}\left( J^\star - 1\right){}^2 \right].\]

Therefore our original hyperelasticity problem can now be reformulated as

\[ \inf_{\mathbf{u}\in\mathcal{K}(\Omega)}\sup_{p} \int_\Omega\Psi^{\star} (\mathbf{u}, p) \, \mathrm{d}\Omega.\]

The total (modified) energy $\Pi^\star$ can then be written as

\[ \Pi^\star(\mathbf{u}, p) = \int_\Omega p (J - J^\star) \ \mathrm{d}\Omega + \int_\Omega \frac{\mu}{2} \left( I_1 - 3\right) \ \mathrm{d}\Omega - \int_\Omega \mu\log(J^\star)\ \mathrm{d}\Omega + \int_\Omega \frac{\lambda}{2}\left( J^\star - 1 \right){}^2\ \mathrm{d}\Omega\]

The Euler-Lagrange equations corresponding to the above energy give us our governing PDEs in the weak form, namely

\[ \int_\Omega \frac{\partial\Psi^\star}{\partial \mathbf{F}}:\delta \mathbf{F} \ \mathrm{d}\Omega = 0\]

and

\[ \int_\Omega \frac{\partial \Psi^\star}{\partial p}\delta p \ \mathrm{d}\Omega = 0,\]

where $\delta \mathrm{F} = \delta \mathrm{grad}_0(\mathbf{u}) = \mathrm{grad}_0(\delta \mathbf{u})$ and $\delta \mathbf{u}$ and $\delta p$ denote arbitrary variations with respect to displacement and pressure (or the test functions). See the references below for a more detailed explanation of the above mathematical trick. Now, in order to apply Newton's method to the above problem, we further need to linearize the above equations and calculate the respective hessians (or tangents), namely, $\partial^2\Psi^\star/\partial \mathbf{F}^2$, $\partial^2\Psi^\star/\partial p^2$ and $\partial^2\Psi^\star/\partial \mathbf{F}\partial p$ which, using Tensors.jl, can be determined conveniently using automatic differentiation (see the code below). Hence we only need to define the above potential. The remaineder of the example follows similarly.

References

  1. A paradigm for higher-order polygonal elements in finite elasticity using a gradient correction scheme, CMAME 2016, 306, 216–251
  2. Approximation of incompressible large deformation elastic problems: some unresolved issues, Computational Mechanics, 2013

Implementation

We now get to the actual code. First, we import the respective packages

using Ferrite, Tensors, ProgressMeter, WriteVTK
using BlockArrays, SparseArrays, LinearAlgebra

and the corresponding struct to store our material properties.

struct NeoHooke
    μ::Float64
    λ::Float64
end

We then create a function to generate a simple test mesh on which to compute FE solution. We also mark the boundaries to later assign Dirichlet boundary conditions

function importTestGrid()
    grid = generate_grid(Tetrahedron, (5, 5, 5), zero(Vec{3}), ones(Vec{3}))
    addfacetset!(grid, "myBottom", x -> norm(x[2]) ≈ 0.0)
    addfacetset!(grid, "myBack", x -> norm(x[3]) ≈ 0.0)
    addfacetset!(grid, "myRight", x -> norm(x[1]) ≈ 1.0)
    addfacetset!(grid, "myLeft", x -> norm(x[1]) ≈ 0.0)
    return grid
end;

The function to create corresponding cellvalues for the displacement field u and pressure p follows in a similar fashion from the incompressible_elasticity example

function create_values(interpolation_u, interpolation_p)
    # quadrature rules
    qr = QuadratureRule{RefTetrahedron}(4)
    facet_qr = FacetQuadratureRule{RefTetrahedron}(4)

    # cell and facetvalues for u
    cellvalues_u = CellValues(qr, interpolation_u)
    facetvalues_u = FacetValues(facet_qr, interpolation_u)

    # cellvalues for p
    cellvalues_p = CellValues(qr, interpolation_p)

    return cellvalues_u, cellvalues_p, facetvalues_u
end;

We now create the function for Ψ*

function Ψ(F, p, mp::NeoHooke)
    μ = mp.μ
    λ = mp.λ
    Ic = tr(tdot(F))
    J = det(F)
    Js = (λ + p + sqrt((λ + p)^2.0 + 4.0 * λ * μ)) / (2.0 * λ)
    return p * (Js - J) + μ / 2 * (Ic - 3) - μ * log(Js) + λ / 2 * (Js - 1)^2
end;

and it's derivatives (required in computing the jacobian and hessian respectively)

function constitutive_driver(F, p, mp::NeoHooke)
    # Compute all derivatives in one function call
    ∂²Ψ∂F², ∂Ψ∂F = Tensors.hessian(y -> Ψ(y, p, mp), F, :all)
    ∂²Ψ∂p², ∂Ψ∂p = Tensors.hessian(y -> Ψ(F, y, mp), p, :all)
    ∂²Ψ∂F∂p = Tensors.gradient(q -> Tensors.gradient(y -> Ψ(y, q, mp), F), p)
    return ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p
end;

The functions to create the DofHandler and ConstraintHandler (to assign corresponding boundary conditions) follow likewise from the incompressible elasticity example, namely

function create_dofhandler(grid, ipu, ipp)
    dh = DofHandler(grid)
    add!(dh, :u, ipu) # displacement dim = 3
    add!(dh, :p, ipp) # pressure dim = 1
    close!(dh)
    return dh
end;

We are simulating a uniaxial tensile loading of a unit cube. Hence we apply a displacement field (:u) in x direction on the right face. The left, bottom and back facets are fixed in the x, y and z components of the displacement so as to emulate the uniaxial nature of the loading.

function create_bc(dh)
    dbc = ConstraintHandler(dh)
    add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myLeft"), (x, t) -> zero(Vec{1}), [1]))
    add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myBottom"), (x, t) -> zero(Vec{1}), [2]))
    add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myBack"), (x, t) -> zero(Vec{1}), [3]))
    add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myRight"), (x, t) -> t * ones(Vec{1}), [1]))
    close!(dbc)
    Ferrite.update!(dbc, 0.0)
    return dbc
end;

Also, since we are considering incompressible hyperelasticity, an interesting quantity that we can compute is the deformed volume of the solid. It is easy to show that this is equal to ∫J*dΩ where J=det(F). This can be done at the level of each element (cell)

function calculate_element_volume(cell, cellvalues_u, ue)
    reinit!(cellvalues_u, cell)
    evol::Float64 = 0.0
    for qp in 1:getnquadpoints(cellvalues_u)
        dΩ = getdetJdV(cellvalues_u, qp)
        ∇u = function_gradient(cellvalues_u, qp, ue)
        F = one(∇u) + ∇u
        J = det(F)
        evol += J * dΩ
    end
    return evol
end;

and then assembled over all the cells (elements)

function calculate_volume_deformed_mesh(w, dh::DofHandler, cellvalues_u)
    evol::Float64 = 0.0
    for cell in CellIterator(dh)
        global_dofs = celldofs(cell)
        nu = getnbasefunctions(cellvalues_u)
        global_dofs_u = global_dofs[1:nu]
        ue = w[global_dofs_u]
        δevol = calculate_element_volume(cell, cellvalues_u, ue)
        evol += δevol
    end
    return evol
end;

The function to assemble the element stiffness matrix for each element in the mesh now has a block structure like in incompressible_elasticity.

function assemble_element!(Ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
    # Reinitialize cell values, and reset output arrays
    ublock, pblock = 1, 2
    reinit!(cellvalues_u, cell)
    reinit!(cellvalues_p, cell)
    fill!(Ke, 0.0)
    fill!(fe, 0.0)

    n_basefuncs_u = getnbasefunctions(cellvalues_u)
    n_basefuncs_p = getnbasefunctions(cellvalues_p)

    for qp in 1:getnquadpoints(cellvalues_u)
        dΩ = getdetJdV(cellvalues_u, qp)
        # Compute deformation gradient F
        ∇u = function_gradient(cellvalues_u, qp, ue)
        p = function_value(cellvalues_p, qp, pe)
        F = one(∇u) + ∇u

        # Compute first Piola-Kirchhoff stress and tangent modulus
        ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p = constitutive_driver(F, p, mp)

        # Loop over the `u`-test functions to calculate the `u`-`u` and `u`-`p` blocks
        for i in 1:n_basefuncs_u
            # gradient of the test function
            ∇δui = shape_gradient(cellvalues_u, qp, i)
            # Add contribution to the residual from this test function
            fe[BlockIndex((ublock), (i))] += (∇δui ⊡ ∂Ψ∂F) * dΩ

            ∇δui∂S∂F = ∇δui ⊡ ∂²Ψ∂F²
            for j in 1:n_basefuncs_u
                ∇δuj = shape_gradient(cellvalues_u, qp, j)

                # Add contribution to the tangent
                Ke[BlockIndex((ublock, ublock), (i, j))] += (∇δui∂S∂F ⊡ ∇δuj) * dΩ
            end
            # Loop over the `p`-test functions
            for j in 1:n_basefuncs_p
                δp = shape_value(cellvalues_p, qp, j)
                # Add contribution to the tangent
                Ke[BlockIndex((ublock, pblock), (i, j))] += (∂²Ψ∂F∂p ⊡ ∇δui) * δp * dΩ
            end
        end
        # Loop over the `p`-test functions to calculate the `p-`u` and `p`-`p` blocks
        for i in 1:n_basefuncs_p
            δp = shape_value(cellvalues_p, qp, i)
            fe[BlockIndex((pblock), (i))] += (δp * ∂Ψ∂p) * dΩ

            for j in 1:n_basefuncs_u
                ∇δuj = shape_gradient(cellvalues_u, qp, j)
                Ke[BlockIndex((pblock, ublock), (i, j))] += ∇δuj ⊡ ∂²Ψ∂F∂p * δp * dΩ
            end
            for j in 1:n_basefuncs_p
                δp = shape_value(cellvalues_p, qp, j)
                Ke[BlockIndex((pblock, pblock), (i, j))] += δp * ∂²Ψ∂p² * δp * dΩ
            end
        end
    end
    return
end;

The only thing that changes in the assembly of the global stiffness matrix is slicing the corresponding element dofs for the displacement (see global_dofsu) and pressure (global_dofsp).

function assemble_global!(
        K::SparseMatrixCSC, f, cellvalues_u::CellValues,
        cellvalues_p::CellValues, dh::DofHandler, mp::NeoHooke, w
    )
    nu = getnbasefunctions(cellvalues_u)
    np = getnbasefunctions(cellvalues_p)

    # start_assemble resets K and f
    fe = BlockedArray(zeros(nu + np), [nu, np]) # local force vector
    ke = BlockedArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix

    assembler = start_assemble(K, f)
    # Loop over all cells in the grid
    for cell in CellIterator(dh)
        global_dofs = celldofs(cell)
        global_dofsu = global_dofs[1:nu] # first nu dofs are displacement
        global_dofsp = global_dofs[(nu + 1):end] # last np dofs are pressure
        @assert size(global_dofs, 1) == nu + np # sanity check
        ue = w[global_dofsu] # displacement dofs for the current cell
        pe = w[global_dofsp] # pressure dofs for the current cell
        assemble_element!(ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
        assemble!(assembler, global_dofs, ke, fe)
    end
    return
end;

We now define a main function solve. For nonlinear quasistatic problems we often like to parameterize the solution in terms of a pseudo time like parameter, which in this case is used to gradually apply the boundary displacement on the right face. Also for definitenessm we consider λ/μ = 10⁴

function solve(interpolation_u, interpolation_p)

    # import the mesh
    grid = importTestGrid()

    # Material parameters
    μ = 1.0
    λ = 1.0e4 * μ
    mp = NeoHooke(μ, λ)

    # Create the DofHandler and CellValues
    dh = create_dofhandler(grid, interpolation_u, interpolation_p)
    cellvalues_u, cellvalues_p, facetvalues_u = create_values(interpolation_u, interpolation_p)

    # Create the DirichletBCs
    dbc = create_bc(dh)

    # Pre-allocation of vectors for the solution and Newton increments
    _ndofs = ndofs(dh)
    w = zeros(_ndofs)
    ΔΔw = zeros(_ndofs)
    apply!(w, dbc)

    # Create the sparse matrix and residual vector
    K = allocate_matrix(dh)
    f = zeros(_ndofs)

    # We run the simulation parameterized by a time like parameter. `Tf` denotes the final value
    # of this parameter, and Δt denotes its increment in each step
    Tf = 2.0
    Δt = 0.1
    NEWTON_TOL = 1.0e-8

    pvd = paraview_collection("hyperelasticity_incomp_mixed")
    for (step, t) in enumerate(0.0:Δt:Tf)
        # Perform Newton iterations
        Ferrite.update!(dbc, t)
        apply!(w, dbc)
        newton_itr = -1
        prog = ProgressMeter.ProgressThresh(NEWTON_TOL; desc = "Solving @ time $t of $Tf;")
        fill!(ΔΔw, 0.0)
        while true
            newton_itr += 1
            assemble_global!(K, f, cellvalues_u, cellvalues_p, dh, mp, w)
            norm_res = norm(f[Ferrite.free_dofs(dbc)])
            apply_zero!(K, f, dbc)
            # Only display output at specific load steps
            if t % (5 * Δt) == 0
                ProgressMeter.update!(prog, norm_res; showvalues = [(:iter, newton_itr)])
            end
            if norm_res < NEWTON_TOL
                break
            elseif newton_itr > 30
                error("Reached maximum Newton iterations, aborting")
            end
            # Compute the incremental `dof`-vector (both displacement and pressure)
            ΔΔw .= K \ f

            apply_zero!(ΔΔw, dbc)
            w .-= ΔΔw
        end

        # Save the solution fields
        VTKGridFile("hyperelasticity_incomp_mixed_$step", grid) do vtk
            write_solution(vtk, dh, w)
            pvd[t] = vtk
        end
    end
    vtk_save(pvd)
    vol_def = calculate_volume_deformed_mesh(w, dh, cellvalues_u)
    print("Deformed volume is $vol_def")
    return vol_def
end;

We can now test the solution using the Taylor-Hood approximation

quadratic_u = Lagrange{RefTetrahedron, 2}()^3
linear_p = Lagrange{RefTetrahedron, 1}()
vol_def = solve(quadratic_u, linear_p)
1.000199977169989

The deformed volume is indeed close to 1 (as should be for a nearly incompressible material).

Plain program

Here follows a version of the program without any comments. The file is also available here: quasi_incompressible_hyperelasticity.jl.

using Ferrite, Tensors, ProgressMeter, WriteVTK
using BlockArrays, SparseArrays, LinearAlgebra

struct NeoHooke
    μ::Float64
    λ::Float64
end

function importTestGrid()
    grid = generate_grid(Tetrahedron, (5, 5, 5), zero(Vec{3}), ones(Vec{3}))
    addfacetset!(grid, "myBottom", x -> norm(x[2]) ≈ 0.0)
    addfacetset!(grid, "myBack", x -> norm(x[3]) ≈ 0.0)
    addfacetset!(grid, "myRight", x -> norm(x[1]) ≈ 1.0)
    addfacetset!(grid, "myLeft", x -> norm(x[1]) ≈ 0.0)
    return grid
end;

function create_values(interpolation_u, interpolation_p)
    # quadrature rules
    qr = QuadratureRule{RefTetrahedron}(4)
    facet_qr = FacetQuadratureRule{RefTetrahedron}(4)

    # cell and facetvalues for u
    cellvalues_u = CellValues(qr, interpolation_u)
    facetvalues_u = FacetValues(facet_qr, interpolation_u)

    # cellvalues for p
    cellvalues_p = CellValues(qr, interpolation_p)

    return cellvalues_u, cellvalues_p, facetvalues_u
end;

function Ψ(F, p, mp::NeoHooke)
    μ = mp.μ
    λ = mp.λ
    Ic = tr(tdot(F))
    J = det(F)
    Js = (λ + p + sqrt((λ + p)^2.0 + 4.0 * λ * μ)) / (2.0 * λ)
    return p * (Js - J) + μ / 2 * (Ic - 3) - μ * log(Js) + λ / 2 * (Js - 1)^2
end;

function constitutive_driver(F, p, mp::NeoHooke)
    # Compute all derivatives in one function call
    ∂²Ψ∂F², ∂Ψ∂F = Tensors.hessian(y -> Ψ(y, p, mp), F, :all)
    ∂²Ψ∂p², ∂Ψ∂p = Tensors.hessian(y -> Ψ(F, y, mp), p, :all)
    ∂²Ψ∂F∂p = Tensors.gradient(q -> Tensors.gradient(y -> Ψ(y, q, mp), F), p)
    return ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p
end;

function create_dofhandler(grid, ipu, ipp)
    dh = DofHandler(grid)
    add!(dh, :u, ipu) # displacement dim = 3
    add!(dh, :p, ipp) # pressure dim = 1
    close!(dh)
    return dh
end;

function create_bc(dh)
    dbc = ConstraintHandler(dh)
    add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myLeft"), (x, t) -> zero(Vec{1}), [1]))
    add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myBottom"), (x, t) -> zero(Vec{1}), [2]))
    add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myBack"), (x, t) -> zero(Vec{1}), [3]))
    add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "myRight"), (x, t) -> t * ones(Vec{1}), [1]))
    close!(dbc)
    Ferrite.update!(dbc, 0.0)
    return dbc
end;

function calculate_element_volume(cell, cellvalues_u, ue)
    reinit!(cellvalues_u, cell)
    evol::Float64 = 0.0
    for qp in 1:getnquadpoints(cellvalues_u)
        dΩ = getdetJdV(cellvalues_u, qp)
        ∇u = function_gradient(cellvalues_u, qp, ue)
        F = one(∇u) + ∇u
        J = det(F)
        evol += J * dΩ
    end
    return evol
end;

function calculate_volume_deformed_mesh(w, dh::DofHandler, cellvalues_u)
    evol::Float64 = 0.0
    for cell in CellIterator(dh)
        global_dofs = celldofs(cell)
        nu = getnbasefunctions(cellvalues_u)
        global_dofs_u = global_dofs[1:nu]
        ue = w[global_dofs_u]
        δevol = calculate_element_volume(cell, cellvalues_u, ue)
        evol += δevol
    end
    return evol
end;

function assemble_element!(Ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
    # Reinitialize cell values, and reset output arrays
    ublock, pblock = 1, 2
    reinit!(cellvalues_u, cell)
    reinit!(cellvalues_p, cell)
    fill!(Ke, 0.0)
    fill!(fe, 0.0)

    n_basefuncs_u = getnbasefunctions(cellvalues_u)
    n_basefuncs_p = getnbasefunctions(cellvalues_p)

    for qp in 1:getnquadpoints(cellvalues_u)
        dΩ = getdetJdV(cellvalues_u, qp)
        # Compute deformation gradient F
        ∇u = function_gradient(cellvalues_u, qp, ue)
        p = function_value(cellvalues_p, qp, pe)
        F = one(∇u) + ∇u

        # Compute first Piola-Kirchhoff stress and tangent modulus
        ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p = constitutive_driver(F, p, mp)

        # Loop over the `u`-test functions to calculate the `u`-`u` and `u`-`p` blocks
        for i in 1:n_basefuncs_u
            # gradient of the test function
            ∇δui = shape_gradient(cellvalues_u, qp, i)
            # Add contribution to the residual from this test function
            fe[BlockIndex((ublock), (i))] += (∇δui ⊡ ∂Ψ∂F) * dΩ

            ∇δui∂S∂F = ∇δui ⊡ ∂²Ψ∂F²
            for j in 1:n_basefuncs_u
                ∇δuj = shape_gradient(cellvalues_u, qp, j)

                # Add contribution to the tangent
                Ke[BlockIndex((ublock, ublock), (i, j))] += (∇δui∂S∂F ⊡ ∇δuj) * dΩ
            end
            # Loop over the `p`-test functions
            for j in 1:n_basefuncs_p
                δp = shape_value(cellvalues_p, qp, j)
                # Add contribution to the tangent
                Ke[BlockIndex((ublock, pblock), (i, j))] += (∂²Ψ∂F∂p ⊡ ∇δui) * δp * dΩ
            end
        end
        # Loop over the `p`-test functions to calculate the `p-`u` and `p`-`p` blocks
        for i in 1:n_basefuncs_p
            δp = shape_value(cellvalues_p, qp, i)
            fe[BlockIndex((pblock), (i))] += (δp * ∂Ψ∂p) * dΩ

            for j in 1:n_basefuncs_u
                ∇δuj = shape_gradient(cellvalues_u, qp, j)
                Ke[BlockIndex((pblock, ublock), (i, j))] += ∇δuj ⊡ ∂²Ψ∂F∂p * δp * dΩ
            end
            for j in 1:n_basefuncs_p
                δp = shape_value(cellvalues_p, qp, j)
                Ke[BlockIndex((pblock, pblock), (i, j))] += δp * ∂²Ψ∂p² * δp * dΩ
            end
        end
    end
    return
end;

function assemble_global!(
        K::SparseMatrixCSC, f, cellvalues_u::CellValues,
        cellvalues_p::CellValues, dh::DofHandler, mp::NeoHooke, w
    )
    nu = getnbasefunctions(cellvalues_u)
    np = getnbasefunctions(cellvalues_p)

    # start_assemble resets K and f
    fe = BlockedArray(zeros(nu + np), [nu, np]) # local force vector
    ke = BlockedArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix

    assembler = start_assemble(K, f)
    # Loop over all cells in the grid
    for cell in CellIterator(dh)
        global_dofs = celldofs(cell)
        global_dofsu = global_dofs[1:nu] # first nu dofs are displacement
        global_dofsp = global_dofs[(nu + 1):end] # last np dofs are pressure
        @assert size(global_dofs, 1) == nu + np # sanity check
        ue = w[global_dofsu] # displacement dofs for the current cell
        pe = w[global_dofsp] # pressure dofs for the current cell
        assemble_element!(ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
        assemble!(assembler, global_dofs, ke, fe)
    end
    return
end;

function solve(interpolation_u, interpolation_p)

    # import the mesh
    grid = importTestGrid()

    # Material parameters
    μ = 1.0
    λ = 1.0e4 * μ
    mp = NeoHooke(μ, λ)

    # Create the DofHandler and CellValues
    dh = create_dofhandler(grid, interpolation_u, interpolation_p)
    cellvalues_u, cellvalues_p, facetvalues_u = create_values(interpolation_u, interpolation_p)

    # Create the DirichletBCs
    dbc = create_bc(dh)

    # Pre-allocation of vectors for the solution and Newton increments
    _ndofs = ndofs(dh)
    w = zeros(_ndofs)
    ΔΔw = zeros(_ndofs)
    apply!(w, dbc)

    # Create the sparse matrix and residual vector
    K = allocate_matrix(dh)
    f = zeros(_ndofs)

    # We run the simulation parameterized by a time like parameter. `Tf` denotes the final value
    # of this parameter, and Δt denotes its increment in each step
    Tf = 2.0
    Δt = 0.1
    NEWTON_TOL = 1.0e-8

    pvd = paraview_collection("hyperelasticity_incomp_mixed")
    for (step, t) in enumerate(0.0:Δt:Tf)
        # Perform Newton iterations
        Ferrite.update!(dbc, t)
        apply!(w, dbc)
        newton_itr = -1
        prog = ProgressMeter.ProgressThresh(NEWTON_TOL; desc = "Solving @ time $t of $Tf;")
        fill!(ΔΔw, 0.0)
        while true
            newton_itr += 1
            assemble_global!(K, f, cellvalues_u, cellvalues_p, dh, mp, w)
            norm_res = norm(f[Ferrite.free_dofs(dbc)])
            apply_zero!(K, f, dbc)
            # Only display output at specific load steps
            if t % (5 * Δt) == 0
                ProgressMeter.update!(prog, norm_res; showvalues = [(:iter, newton_itr)])
            end
            if norm_res < NEWTON_TOL
                break
            elseif newton_itr > 30
                error("Reached maximum Newton iterations, aborting")
            end
            # Compute the incremental `dof`-vector (both displacement and pressure)
            ΔΔw .= K \ f

            apply_zero!(ΔΔw, dbc)
            w .-= ΔΔw
        end

        # Save the solution fields
        VTKGridFile("hyperelasticity_incomp_mixed_$step", grid) do vtk
            write_solution(vtk, dh, w)
            pvd[t] = vtk
        end
    end
    vtk_save(pvd)
    vol_def = calculate_volume_deformed_mesh(w, dh, cellvalues_u)
    print("Deformed volume is $vol_def")
    return vol_def
end;

quadratic_u = Lagrange{RefTetrahedron, 2}()^3
linear_p = Lagrange{RefTetrahedron, 1}()
vol_def = solve(quadratic_u, linear_p)

This page was generated using Literate.jl.