Phase field modeling

New to Julia
1

Hi everyone. I am doing a phase field modeling for fracture in julia using gridap. I have done two different method:1- is the phase field modeling by applying displacement on a boundary incrementally which works. 2- i am trying to do the same but not defining displacement on boundary. I have defined force implemented on boundary incrementally which is not working. I couldn’t understand why it is not working as every thing is the same. could you please tell me if you have any recommendation. thank you.

function   stepPhaseField(fem_params, uh_in,ΨPlusPrev_in)
    A_PF(s,φ,uh_in,ΨPlusPrev_in) =Gc*ls*∇(φ)⋅∇(s)+ (2*ΨPlusPrev_in*s*φ) + (Gc/ls)*s*φ
    a_PF(s,φ) = ∫(A_PF(s,φ,uh_in,ΨPlusPrev_in))fem_params.dΩ
    b_PF(φ) =∫((Gc/ls)*φ)fem_params.dΩ
    op_PF = AffineFEOperator(a_PF,b_PF ,U_PF ,V0_PF)
    sh_out = solve(op_PF)
    return  get_free_dof_values(sh_out)
end
function  stepDisp(fem_params,pth, uh_in,sh_in, vApp)
    g = VectorValue(0, -vApp)
    A_Disp(u,v,pth,uh_in,sh_in) =  ((p->Em(p))∘pth) * (ε(v) ⊙ (σfun∘(ε(u), ε(uh_in), sh_in)))
    a_Disp(u,v) =∫(A_Disp(u,v,pth,uh_in,sh_in))fem_params.dΩ
    b_Disp(v) = ∫(v ⋅ g)fem_params.dΓ_Load
    op_Disp = AffineFEOperator(a_Disp ,b_Disp ,fem_params.U0_Disp ,fem_params.V0_Disp)
    uh_out = solve(op_Disp)
    return  get_free_dof_values(uh_out)
end
vApp = 0.0
delv = 1.5e-2
const vAppMax = 1.0
vAppmid = 5e-1
delvmid = 1.5e-2
co = Int(round(((vAppmid - vApp)/ delv) + ((vAppMax - vAppmid)/ delvmid)))
innerMax = 30
count = 0
ψPlusPrev = CellState(0.0,dΩ)
p0 = ones(fem_params.np)
pf_vec = Filter(p0;r, fem_params)
pfh = FEFunction(fem_params.Pf, pf_vec)
pth = (pf -> Threshold(pf; β, η)) ∘ pfh
sPrev = CellState(1.0,dΩ)
sh = project(sPrev, model ,dΩ,order)

while vApp .< vAppMax
    count = count .+ 1
    if vApp >= vAppmid
        delv = delvmid
    end
    vApp = vApp .+ delv
    print("\n step", float(vApp))
   for  inner = 1: innerMax
        ψhPlusPrev = project(ψPlusPrev ,model ,dΩ,order)
        RelErr = abs(sum(∫(Gc*ls*∇(sh)⋅∇(sh) + 2*ψhPlusPrev*sh*sh + (Gc/ls)*sh*sh)fem_params.dΩ-∫( (Gc/ls)*sh)fem_params.dΩ))/abs(sum(∫( (Gc/ls)*sh)fem_params.dΩ))
        sh_vec = stepPhaseField(fem_params,uh,ψhPlusPrev)
        sh = FEFunction(fem_params.U_PF, sh_vec)
        u_vec = stepDisp(fem_params,pth, uh, sh, vApp)
        uh = FEFunction(fem_params.U0_Disp, u_vec)
        ψhPos_in = ψPos∘(ε(uh))
        update_state!( new_EnergyState ,ψPlusPrev ,ψhPos_in)
        if   RelErr  < 1e-8
            break
        end
  end
end
writevtk(Ω,"resultstph",cellfields=["uh"=>uh, "sh"=>sh])