ArityMismatch error in 3D elasticity + scalar field

Errors
1

I am doing elasticity assisted phase field calculations in FEniCS. The system is an anisotropic (single domain) sphere compressed uniaxially. I am using a mixed function space (a vector element for displacement and a scalar element for field variable. The variable can assume any value but only between 0 and 1 (inclusive).

Here is the minimal code of my problem

## TETRAGONAL
t11 = 361.0
t12 = 100.0
t13 = 62.0
t33 = 264.0
t44 = 59.0
t66 = 64.0

## MONOCLINIC
m11 = 327
m22 = 408
m33 = 258
m44 = 100
m55 = 81
m66 = 126
m12 = 142
m13 = 55
m16 = -21
m23 = 196
m26 = 31
m36 = -18
m45 = -23

Compliance_tetragonal = [
    [t11, t12, t13, 0.0, 0.0, 0.0],
    [t12, 0.0, 0.0, 0.0, 0.0, 0.0],
    [t13, 0.0, t33, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, t44, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, t66]
    ]

Compliance_tetragonal = as_tensor(Compliance_tetragonal)


Compliance_monoclinic = [
    [m11, m12, m13, 0.0, 0.0, m16],
    [m12, m22, 0.0, 0.0, 0.0, m26],
    [m13, 0.0, m33, 0.0, 0.0, m36],
    [0.0, 0.0, 0.0, m44, m45, 0.0],
    [0.0, 0.0, 0.0, m45, m55, 0.0],
    [m16, m26, m36, 0.0, 0.0, m66]
    ]

Compliance_monoclinic = as_tensor(Compliance_monoclinic)

N = 20
R = 1.0

sphr = Sphere(Point(0,0,0), R)

cb1 = Box(Point(0, 2*R, 2*R), Point(2*R, -2*R, -2*R))
cb2 = Box(Point(0, 2*R, -2*R), Point(-2*R, -2*R, 0))
cb3 = Box(Point(2*R, 0, 0), Point(-2*R, -2*R, 2*R))

domain = sphr-cb1-cb2-cb3

# nucleation volume
tol = 1E-14
nucleus = '(x[0] - 0.5)*(x[0] - 0.5) + (x[1] - 0.5)*(x[1] - 0.5) + (x[2] - 0.5)*(x[2] - 0.5) <= 0.1 + tol'

mesh = generate_mesh(domain, N)

X = SpatialCoordinate(mesh)
h = CellDiameter(mesh)

subdomain_0 = CompiledSubDomain(nucleus, tol=tol)

P1 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P2 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)

element = MixedElement([P1, P2])#, Re])

W = FunctionSpace(mesh, element)

w = Function(W)
v, p0 = TestFunctions(W)

u, c0 = split(w)

bdry_disp = Constant(0.1)

bcs = [DirichletBC(W.sub(0).sub(0), Constant(0), "x[0]<DOLFIN_EPS"),
       DirichletBC(W.sub(0).sub(1), Constant(0), "x[1]<DOLFIN_EPS"),
       DirichletBC(W.sub(0).sub(2), bdry_disp, "x[2]<DOLFIN_EPS"),
       DirichletBC(W.sub(1), Constant(1.0), subdomain_0)
       ]

'''
VOIGT FORMAT DEFINITIONS:
'''

def strain3voigt(e):
    return as_vector(
        [e[0, 0], 
         e[1, 1], 
         e[2, 2], 
         e[1, 2] * 2, 
         e[0, 2] * 2, 
         e[0, 1] * 2]
                    )


def voigt3stress(s):
    return as_tensor(
        [[s[0], s[5], s[4]], 
         [s[5], s[1], s[3]], 
         [s[4], s[3], s[2]]]
                     )

'''
CONSTITUTIVE LAW: FUNCTIONS AND DEFINITIONS
'''

# 3*3 matrix
def epsilon(v):
    return sym(grad(v))


def trans_strain(c0, e_v1):
    return (1-c0)*e_v1


