Plane stress hyperelasticity

[Aminofa70]

Dear Ferrite developer
The package has FEM for hyperelasticty in 3D case.
However, the challenge here is plane stress case. For this case C_33 the component of right Cauchy-Green deformation tensor is not equal to 1 anymore and should be obtained when S (second Piola) is equal to zero.
I implemented the following code for neo-Hookean and plane stress. But there is error:
I would appreciate it if you could advice me.

ERROR: DomainError with -1223.1360357656133:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)

using Ferrite, Tensors, TimerOutputs, ProgressMeter, IterativeSolvers
# Define the NeoHookean material structure
struct NeoHooke
    μ::Float64
    λ::Float64
end

# Define the strain energy function
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

# Compute the second Piola-Kirchhoff stress tensor S and its derivative
function constitutive_driver(C, mp::NeoHooke)
    ∂²Ψ∂C², ∂Ψ∂C = Tensors.hessian(y -> Ψ(y, mp), C, :all)
    S = 2.0 * ∂Ψ∂C
    ∂S∂C = 2.0 * ∂²Ψ∂C²
    return S, ∂S∂C
end

# Newton-Raphson method to find C_33 such that S_33 = 0
# Newton-Raphson method to find C_33 such that S_33 = 0
function newton_raphson_C33(C, mp::NeoHooke; tol=1e-8, max_iter=50)
    C33 = C[3,3]  # Initial guess
    for iter in 1:max_iter
        # Create a new tensor with updated C33 (since Tensors.jl objects are immutable)
        C_new = Tensor{2, 3, Float64}([
            C[1,1] C[1,2] 0.0;
            C[2,1] C[2,2] 0.0;
            0.0    0.0    C33
        ])

        # Compute S and its derivative
        S, ∂S∂C = constitutive_driver(C_new, mp)

        S33 = S[3,3]
        dS33_dC33 = ∂S∂C[3,3,3,3]  # Extract derivative ∂S_33/∂C_33

        # Check convergence
        if abs(S33) < tol
            println("Converged in $iter iterations: C33 = $C33")
            return C33
        end

        # Newton-Raphson update
        C33 -= S33 / dS33_dC33
    end

    error("Newton-Raphson did not converge after $max_iter iterations.")
end

function assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN, tn)

    # Reinitialize cell values, and reset output arrays
    reinit!(cv, cell)
    fill!(ke, 0.0)
    fill!(ge, 0.0)

    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)

        F2D =  one(∇u) + ∇u
        C2D = tdot(F2D)  # C = F' ⋅ F
        # Extend C to 3D with initial guess for C33
        C33_Guess = 1.11
        

        # C33_Guess = max(1.0, 1.0 / (tr(C2D) + 1e-3))  # Ensure it's reasonable

        C3D = Tensor{2,3,Float64}([
            C2D[1,1] C2D[1,2] 0.0;
            C2D[2,1] C2D[2,2] 0.0;
            0.0      0.0      C33_Guess
        ])

        # Solve for C33 using Newton-Raphson
        C33_solution = newton_raphson_C33(C3D, mp)

        # Update C3D with computed C33
        C_updated = Tensor{2,3,Float64}([
            C2D[1,1] C2D[1,2] 0.0;
            C2D[2,1] C2D[2,2] 0.0;
            0.0      0.0      C33_solution
        ])

        F = Tensor{2,3,Float64}([
            F2D[1,1] F2D[1,2] 0.0;
            F2D[2,1] F2D[2,2] 0.0;
            0.0      0.0      sqrt(C33_solution)
        ])
        # F = one(∇u) + ∇u

        # C = tdot(F) # F' ⋅ F
        # # Compute stress and tangent
        # S, ∂S∂C = constitutive_driver(C, mp)
        S, ∂S∂C = constitutive_driver(C_updated, mp)
        P = F ⋅ S
        I = one(S)
        ∂P∂F =  otimesu(I, S) + 2 * otimesu(F, I) ⊡ ∂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_2D = shape_gradient(cv, qp, i)
         
            ∇δui = Tensor{2,3,Float64}([
                ∇δui_2D[1,1] ∇δui_2D[1,2] 0.0;
                ∇δui_2D[2,1] ∇δui_2D[2,2] 0.0;
                0.0      0.0      0.0;
            ])
            # Add contribution to the residual from this test function
            ge[i] += ( ∇δui ⊡ P) * dΩ

            ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation
            for j in 1:ndofs
                ∇δuj_2D = shape_gradient(cv, qp, j)
                ∇δuj = Tensor{2,3,Float64}([
                    ∇δuj_2D[1,1] ∇δuj_2D[1,2] 0.0;
                    ∇δuj_2D[2,1] ∇δuj_2D[2,2] 0.0;
                    0.0      0.0      0.0;
                ])
                # 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)
                t = tn 
                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

    
