Dof ordering for multi-field dof handlers

[bhaveshshrimali]

Hi,

In trying to debug #364, I am trying to understand the dof ordering for dof handlers that have more than one field. Consider the classical Taylor-Hood approximation for incompressible elasticity (displacement-pressure). From what I understand by taking a look at the incompressible elasticity and DOFs the dofs for each element are ordered as, all displacement (u) and then all pressure (p) in the order: nodes -> facets -> cells, which also explains the block structure of the element matrices.

Since my problem is nonlinear I need the dofs at each newton iterate to reassemble the jacobian and residual. My current idea is: get all the dofs for the corresponding cell by celldofs(cell) and then index into the dofs based on the ordering explained above, namely

@timeit "assemble" 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] # element dofs displacement
        pe = w[global_dofsp] # element dofs pressure
        @timeit "element assemble" assemble_element!(ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
        assemble!(assembler, global_dofs, fe, ke)
end

where w is my global dofs-vector. Is this correct? Thanks!

[kimauth]

Hi Bhavesh, this indeed sounds correct to me! You can trust the dofs of one field to be a "block", there even is a function to get the dof-offset to a given field here. (Though considering that this goes looking for the fieldname every time, it might be worth to compute it once and store it. If you're working with the DofHandler it won't ever change.) Regarding the ordering within each field: You're correct, but in 3D there also is edgedofs, which come after vertexdofs (so in general we have vertexdofs - edgedofs - facedofs- celldofs). You can find the vertex/edge/face definitions for the defined celltypes here.

[bhaveshshrimali]

Thanks Kim! Things seem to be working now. I had to print the assembled jacobian to notice that I wasn't assembling the jacobian properly and hence the singular system in #364. In fact except for this bit (element dofs for displacement and pressure) I don't think I had to change anything to merge together two examples to come up with the code.

Here is the code for incompressible hyperelasticity (and in the process my first successful attempt at using Ferrite :) ), for anyone who might find it useful to start off with in the future.

using Ferrite, Tensors, TimerOutputs, ProgressMeter
import KrylovMethods, IterativeSolvers
using BlockArrays, SparseArrays, LinearAlgebra
using FerriteGmsh, BenchmarkTools

struct NeoHooke
    μ::Float64
    λ::Float64
end

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

function create_values(interpolation_u, interpolation_p)
    # quadrature rules
    qr      = QuadratureRule{3,RefTetrahedron}(4)
    face_qr = QuadratureRule{2,RefTetrahedron}(4)

    # geometric interpolation
    interpolation_geom = Lagrange{3,RefTetrahedron,1}()

    # cell and facevalues for u
    cellvalues_u = CellVectorValues(qr, interpolation_u, interpolation_geom)
    facevalues_u = FaceVectorValues(face_qr, interpolation_u, interpolation_geom)

    # cellvalues for p
    cellvalues_p = CellScalarValues(qr, interpolation_p, interpolation_geom)

    return cellvalues_u, cellvalues_p, facevalues_u
end;

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

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

function create_bc(dh)
    dbc = ConstraintHandler(dh)
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myLeft"), (x,t) -> zero(Vec{1}), [1]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myBottom"), (x,t) -> zero(Vec{1}), [2]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myBack"), (x,t) -> zero(Vec{1}), [3]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myRight"), (x,t) -> t*ones(Vec{1}), [1]))
    close!(dbc)
    t = 1.0
    Ferrite.update!(dbc, t)
    return dbc
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 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)

    @inbounds for qp in 1:getnquadpoints(cellvalues_u)
        dΩ = getdetJdV(cellvalues_u, qp)

        ## Compute deformation gradient F and right Cauchy-Green tensor C
        ∇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
        ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p = constitutive_driver(F, p, mp)

        ## Loop over test functions
        @inbounds for i in 1:n_basefuncs_u
            ## Test function gradient
            ∇δ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² # Hoisted computation
            @inbounds 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

            @inbounds 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

        @inbounds for i in 1:n_basefuncs_p
            δp = shape_value(cellvalues_p, qp, i)
            fe[BlockIndex((pblock), (i))] += ( δp * ∂Ψ∂p) * dΩ

            @inbounds 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
            @inbounds 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
    # symmetrize_lower!(Ke)
end;

function symmetrize_lower!(K)
    for i in 1:size(K,1)
        for j in i+1:size(K,1)
            K[i,j] = K[j,i]
        end
    end
end;