def epsilon_elastic(v, c0):
    return epsilon(v)- trans_strain(c0, e_v1)


def Vstrain_elastic(v, c0):
    return strain3voigt(epsilon_elastic(v, c0))
  
  
def mixture_compliance(c0):
    global Compliance_monoclinic, Compliance_tetragonal
    return (c0 * Compliance_tetragonal) + ((1 - c0) * Compliance_monoclinic)


def sigma(v, c0):
    return voigt3stress(dot(mixture_compliance(c0), Vstrain_elastic(v, c0)))

'''
PENALTY FORCE:
''' 
x = X + u
R = Constant(R)

k = Constant(1.0)
G = Constant(1.0)

E = 9*k*G/(3*k+G)
C_pen = Constant(1e3) # (Dimensionless constant)
k_pen = C_pen*E/h # (Typical choice of contact penalty)
chi_pen = conditional(gt(x[2], R), 1, Constant(0))
t_pen = k_pen*chi_pen*as_vector([Constant(0), Constant(0), R-x[2]])


'''
RESIDUAL OF THE VARIATIONAL PROBLEM:
''' 


res = 0.5*inner(sigma(v, p0), epsilon_elastic(v, p0))*dx + 13.0*p0*(1 - p0)*dx + 18.0*(1 - p0)*dx - dot(t_pen, v)*ds

dw = Function(W)

J = derivative(res, w, dw)


'''
SOLUTION:
'''

w = Function(W)
problem = NonlinearVariationalProblem(res, w, bcs, J)
solver  = NonlinearVariationalSolver(problem)

u_, c0_ = w.split()

I am getting the following error consistently and I am not able to figure out how to proceed from here. I see that there is some issue with the solver function as it raises arity mismatch. The exact displayed message is -

Generating mesh with CGAL 3D mesh generator