end;

function assemble_global!(K, g, dh, cv, fv, mp, u, ΓN, tn)
    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
    @timeit "assemble" for cell in CellIterator(dh)
        global_dofs = celldofs(cell)
        ue = u[global_dofs] # element dofs
        @timeit "element assemble" assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN,tn)
        assemble!(assembler, global_dofs, ke, ge)
    end
end;

function create_grid(Lx, Ly, nx, ny)
    corners = [
        Ferrite.Vec{2}((0.0, 0.0)), Ferrite.Vec{2}((Lx, 0.0)),
        Ferrite.Vec{2}((Lx, Ly)), Ferrite.Vec{2}((0.0, Ly))
    ]
    grid = Ferrite.generate_grid(Ferrite.Quadrilateral, (nx, ny), corners)
    addnodeset!(grid, "clamped", x -> x[1] ≈ 0.0)
    addfacetset!(grid, "traction", x -> x[1] ≈ Lx && norm(x[2] - 0.5) <= 0.05)
    return grid
end
# Function to create CellValues and FacetValues
function create_values()
    dim, order = 2, 1
    ip = Ferrite.Lagrange{Ferrite.RefQuadrilateral,order}()^dim
    qr = Ferrite.QuadratureRule{Ferrite.RefQuadrilateral}(2)
    qr_face = Ferrite.FacetQuadratureRule{Ferrite.RefQuadrilateral}(1)
    cell_values = Ferrite.CellValues(qr, ip)
    facet_values = Ferrite.FacetValues(qr_face, ip)
    return cell_values, facet_values
end

# Function to create DofHandler
function create_dofhandler(grid)
    dh = Ferrite.DofHandler(grid)
    Ferrite.add!(dh, :u, Ferrite.Lagrange{Ferrite.RefQuadrilateral,1}()^2)
    Ferrite.close!(dh)
    return dh
end

# Function to create Dirichlet boundary conditions
function create_bc(dh)
    ch = Ferrite.ConstraintHandler(dh)
    add!(ch, Dirichlet(:u, getnodeset(dh.grid, "clamped"), (x, t) -> [0.0, 0.0], [1, 2]))
    Ferrite.close!(ch)
    return ch
end

function solve()

    reset_timer!()
    # Define parameters for the plate and mesh
    Lx, Ly = 2.0, 1.0  # Plate dimensions
    nx, ny = 30, 20   # Number of elements along x and y
    grid = create_grid(Lx, Ly, nx, ny)  # Generate the grid

    # Material parameters
    E = 10.0
    ν = 0.3
    μ = E / (2(1 + ν))
    λ = (E * ν) / ((1 + ν) * (1 - 2ν))
    mp = NeoHooke(μ, λ)

    # Create DOF handler and constraints
    dh = create_dofhandler(grid)
    dbcs = create_bc(dh)
    cv, fv = create_values()

    ΓN = getfacetset(grid, "traction")


    # 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, dbcs)

    # Create sparse matrix and residual vector
    K = allocate_matrix(dh)
    g = zeros(_ndofs)
    n_timesteps = 10
    # traction_magnitude = 1.e7 * range(0.5, 1.0, length=n_timesteps)
    traction_magnitude =  range(0.1, 1.0, length=n_timesteps)

    
    for timestep in 1:n_timesteps

        t = timestep # actual time (used for evaluating d-bndc)
        tn = Vec((0.0, -traction_magnitude[timestep]))


        # Perform Newton iterations
        newton_itr = -1
        NEWTON_TOL = 1e-8
        NEWTON_MAXITER = 50
        prog = ProgressMeter.ProgressThresh(NEWTON_TOL; desc = "Solving:")

        Ferrite.update!(dbcs, t) # evaluates the D-bndc at time t
        Ferrite.apply!(u, dbcs)  # set the prescribed values in the solution vector
        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,tn)
            # Apply boundary conditions
            apply_zero!(K, g, dbcs)
            # Compute the residual norm and compare with tolerance
            normg = norm(g)
            ProgressMeter.update!(prog, normg; showvalues = [(:iter, newton_itr)])
            if normg < NEWTON_TOL
                break
            elseif newton_itr > NEWTON_MAXITER
                error("Reached maximum Newton iterations, aborting")
            end

            # Compute increment using conjugate gradients
            @timeit "linear solve" IterativeSolvers.cg!(ΔΔu, K, g; maxiter=1000)

            apply_zero!(ΔΔu, dbcs)
            Δu .-= ΔΔu
        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)

    end

    return u

end 
u = solve();