# Assembling global residual and tangent
function assemble_global!(K::SparseMatrixCSC, f, cellvalues_u::CellVectorValues{dim}, 
    cellvalues_p::CellScalarValues{dim}, dh::DofHandler, mp::NeoHooke, w) where {dim}

    nu = getnbasefunctions(cellvalues_u)
    np = getnbasefunctions(cellvalues_p)

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

    assembler = start_assemble(K, f)
    ## Loop over all cells in the grid
    @timeit "assemble" 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] # element dofs displacement
        pe = w[global_dofsp] # element dofs pressure
        @timeit "element assemble" assemble_element!(ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
        assemble!(assembler, global_dofs, fe, ke)
        # return K, f;
        # break
    end
end;

# Define a main function, with a loop for Newton iterations

function solve(interpolation_u, interpolation_p)
    reset_timer!()
    grid = importTestGrid()

    ## Material parameters
    μ = 1.
    λ = 1.E3 * μ
    mp = NeoHooke(μ, λ)

    ## Finite element base
    dh = create_dofhandler(grid, interpolation_u, interpolation_p)
    cellvalues_u, cellvalues_p, facevalues_u = create_values(interpolation_u, interpolation_p)

    dbc = create_bc(dh)

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

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

    ## Perform Newton iterations
    newton_itr = -1
    NEWTON_TOL = 1e-12
    prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:")

    while true; newton_itr += 1
        w .= wn .+ Δw # Current guess
        assemble_global!(K, f, cellvalues_u, cellvalues_p, dh, mp, w)
        # return K, f;
        # break;
        norm_res = norm(f[Ferrite.free_dofs(dbc)])
        apply_zero!(K, f, dbc)
        ProgressMeter.update!(prog, norm_res; showvalues = [(:iter, newton_itr)])

        if norm_res < NEWTON_TOL
            break
        elseif newton_itr > 30
            error("Reached maximum Newton iterations, aborting")
        end

        ## Compute increment using cg! from IterativeSolvers.jl
        # @timeit "linear solve (KrylovMethods.cg)" ΔΔw′, flag, relres, iter, resvec = KrylovMethods.cg(K, f; maxIter = 1000)
        # @assert flag == 0
        # @timeit "linear solve (IterativeSolvers.cg!)" IterativeSolvers.cg!(ΔΔw, K, f; maxiter=1000)
        ΔΔw .= K\f;

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

    ## Save the solution
    @timeit "export" begin
        vtk_grid("hyperelasticity_incomp_mixed", dh) do vtkfile
            vtk_point_data(vtkfile, dh, w)
        end
    end

    print_timer(title = "Analysis with $(getncells(grid)) elements", linechars = :ascii)
    return w
end

quadratic = Lagrange{3, RefTetrahedron, 2}()
linear = Lagrange{3, RefTetrahedron, 1}()
dofs = solve(quadratic, linear)

Contour

[termi-official]

Hi @bhaveshshrimali , do you mind if I recycle your code either for a new Ferrite example or to update the current Hyperelasticity example?

[bhaveshshrimali]

Hi @termi-official. Sure, please feel free to. I just need to add some comments on the formulation, and a meaningful test (like computing the deformed volume). Will do so in this tomorrow, if that suits you. In fact, much of it would carry in a straightforward manner from ex36 in scikit-fem which I wrote quite sometime back.

[bhaveshshrimali]

Hi @termi-official . Here are the two helper functions

function calculate_element_volume(cell, cellvalues_u, ue)
    reinit!(cellvalues_u, cell)
    evol=0.0;
    @inbounds 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
    @inbounds 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

and doing

isapprox(vol_def, 1.0, atol=1E-3)

gives true. Let me check things again just to be completely sure about it. And rewrite the formulation a little bit. It should be good to go from thereon. If you want, I can also open a PR for the same.

[bhaveshshrimali]

^ The mesh is that of a unit cube (undeformed volume = 1.0). And the formulation is "nearly incompressible" hyperelasticity. I will write things down and post here for explicitness.

[termi-official]

Thanks for the quick anwser @bhaveshshrimali ! Sure, if you want then just go ahead and open a PR. Then we have the code saved in the docs, which should be easier to find for everyone. :)

[bhaveshshrimali]

Thanks. Sounds good. Let me put together a rough template and then we can just take it from there.

[bhaveshshrimali]

This is a rough template of the code @termi-official . Still some cleaning required. I was traveling over the last week, will update with a PR this weekend, with proper comments and a formulation. I added the bit with incrementally solving the equations (with a pseudo time-step like parameter) as is usually done with nonlinear problems.

Code

using Ferrite, Tensors, TimerOutputs, ProgressMeter
import KrylovMethods, IterativeSolvers
using BlockArrays, SparseArrays, LinearAlgebra
using BenchmarkTools

struct NeoHooke
    μ::Float64
    λ::Float64
end

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

function create_values(interpolation_u, interpolation_p)
    # quadrature rules
    qr      = QuadratureRule{3,RefTetrahedron}(4)
    face_qr = QuadratureRule{2,RefTetrahedron}(4)

    # geometric interpolation
    interpolation_geom = Lagrange{3,RefTetrahedron,1}()

    # cell and facevalues for u
    cellvalues_u = CellVectorValues(qr, interpolation_u, interpolation_geom)
    facevalues_u = FaceVectorValues(face_qr, interpolation_u, interpolation_geom)

    # cellvalues for p
    cellvalues_p = CellScalarValues(qr, interpolation_p, interpolation_geom)

    return cellvalues_u, cellvalues_p, facevalues_u
end;

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

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