Output from spyder call 'get_cwd':
                              Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/spyder_kernels/py3compat.py", line 356, in compat_exec
    exec(code, globals, locals)

  File "/Users/dhariwal/drive backup/Lu group/Zirconia/elasticity/aniso-3D-without-loop.py", line 333, in <module>
    problem = NonlinearVariationalProblem(res, w, bcs, J)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/dolfin/fem/problem.py", line 90, in __init__
    F = Form(F, form_compiler_parameters=form_compiler_parameters)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/jitcompiler.py", line 65, in jit_generate
    code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/compiler.py", line 142, in compile_form
    return compile_ufl_objects(forms, "form", object_names,

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
    form_datas = tuple(_analyze_form(form, parameters)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/analysis.py", line 89, in <genexpr>
    form_datas = tuple(_analyze_form(form, parameters)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/analysis.py", line 169, in _analyze_form
    form_data = compute_form_data(form,

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/algorithms/compute_form_data.py", line 418, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress)

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/corealg/map_dag.py", line 86, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])

  File "/Users/dhariwal/opt/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/algorithms/check_arities.py", line 48, in sum
    raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))

ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_0',).
2

As p0 is a test function, it should never be non-linear with respect to that function.

3

Thanks @dokken. I understand that this term is introducing nonlinearity in the residual, but it is necessary for computation. I would grateful if you know of any improvements to solve this nonlinear variational form.

Edit 1: I tried replacing 1 - c0 or 1 - p0 with c1 or q1 so that the residual becomes linear in test functions. c1 and q1 would be first defined in the mixed function space in the beginning of the code. The rationale is that c1 and c0 are conjugate scalar fields. Everywhere in the domain the sum of c1 and c0 would be 1. It does not work. The same error again.

4

You need to write down your variational form carefully. A weak formulation cannot be non linear with respect to the test function.

1 Like
5

Thanks @dokken. I have written my variational form correctly and the syntax is working now. Since I have a nonlinear residual (defined as res, its Jacobian as Dres), using the suggestions from demo_elasticity.py was done. The problem is now attempted by PETScKrylovSolver where I have modified parameters to my need, however the convergence is still not achieved.

A big thanks to @kamensky (especially elastic-sphere-compression example) @nate and @bhaveshshrimali whose informative answers to various queries on this platform have helped me, given that I am learning the math of 3D FEA concurrently, for the first time.

I am reproducing the error here and seeking the root cause in the code which is causing it and why. Prima facie observation is that for some reason residual norm of the solver is increasing after a certain point, but it does decrease initially up to a point.

Solving linear system of size 362560 x 362560 (PETSc Krylov solver).
PETSc Krylov solver starting to solve 362560 x 362560 system.
  0 KSP preconditioned resid norm 3.886660974961e+02 true resid norm 1.146557191213e+02 ||r(i)||/||b|| 1.000000000000e+00
  1 KSP preconditioned resid norm 2.966810222835e+02 true resid norm 1.210937155419e+01 ||r(i)||/||b|| 1.056150678483e-01
  2 KSP preconditioned resid norm 9.983863312690e+01 true resid norm 2.371276866593e+00 ||r(i)||/||b|| 2.068171465642e-02
  3 KSP preconditioned resid norm 7.370166281704e+01 true resid norm 1.237032761849e+00 ||r(i)||/||b|| 1.078910647746e-02
  4 KSP preconditioned resid norm 5.185251085415e+01 true resid norm 7.878716958911e-01 ||r(i)||/||b|| 6.871630145703e-03
  5 KSP preconditioned resid norm 3.392352954122e+01 true resid norm 5.064728035698e-01 ||r(i)||/||b|| 4.417335719940e-03
  6 KSP preconditioned resid norm 2.589494739501e+01 true resid norm 3.391618128047e-01 ||r(i)||/||b|| 2.958088923989e-03
  7 KSP preconditioned resid norm 2.146935095373e+01 true resid norm 2.571871910014e-01 ||r(i)||/||b|| 2.243125706877e-03
  8 KSP preconditioned resid norm 1.712512542526e+01 true resid norm 2.156247366577e-01 ||r(i)||/||b|| 1.880627833572e-03
  9 KSP preconditioned resid norm 1.366800264355e+01 true resid norm 1.649115669326e-01 ||r(i)||/||b|| 1.438319590130e-03
 10 KSP preconditioned resid norm 1.153654112645e+01 true resid norm 1.293354567489e-01 ||r(i)||/||b|| 1.128033191368e-03
 11 KSP preconditioned resid norm 1.016611567463e+01 true resid norm 1.165052572962e-01 ||r(i)||/||b|| 1.016131233480e-03
 12 KSP preconditioned resid norm 8.810891732436e+00 true resid norm 1.042247183336e-01 ||r(i)||/||b|| 9.090232840746e-04
 13 KSP preconditioned resid norm 7.932552186807e+00 true resid norm 8.750884815915e-02 ||r(i)||/||b|| 7.632314273532e-04
 14 KSP preconditioned resid norm 7.475213967317e+00 true resid norm 8.018456388011e-02 ||r(i)||/||b|| 6.993507562871e-04
 15 KSP preconditioned resid norm 6.996974716345e+00 true resid norm 7.882007037301e-02 ||r(i)||/||b|| 6.874499674076e-04
 16 KSP preconditioned resid norm 6.646076221178e+00 true resid norm 7.273771832433e-02 ||r(i)||/||b|| 6.344011348216e-04
 17 KSP preconditioned resid norm 6.438538832646e+00 true resid norm 6.747066053507e-02 ||r(i)||/||b|| 5.884631054793e-04
 18 KSP preconditioned resid norm 6.300462731612e+00 true resid norm 6.761889543540e-02 ||r(i)||/||b|| 5.897559751369e-04
 19 KSP preconditioned resid norm 6.255534477999e+00 true resid norm 6.833561482641e-02 ||r(i)||/||b|| 5.960070317480e-04
 20 KSP preconditioned resid norm 6.346841647845e+00 true resid norm 6.656854720232e-02 ||r(i)||/||b|| 5.805950868609e-04
 21 KSP preconditioned resid norm 6.505075316143e+00 true resid norm 6.820472298630e-02 ||r(i)||/||b|| 5.948654241499e-04
 22 KSP preconditioned resid norm 6.733835353448e+00 true resid norm 7.220615289820e-02 ||r(i)||/||b|| 6.297649471966e-04
 23 KSP preconditioned resid norm 6.910707797664e+00 true resid norm 7.360800943830e-02 ||r(i)||/||b|| 6.419916075921e-04
 24 KSP preconditioned resid norm 7.138478699044e+00 true resid norm 7.432036538063e-02 ||r(i)||/||b|| 6.482046072383e-04
 25 KSP preconditioned resid norm 7.505562830334e+00 true resid norm 7.925796254416e-02 ||r(i)||/||b|| 6.912691591103e-04
 26 KSP preconditioned resid norm 7.877118582402e+00 true resid norm 8.502910733497e-02 ||r(i)||/||b|| 7.416037157730e-04
 27 KSP preconditioned resid norm 8.271484235744e+00 true resid norm 8.626598179733e-02 ||r(i)||/||b|| 7.523914415999e-04
 28 KSP preconditioned resid norm 8.664723653491e+00 true resid norm 8.992825726693e-02 ||r(i)||/||b|| 7.843329400063e-04
 29 KSP preconditioned resid norm 9.076958595368e+00 true resid norm 9.686283838074e-02 ||r(i)||/||b|| 8.448147124547e-04
 30 KSP preconditioned resid norm 9.591042860646e+00 true resid norm 1.012577008742e-01 ||r(i)||/||b|| 8.831456612038e-04
 31 KSP preconditioned resid norm 1.011541324100e+01 true resid norm 1.045817270660e-01 ||r(i)||/||b|| 9.121370296007e-04
 32 KSP preconditioned resid norm 1.065799135605e+01 true resid norm 1.121123240703e-01 ||r(i)||/||b|| 9.778171113446e-04
 33 KSP preconditioned resid norm 1.133828428452e+01 true resid norm 1.201465714236e-01 ||r(i)||/||b|| 1.047889912029e-03
 34 KSP preconditioned resid norm 1.224316816849e+01 true resid norm 1.275651490407e-01 ||r(i)||/||b|| 1.112592987235e-03
 35 KSP preconditioned resid norm 1.322056537205e+01 true resid norm 1.371923082984e-01 ||r(i)||/||b|| 1.196558787907e-03
 36 KSP preconditioned resid norm 1.446967252353e+01 true resid norm 1.520718435427e-01 ||r(i)||/||b|| 1.326334566720e-03
 37 KSP preconditioned resid norm 1.576779473671e+01 true resid norm 1.656446086172e-01 ||r(i)||/||b|| 1.444713005917e-03
 38 KSP preconditioned resid norm 1.744525513622e+01 true resid norm 1.801758332543e-01 ||r(i)||/||b|| 1.571450902189e-03
 39 KSP preconditioned resid norm 1.973013532154e+01 true resid norm 2.042661032256e-01 ||r(i)||/||b|| 1.781560525642e-03
 40 KSP preconditioned resid norm 2.248861373914e+01 true resid norm 2.366658172952e-01 ||r(i)||/||b|| 2.064143150547e-03
 41 KSP preconditioned resid norm 2.604850488776e+01 true resid norm 2.719379789703e-01 ||r(i)||/||b|| 2.371778582477e-03
 42 KSP preconditioned resid norm 3.138629661622e+01 true resid norm 3.257768136305e-01 ||r(i)||/||b|| 2.841348134460e-03
 43 KSP preconditioned resid norm 3.926255297759e+01 true resid norm 4.114872866725e-01 ||r(i)||/||b|| 3.588894560395e-03
 44 KSP preconditioned resid norm 5.228432521344e+01 true resid norm 5.529776416210e-01 ||r(i)||/||b|| 4.822939892219e-03
 45 KSP preconditioned resid norm 7.891588393042e+01 true resid norm 8.178337766158e-01 ||r(i)||/||b|| 7.132952310479e-03
 46 KSP preconditioned resid norm 1.597900707279e+02 true resid norm 1.649405429891e+00 ||r(i)||/||b|| 1.438572312425e-02
*** Warning: Krylov solver did not converge in 47 iterations (PETSc reason DIVERGED_INDEFINITE_MAT, residual norm ||r|| = 1.597901e+02).
PETSc Krylov solver (cg, gamg) failed to converge in 47 iterations.
Elapsed wall, usr, sys time: 28.3261, 27.13, 1.23 (PETSc Krylov solver)

Apologies for not being able to shrink the code, however I tried to make at decently commented and spacious.

from dolfin import *
from mshr import *
from matplotlib import pyplot as plt
from vedo.dolfin import *

parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
              "eliminate_zeros": True, \
              "precompute_basis_const": True, \
              "precompute_ip_const": True}

# Test for PETSc
if not has_linear_algebra_backend("PETSc"):
   print("DOLFIN has not been configured with PETSc. Exiting.")
   exit()

# Set backend to PETSC
parameters["linear_algebra_backend"] = "PETSc"



'''
MATERIAL PROPERTIES:
'''
# Define elastic constants

# Create fourth-order stiffness tensors for each phase (in Voigt notation)

Compliance_tetragonal = [
   [169.0, 138.0, 138.0, 0.0, 0.0, 0.0],
   [138.0, 169.0, 138.0, 0.0, 0.0, 0.0],
   [138.0, 138.0, 169.0, 0.0, 0.0, 0.0],
   [0.0, 0.0, 0.0, 40.0, 0.0, 0.0],
   [0.0, 0.0, 0.0, 0.0, 40.0, 0.0],
   [0.0, 0.0, 0.0, 0.0, 0.0, 40.0]
   ]

Compliance_tetragonal = as_tensor(Compliance_tetragonal)


Compliance_monoclinic = [
   [223.0, 129.0, 99.0, 0.0, 0.0, 27.0],
   [129.0, 241.0, 125.0, 0.0, 0.0, -9.0],
   [99.0, 125.0, 200.0, 0.0, 0.0, 4.0],
   [0.0, 0.0, 0.0, 76.0, -4.0, 0.0],
   [0.0, 0.0, 0.0, -4.0, 45.0, 0.0],
   [27.0, -9.0, 4.0, 0.0, 0.0, 90.0]
   ]

Compliance_monoclinic = as_tensor(Compliance_monoclinic)




alpha = 1.0243
beta = -0.0247
gamma = 0.9563
delta = 0.058


e_v1 = as_tensor([[gamma - 1, beta, beta],
       [beta, alpha - 1, delta],
       [beta, delta, alpha - 1]
       ])


e_v2 = as_tensor([[gamma - 1, - beta, - beta],
       [- beta, alpha - 1, delta],
       [- beta, delta, alpha - 1]
       ])


e_v3 = as_tensor([[gamma - 1, - beta, beta],
       [- beta, alpha - 1, - delta],
       [beta, - delta, alpha - 1]
       ])


e_v4 = as_tensor([[gamma - 1, beta, - beta],
       [beta, alpha - 1, - delta],
       [- beta, - delta, alpha - 1]
       ])


e_v5 = as_tensor([[alpha - 1, beta, delta],
       [beta, gamma - 1, beta],
       [delta, beta, alpha - 1]
       ])

e_v6 = as_tensor([[alpha - 1, - beta, delta],
       [- beta, gamma - 1, - beta],
       [delta, - beta, alpha - 1]
       ])

e_v7 = as_tensor([[alpha - 1, - beta, - delta],
       [- beta, gamma - 1, beta],
       [- delta, beta, alpha - 1]
       ])

e_v8 = as_tensor([[alpha - 1, beta, - delta],
       [beta, gamma - 1, - beta],
       [- delta, - beta, alpha - 1]
       ])


e_v9 = as_tensor([[alpha - 1, delta, beta],
       [delta, alpha - 1, beta],
       [beta, beta, gamma - 1]
       ])


e_v10 = as_tensor([[alpha - 1, delta, - beta],
       [delta, alpha - 1, - beta],
       [- beta, - beta, gamma - 1]
       ])


e_v11 = as_tensor([[alpha - 1, - delta, beta],
       [- delta, alpha - 1, - beta],
       [beta, - beta, gamma - 1]
       ])


e_v12 = as_tensor([[alpha - 1, - delta, - beta],
       [- delta, alpha - 1, beta],
       [- beta, beta, gamma - 1]
       ])


var_list = [e_v1, e_v2, e_v3, e_v4, e_v5, e_v6, e_v7, e_v8, e_v9, e_v10, e_v11, e_v12]




'''
DECLARE GEOMETRY:
'''
# Create mesh of quarter circle 
# (to takeadvantage of symmetry).

N = 25
R = 1.0

center = Point(0.0, 0.0, 0.0)

sphr = Sphere(center, R)

cb1 = Box(Point(-2*R, -2*R, 0), Point(2*R, 2*R, -2*R))

#cb1 = Box(Point(0, 2*R, 2*R), Point(2*R, -2*R, -2*R))
#cb2 = Box(Point(0, 2*R, -2*R), Point(-2*R, -2*R, 0))
#cb3 = Box(Point(2*R, 0, 0), Point(-2*R, -2*R, 2*R))


# martensite nucleus
nucleus = Point(-0.123, 0.124, 0.125)


domain = sphr-cb1 #-cb2-cb3

mesh = generate_mesh(domain, N)

#plot(mesh)


X = SpatialCoordinate(mesh)
h = CellDiameter(mesh)



print(mesh.num_cells())


markers = MeshFunction('size_t', mesh, mesh.topology().dim(), mesh.domains())
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1, mesh.domains())

