Thank you very much. This seems to have resolved my issue. I also had to create the function p using `conditional`

as opposed to `expression`

.

```
#p_Ex = Expression('x[0]-x[1] >= 0 ? -2*(exp(x[0]-x[1])-1)/(exp(x[0]-x[1])-1): -(x[0]-x[1])',
# degree = 3)
x0, x1 = MeshCoordinates(mesh)
p_Ex = conditional(ge(0, x0-x1), -2*(exp(x0-x1)-1)/(exp(x0-x1)-1), -(x0-x1))
```

My functioning code in full now reads as follows:

```
from fenics import *
from dolfin import *
from mshr import *
import numpy as np
#import ufl #from https://fenicsproject.org/qa/10763/how-can-create-an-expression-is-dependent-on-function-value/
#ufl.algorithms.apply_derivatives.CONDITIONAL_WORKAROUND = True #from https://fenicsproject.org/qa/10763/how-can-create-an-expression-is-dependent-on-function-value/
mesh = UnitSquareMesh(16, 16)
#***************Defining variables and corresponding function spaces**********
W = FiniteElement('DG', mesh.ufl_cell(), 0) #Trial/test space of Discontinuous Lagrange elements
V = FiniteElement('RT', mesh.ufl_cell(), 1) #Trial/test space of Raviart-Thomas elements
element = MixedElement([W,V])
M = FunctionSpace(mesh, element) #Creating mixed finite element space
W_0 = M.sub(0).collapse() #Retrieving trial/test space of constant elements
V_0 = M.sub(1).collapse() #Retrieving trial/test space of vector elements
#**************************boundary condition*********************************
#p_Ex = Expression('x[0]-x[1] >= 0 ? -2*(exp(x[0]-x[1])-1)/(exp(x[0]-x[1])-1): -(x[0]-x[1])',
# degree = 3)
x0, x1 = MeshCoordinates(mesh)
p_Ex = conditional(ge(0, x0-x1), -2*(exp(x0-x1)-1)/(exp(x0-x1)-1), -(x0-x1))
def boundary(x, on_boundary):
return on_boundary
#bc = DirichletBC(M.sub(0), u_Ex, boundary)
#bc = DirichletBC(W_0, u_Ex, boundary)
#*******************Initial Condition******************************************
p = project(p_Ex, W_0) #Projection of u_Ex onto scalar component of mixed space.
bc = DirichletBC(M.sub(0), p, boundary)
#*****************Setting values of certain parameters************************
k = 2.0 # Linearisation constant
T = 0.14 #Final time in corresponding parabolic problem
N = 10 #Number of steps in parabolic problem
tau = T/N #Time step size in parabolic problem
#*****************Defining variational form************************************
#uq = Function(M)
uq = TrialFunction(M) #Defining trial functions of the mixed formulation
#u, q = uq.split(True)
u,q = split(uq)
w,v = TestFunctions(M) #Defining test functions of the mixed formulation
bn_ = conditional(le(0, p), pow(pi,2)/2 - pow(p,2)/2, pow(pi,2)/2) #The function b acting on the projection u_p of u_Ex
a1 = k*(u*w*dx) + tau*((div(q))*w*dx) #bilinear form in equation for u
L1 = (bn_ + k*p)*w*dx #linear form in equation for u
a2 = dot(q, v)*dx - u*(div(v))*dx #bilinear form in equation for q
L2 = 0 #linear form in equation for q
a = a1 + a2
L = L1 + L2
#*****************Solving mixed problem****************************************
uq = Function(M) #Solution
solve(a==L, uq, bc) #Solving mixed variational form
u, q = uq.split(True) #Splitting solution into functions from the respective components of the mixed space
print(u)
unodes = u.compute_vertex_values(mesh) #Storing vertex values of numerical solution
print(unodes) #Displaying vertex values of numerical solution
```

Thanks again Nate and David .