function create_bc(dh)
    dbc = ConstraintHandler(dh)
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myLeft"), (x,t) -> zero(Vec{1}), [1]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myBottom"), (x,t) -> zero(Vec{1}), [2]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myBack"), (x,t) -> zero(Vec{1}), [3]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myRight"), (x,t) -> t*ones(Vec{1}), [1]))
    close!(dbc)
    Ferrite.update!(dbc, 0.0)
    return dbc
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 calculate_element_volume(cell, cellvalues_u, ue)
    reinit!(cellvalues_u, cell)
    evol=0.0;
    @inbounds 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
    @inbounds 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)

    @inbounds for qp in 1:getnquadpoints(cellvalues_u)
        dΩ = getdetJdV(cellvalues_u, qp)

        ## Compute deformation gradient F and right Cauchy-Green tensor C
        ∇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
        ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p = constitutive_driver(F, p, mp)

        ## Loop over test functions
        @inbounds for i in 1:n_basefuncs_u
            ## Test function gradient
            ∇δ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² # Hoisted computation
            @inbounds 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

            @inbounds 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

        @inbounds for i in 1:n_basefuncs_p
            δp = shape_value(cellvalues_p, qp, i)
            fe[BlockIndex((pblock), (i))] += ( δp * ∂Ψ∂p) * dΩ

            @inbounds 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
            @inbounds 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
    # symmetrize_lower!(Ke)
end;

# Assembling global residual and tangent
function assemble_global!(K::SparseMatrixCSC, f, cellvalues_u::CellVectorValues{dim}, 
    cellvalues_p::CellScalarValues{dim}, dh::DofHandler, mp::NeoHooke, w) where {dim}

    nu = getnbasefunctions(cellvalues_u)
    np = getnbasefunctions(cellvalues_p)

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

    assembler = start_assemble(K, f)
    ## Loop over all cells in the grid
    @timeit "assemble" 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] # element dofs displacement
        pe = w[global_dofsp] # element dofs pressure
        @timeit "element assemble" assemble_element!(ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
        assemble!(assembler, global_dofs, fe, ke)
        # return K, f;
        # break
    end
end;

# Define a main function, with a loop for Newton iterations

function solve(interpolation_u, interpolation_p)
    reset_timer!()
    grid = importTestGrid()

    ## Material parameters
    μ = 1.
    λ = 1.E3 * μ
    mp = NeoHooke(μ, λ)

    ## Finite element base
    dh = create_dofhandler(grid, interpolation_u, interpolation_p)
    cellvalues_u, cellvalues_p, facevalues_u = create_values(interpolation_u, interpolation_p)

    dbc = create_bc(dh)

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

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

    Tf = 2.0;
    Δt = 0.01;
    NEWTON_TOL = 1e-8

    pvd = paraview_collection("hyperelasticity_incomp_mixed.pvd");
    for t ∈ 0.0:Δt:Tf
        ## Perform Newton iterations
        print("Time is $t")
        Ferrite.update!(dbc, t)
        apply!(w, dbc)
        newton_itr = -1
        prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:")
        fill!(ΔΔw, 0.0);
        while true; newton_itr += 1
            # w .= wn .+ Δw # Current guess
            assemble_global!(K, f, cellvalues_u, cellvalues_p, dh, mp, w)
            # return K, f;
            # break;
            norm_res = norm(f[Ferrite.free_dofs(dbc)])
            apply_zero!(K, f, dbc)
            ProgressMeter.update!(prog, norm_res; showvalues = [(:iter, newton_itr)])

            if norm_res < NEWTON_TOL
                break
            elseif newton_itr > 30
                error("Reached maximum Newton iterations, aborting")
            end

            ## Compute increment using cg! from IterativeSolvers.jl
            # @timeit "linear solve (KrylovMethods.cg)" ΔΔw′, flag, relres, iter, resvec = KrylovMethods.cg(K, f; maxIter = 1000)
            # @assert flag == 0
            # @timeit "linear solve (IterativeSolvers.cg!)" IterativeSolvers.cg!(ΔΔw, K, f; maxiter=1000)
            ΔΔw .= K\f;

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

        ## Save the solution
        @timeit "export" begin
            vtk_grid("hyperelasticity_incomp_mixed-$t", dh) do vtkfile
                vtk_point_data(vtkfile, dh, w)
                vtk_save(vtkfile)
                pvd[t] = vtkfile
            end
        end

        print_timer(title = "Analysis with $(getncells(grid)) elements", linechars = :ascii)
    end;
    vtk_save(pvd);
    vol_def = calculate_volume_deformed_mesh(w, dh, cellvalues_u);
    print("Deformed volume is $vol_def")
end

quadratic = Lagrange{3, RefTetrahedron, 2}()
linear = Lagrange{3, RefTetrahedron, 1}()
solve(quadratic, linear);
# vol_def = calculate_volume_deformed_mesh(dof_vector, dh, cellvalues_u)

isapprox(vol_def, 1.0, atol=1E-3)

[bhaveshshrimali]

Sort of like this (there is some rendering issue but I think the gist is clear):

[termi-official]

Thanks @bhaveshshrimali , looks already pretty nice! No need to rush if you do not have the time atm.