#markers = MeshFunction('size_t', mesh, 0)
#boundaries = FacetFunction('size_t', mesh, 2)

#markers.set_all(0)
#boundaries.set_all(0)

for cell in cells(mesh):
   if cell.collides(nucleus):
       print('Hooray, we found the nucleating cell')
       markers[cell] = 3
       for f in facets(cell):
           boundaries[f] = 5
           
           
'''
WEAK FORMULATION:
'''

P1 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P2 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)

martensite_element = MixedElement([P2, P2])

element = MixedElement([P1, martensite_element])#, P3])#, Re])

W = FunctionSpace(mesh, element)


w = Function(W)

u, c1 = split(w) # c1 is volume fraction of martensite; c0 is volume fraction of austenite

c11, c12 = split(c1)


v, p1 = TestFunctions(W)

p11, p12 = split(p1)


'''
BOUNDARY CONDITIONS:
'''
# Quarter-symmetry BCs, with a constant vertical displacement
# applied to the bottom.

bdry_disp = Constant(0.22)

bc1 = DirichletBC(W.sub(0).sub(2), bdry_disp, "x[2]<DOLFIN_EPS")

bc2 = DirichletBC(W.sub(1).sub(0), Constant(0.99), boundaries, 5)

#bc3 = DirichletBC(W.sub(0).sub(0), Constant(0.0), "near(x[0], 0) and near(x[1], 0)")

