DOLFIN error in FEM with lagrange multiplier problem

General
1

Hello,

I made a code that solves for and plot u (displacement) and p (lagrange multiplier) within 2D body.
I think I set the weak form right, but keep getting DOLFIN related error.
Could anyone help me out with this? I attached my code below:

# Given constants (units are in GPa)
mu = 69.2308

# Define strain
def epsilon(u):
    return 0.5*(grad(u) + grad(u).T)

# Define stress
def sigma(u,p):
    return 2*mu*epsilon(u) + p*Identity(d)

# Read the mesh file into domain
mesh = Mesh()
XDMFFile('mymesh.xdmf').read(mesh)

# Define function space
V1 = VectorFunctionSpace(mesh, 'CG', 1)   # for displacement u
V2 = FunctionSpace(mesh, 'R', 0)   # for lagrange multiplier p
W = MixedFunctionSpace(V1,V2)

# Assign left/right boundaries
def left(x,on_boundary):
    return near(x[0],0) and on_boundary
def right(x,on_boundary):
    return near(x[0],3) and on_boundary

# Mark the boundary on the right
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(left).mark(boundary_subdomains, 1)
AutoSubDomain(right).mark(boundary_subdomains, 2)     # This is the boundary where the traction is applied
dss = ds(subdomain_data=boundary_subdomains)

# Dirichlet BC: u_x = u_y = 0 at the left
bc = DirichletBC(V1,Constant((0,0)),left)

# Define variational problem 
u,p = TrialFunctions(W)
v,q = TestFunctions(W)
d = u.geometric_dimension()  # space dimension of u
T = Constant((1,0))   # Traction BC on the right

# Set the weak form with the traction BC on the right edge (t = (1,0) GPa)
a = inner(sigma(u,p),epsilon(v))*dx - dot(q,div(u))*dx - dot(grad(q),div(sigma(u,p)))*dx
L = dot(T,v)*dss(2)

# Solve for u
w = Function(W)
solve(a == L, w, bc)
(u,p) = w.split()

# Plot u and p
plot(u)
p_project = project(p,V2)
plot(p)
2

The error message:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[16], line 51
     49 # Solve for u
     50 w = Function(W)
---> 51 solve(a == L, w, bc)
     52 (u,p) = w.split()
     54 # Plot u and p

File /usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py:233, in solve(*args, **kwargs)
    230 # Call variational problem solver if we get an equation (but not a
    231 # tolerance)
    232 elif isinstance(args[0], ufl.classes.Equation):
--> 233     _solve_varproblem(*args, **kwargs)
    235 # Default case, just call the wrapped C++ solve function
    236 else:
    237     if kwargs:

File /usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py:264, in _solve_varproblem(*args, **kwargs)
    262     solver = MixedLinearVariationalSolver(problem)
    263     solver.parameters.update(solver_parameters)
--> 264     solver.solve()
    266 else:
    267     # Create problem
    268     problem = LinearVariationalProblem(eq.lhs, eq.rhs, u, bcs,
    269                                        form_compiler_parameters=form_compiler_parameters)

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 10000 iterations (PETSc reason DIVERGED_ITS, residual norm ||r|| = 9.099157e-04).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------