Problem with Mixed-dimensional branch - Error: Unable to evaluate function at point

2

Please follow the guidelines on how to make a minimal working example.
There are many things that can be removed from your code, with still giving this error,

from fenics import *
from mshr import *
import logging
logging.getLogger('FFC').setLevel(logging.WARNING)

#Creation of mesh
background = Rectangle(Point(0,0), Point(0.3,0.07))
steel = Rectangle(Point(0,0), Point(0.3,0.01))
anode = Rectangle(Point(0.025,0.025), Point(0.125,0.035))
cathode = Rectangle(Point(0.175,0.025), Point(0.275,0.035))
electrolyte = background - steel - anode - cathode
domain = steel + electrolyte
domain.set_subdomain(1, electrolyte)
mesh = generate_mesh(domain, 5)

#Define subdomains
class Gamma_a(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and between(x[0], [0.025, 0.125]) and between(x[1], [0.025, 0.035])

class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 0.01

#Initializate subdomains instances
gamma_a = Gamma_a()
omega_1 = Omega_1()

#Creation of facet and cell markers and marking of mesh elements
facet_marker = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
gamma_a.mark(facet_marker, 1)

cell_marker = MeshFunction("size_t", mesh, mesh.topology().dim())
omega_1.mark(cell_marker, 1)

xdmffile = XDMFFile('facet_function.xdmf')
xdmffile.write(facet_marker)

xdmffile = XDMFFile('mesh_function.xdmf')
xdmffile.write(cell_marker)


#Definition of new meshes
mesh_1 = MeshView.create(cell_marker, 1)
mesh_a = MeshView.create(facet_marker, 1)

#Definition of measures
dx1 = Measure("dx", domain=mesh_1)
dxa = Measure("dx", domain=mesh_a)

#Definition of function space, trial and test functions
V1 = FunctionSpace(mesh_1, "CG", 2)
V2 = FunctionSpace(mesh_a, "CG", 2)
W = MixedFunctionSpace(V1, V2)
u = Function(W)
(u1, ia) = u.split()
(v1, ea) = TestFunctions(W)

xdmffile = XDMFFile("u1.xdmf").write(u1)
xdmffile = XDMFFile("ia.xdmf").write(ia)

Res = inner(u1,v1)*dx1 + inner(u1,ea)*dxa

#Solution of the system
solve(Res == 0, u)
(u0, u1) = u.split()