# 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 *
import warnings

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)
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)

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

n = FacetNormal(mesh)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
def F_ext(v):

# -----------------------------
# Constitutive relation update
# -----------------------------
def eps(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
# --------------------------------------------

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
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?