#bc4 = DirichletBC(W.sub(0).sub(1), Constant(0.0), "near(x[0], 0) and near(x[1], 0)")



bcs = [bc1, bc2]#, bc3, bc4]


'''
VOIGT FORMAT DEFINITIONS:
'''

def strain3voigt(e):
   return as_vector(
       [e[0, 0], 
        e[1, 1], 
        e[2, 2], 
        e[1, 2] * 2, 
        e[0, 2] * 2, 
        e[0, 1] * 2]
                   )


def voigt3stress(s):
   return as_tensor(
       [[s[0], s[5], s[4]], 
        [s[5], s[1], s[3]], 
        [s[4], s[3], s[2]]]
                    )


'''
CONSTITUTIVE LAW: FUNCTIONS AND DEFINITIONS
'''

# 3*3 matrix
def epsilon(v):
   return sym(grad(v))


def c(c1):
   c11, c12 = split(c1)
   return c11 + c12


def c0(c1):
   return 1.0 - c(c1)


def trans_strain(c1):
   global e_v1, e_v2
   c11, c12 = split(c1)
   return c11*e_v1 + c12*e_v2


def epsilon_elastic(v, c1):
   return epsilon(v)- trans_strain(c1)


def Vstrain_elastic(v, c1):
   return strain3voigt(epsilon_elastic(v, c1))


