Not converging when applying dirichlet BC instead of Force

Hello everyone,

in the following MWE I just want to solve the displacement of a rectangle and is based on the following example: Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation

The code works fine when I impose a force on the top edge. However the same code does not converge when I try to define a displacement controlling dirichlet BC (here called “bctopy”) instead of the force.

MWE:

from dolfin import *
from mshr import *
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)

domain = Rectangle(Point(0,0), Point(1,1))
nxy = 40
mesh = generate_mesh(domain, nxy)

# elastic parameters
kappa = Constant(71245)
mu    = Constant(27120)
lambda_= Constant(kappa - 2./3.*mu)

V = VectorFunctionSpace(mesh, "CG", 2)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=2, quad_scheme='default')
W = FunctionSpace(mesh, We)

sig = Function(W)
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
du_ = TrialFunction(V)
v = TestFunction(V)

# Boundary conditions for mesh D2
class bottom(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1e-13
        return abs(x[1]) < tol and on_boundary
    
class top(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1e-13
        return abs(x[1]-1) < tol and on_boundary

Top = top()
Bottom = bottom()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
Bottom.mark(boundaries,1)
Top.mark(boundaries,2)

loading  = Expression("t", t = 0.0, degree=1)
bcbot = DirichletBC(V, Constant((0.0,0.0)), boundaries, 1)
bctopx = DirichletBC(V.sub(0), Constant(0.0), boundaries, 2)
bctopy = DirichletBC(V.sub(1), loading, boundaries, 2)
bc_u  = [bcbot, bctopx, bctopy]

n = FacetNormal(mesh)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
def F_ext(v):
    return loading*dot(n, v)*ds(2)

# -----------------------------
# Constitutive relation update
# -----------------------------
def eps(v):
    return sym(grad(v))

def sigma(eps_el):
    return lambda_*tr(eps_el)*Identity(2) + 2*mu*eps_el

sig = sigma(eps(u))

# --------------------------------------------
# Global problem and Newton-Raphson procedure
# --------------------------------------------
metadata = {"quadrature_degree": 2, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)

a_Newton = inner(eps(v), sigma(eps(du_)))*dxm
res = -inner(eps(v), sig)*dxm# + F_ext(v)

file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

Nitermax, tol = 50, 1e-8  # parameters of the Newton-Raphson procedure
Uincr = 9E-2
Nincr = 10
j = 0
t = 0

while j < Nincr:
    j += 1
    t += Uincr
    loading.t = t
    A, Res = assemble_system(a_Newton, res, bc_u)
    nRes0 = Res.norm("l2")
    nRes = nRes0
    print("Increment:", str(j))
    niter = 0
    while nRes/nRes0 > tol and niter < Nitermax:
        solve(A, du.vector(), Res, "mumps")
        u.assign(u+du)
        A, Res = assemble_system(a_Newton, res, bc_u)
        nRes = Res.norm("l2")
        print("    Residual:", nRes)
        niter += 1
    
    file_results.write(u, t)

Dolfin version (installed via Docker): 2019.1.0
Ubunutu version: 20.04

Does anyone have an idea what the mistake is?
Your help is much appreciated.

Best,
Dejan

1 Like

Dear future readers,

the solution for me was to switch to Fenicsx and to use the following custom Newton Solver
https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html