Nearly-Incompressible Hyperelasticity throwing error

[bhaveshshrimali]

Hi,

I came back to an old code that I was working on a while ago. It is nearly-incompressible hyperelasticity (λ/μ >> 1). But it is failing with the stacktrace pointing to a problem in constitutive relation (the function aJ below).. here is the code

using Gridap
using LineSearches: BackTracking
using BenchmarkTools

const λ = 1000.0
const μ = 1.0

qmodel = CartesianDiscreteModel((0, 1, 0, 1, 0, 1), (5, 5, 5))
model = simplexify(qmodel)
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"rightFace",[2, 4, 6, 8, 14, 16, 18, 20, 26])
add_tag_from_tags!(labels,"leftFace",[1, 3, 5, 7, 13, 15, 17, 19, 25])
add_tag_from_tags!(labels,"bottomFace",[1, 2, 5, 6, 9, 11, 17, 18, 23])
add_tag_from_tags!(labels,"backFace",[1, 2, 3, 4, 9, 10, 13, 14, 21])


degree=3
trian = Triangulation(model)
dΩ = Measure(trian, degree)

writevtk(trian, "domainGridap")

# helper functions
F(∇u) = one(∇u) + ∇u
J(F) = det(F)

function S(∇u, p)
    Fs = F(∇u)
    Js = J(Fs)
    μ * Fs + p * Js * inv(Fs)'
end


# this is the problematic function
function aJ(∇u, p)
    Js = (λ + p + sqrt((λ + p)^2 + 4 * λ * μ )) / (2 * λ)   # specifically this 
    dJdp = (λ + p + sqrt(4*λ*μ + (λ+ p)^2))/(2*λ*sqrt(4*λ*μ + (λ + p)^2))
    J(F(∇u)) - (Js + (p + μ/Js - λ * (Js - 1))* dJdp)
end

function res(Wh, Vh)
    u, p = Wh
    v, q = Vh
    (S∘(∇(u), p) ⊙ ∇(v) + aJ(∇(u), p) * q) * dΩ
end

nls = NLSolver(
    show_trace=true,
    method=:newton,
    linesearch=BackTracking()
)

reffe_u = ReferenceFE(lagrangian, VectorValue{3,Float64},2)
reffe_p = ReferenceFE(lagrangian,Float64,1)

Vu = TestFESpace(
    model, reffe_u, conformity=:H1,
    dirichlet_tags = ["leftFace", "rightFace", "bottomFace", "backFace"],
    dirichlet_masks = [
        (true, false, false), (true, false, false), (false, true, false), (false, false, true)
    ]
)

Vp = TestFESpace(
    model,reffe_p,conformity=:H1
)

solver = FESolver(nls)

function main(init_u, init_p, disp_x)

    g0 = VectorValue(
        0.0, 0.0, 0.0
    )

    g1 = VectorValue(
        disp_x, 0.0, 0.0
    )

    Uh = TrialFESpace(Vu, [g0, g1, g0, g0])
    Qh = TrialFESpace(Vp)

    W1 = MultiFieldFESpace([Uh, Qh])
    W2 = MultiFieldFESpace([Vu, Vp])

    op = FEOperator(res, W1, W2)

    uh = FEFunction(Uh, init_u)
    ph = FEFunction(Qh, init_p)

    uh, ph = solve(solver, op)
    #  = wh;
    # @assert(ph.free_values[1] ≈ -0.5) #exact solution
    return get_free_dof_values(uh), get_free_dof_values(ph)
end

function runs()
    nsteps = 10;
    disp_max = 2;
    u0 = zeros(Float64, num_free_dofs(Vu));
    p0 = ones(Float64, num_free_dofs(Vp));
    cache=nothing;
    for step in 1:nsteps
        disp_x = step * disp_max / nsteps
        println("\n+++ Solving for disp_x $disp_x")
        u0, p0 = main(u0, p0, disp_x)
    end
end

runs()

The stacktrace is

ERROR: MethodError: no method matching ^(::Gridap.CellData.OperationCellField{ReferenceDomain}, ::Int64)

Closest candidates are:
  ^(::Union{AbstractChar, AbstractString}, ::Integer)
   @ Base strings\basic.jl:733
  ^(::Rational, ::Integer)
   @ Base rational.jl:491
  ^(::Complex{<:AbstractFloat}, ::Integer)
   @ Base complex.jl:869
  ...

Stacktrace:
  [1] literal_pow
    @ .\intfuncs.jl:338 [inlined]
  [2] aJ(∇u::Gridap.CellData.GenericCellField{ReferenceDomain}, p::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}})
    @ Main d:\juliaPractice\Gridap\hyperelasticity.jl:34

The weak form (with a slight abuse of notation) looks like the following

∫ (S(∇u, p) ⊙ ∇(v) + fJ⋆(∇u, p) * q) dΩ = 0 ∀ (v, q)

where S is the first piola kirchhoff stress and fJ⋆(∇u, p) is a nonlinear function in the constitutive relation as defined above, namely

function aJ(∇u, p)
    Js = (λ + p + sqrt((λ + p)^2 + 4 * λ * μ )) / (2 * λ)
    dJdp = (λ + p + sqrt(4*λ*μ + (λ+ p)^2))/(2*λ*sqrt(4*λ*μ + (λ + p)^2))
    J(F(∇u)) - (Js + (p + μ/Js - λ * (Js - 1))* dJdp)
end

Could anyone point me to the cause of the problem and how should I define the above function ?

[bhaveshshrimali]

It seems it was

∫(S∘(∇(u), p) ⊙ ∇(v) + aJ∘(∇(u), p) ⊙ q)*dΩ