def mixture_compliance(c1):
   global Compliance_monoclinic, Compliance_tetragonal
   return (c0(c1) * Compliance_tetragonal) + (c(c1) * Compliance_monoclinic)


def sigma(w):
   u, c1 = split(w)
   return voigt3stress(dot(mixture_compliance(c1), Vstrain_elastic(u, c1)))


'''
PENALTY FORCE:
''' 
x = X + u
R = Constant(R)

#k = Constant(210.0) # GPa
k = Constant(1.0)

#G = Constant(86.4) # GPa
G = Constant(1.0)

E = 9*k*G/(3*k+G)

C_pen = Constant(1e3) # (Dimensionless constant)

k_pen = C_pen*E/h # (Typical choice of contact penalty)

chi_pen = conditional(gt(x[2], R), 1, Constant(0))

t_pen = k_pen*chi_pen*as_vector([Constant(0), Constant(0), R-x[2]])



'''
RESIDUAL OF THE VARIATIONAL PROBLEM:
''' 

res = 0.5*inner(sigma(w), epsilon_elastic(v, p1))*dx + 13.0*(1e-3)*c(p1)*c0(c1)*dx + 18.0*(1e-3)*c(p1)*dx - dot(t_pen, v)*ds

#a, l = 0.5*inner(sigma(v, p0), epsilon_elastic(v, p0))*dx + 13.0*p0*d0(p0)*dx + 18.0*d0(p0)*dx, dot(t_pen, v)*ds

