Plasticity in fenics

43

Hi, I am also trying the same thing and I changed the Newton loop accordingly like this -


Nitermax, tol = 200, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))

for (i, t) in enumerate(load_steps):
    loading.t = t
    A, Res = assemble_system(a, L, bc_u)
    for bc in bc_u:
        bc.homogenize()
    nRes0 = Res.norm("l2")
    nRes = nRes0
    print("Increment:", str(i+1))    
    niter = 0
    while nRes/nRes0 > tol and niter < Nitermax:
        solve(A, du.vector(), Res, "mumps")
        Du.assign(Du+du)
        deps = eps(Du)
        sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
        local_project(sig_, W, sig)
        local_project(n_elas_, W, n_elas)
        local_project(beta_, W0, beta)
        A, Res = assemble_system(a, L, bc_u)
        nRes = Res.norm("l2")
        print("    Residual:", nRes)
        niter += 1
    u.assign(u+Du)
    p.assign(p+local_project(dp_, W0))
    sig_old.assign(sig)    
    p_avg.assign(project(p, P0))
    results[i+1, :] = (u(0.5, 1.)[1], t)

where my boundary conditions are -

# Boundary conditions
class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1., 1e-8)
class bot(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0., 1e-8)

boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)

top().mark(boundaries, 1)
bot().mark(boundaries, 2)

loading = Expression("t", t = 0.0, degree=1)

bc_top = DirichletBC(V.sub(1), loading, boundaries, 1)
bc_bottom = DirichletBC(V, Constant((0., 0.)), boundaries, 2)

bc_u = [bc_bottom, bc_top]

But the solver is diverging so I don’t know where am I making mistake. Can you please tell me how were you able to do it ? I have made a separate post regarding the same problem Iterative solver diverging . Any pointers would be highly appreciable.
Thanks.