#dw = Function(W)

Dres = derivative(res, w)

#W2 = FunctionSpace(mesh, P2)

#ci_min, ci_max = Function(W2), Function(W2)
#ci_min.x.array[:] = 0.0
#ci_max.x.array[:] = 1.5


A, b = assemble_system(Dres, res, bcs=bcs)
#apply_lifting(r, Dres, bcs, x0=[w.vector], scale=-1.0)

set_log_level(1)


# Create PETSC smoothed aggregation AMG preconditioner and attach near
# null space
pc = PETScPreconditioner("petsc_amg")

#PETScOptions.set("KSPConvergedSkip") #this has no effect

# Use Chebyshev smoothing for multigrid
PETScOptions.set("mg_levels_ksp_type", "chebyshev")
PETScOptions.set("mg_levels_pc_type", "jacobi")


# Improve estimate of eigenvalues for Chebyshev smoothing
PETScOptions.set("mg_levels_esteig_ksp_type", "cg")
PETScOptions.set("mg_levels_ksp_chebyshev_esteig_steps", 50)


# Create CG Krylov solver and turn convergence monitoring on
solver = PETScKrylovSolver("cg", pc)
solver.parameters["monitor_convergence"] = True
solver.parameters["error_on_nonconvergence"] = False
solver.parameters["nonzero_initial_guess"] = True
solver.parameters["report"] = True

# Set matrix operator
solver.set_operator(A)


# Compute solution
solver.solve(w.vector(), b)
6

I am facing the same kind of error with my code.
Attaching the link to my code.
Could anyone help me correct it?