Hi, I am a student trying to solve a Stokes system that has a liquid phase bubble within a solid phase, which I represent with two separate viscosities. I am following the example shown here at the bottom example with varying permeabilities.

My code is here:

```
# Set domain geometries
L = 400
H = 200
resolution = 32
domain = Rectangle(Point(0.0,0.0), Point(L,H))
liquid = Circle(Point(100, 0.5*H), 50)
solid = domain - liquid
domain.set_subdomain(1, liquid)
domain.set_subdomain(2, solid)
# Export mesh as xml file, read the xml file
domain_mesh = generate_mesh(domain, resolution)
File("meshes/twomaterials/mesh.xml") << domain_mesh
mesh = Mesh("meshes/twomaterials/mesh.xml")
# Define subdomains ----------------------------
# Sub domain for no-slip (mark whole boundary, inflow and outflow will overwrite)
class Noslip(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
# Sub domain for inflow (left)
class Inflow(SubDomain):
def inside(self, x, on_boundary):
return x[0] < DOLFIN_EPS and on_boundary
# Sub domain for outflow (right)
class Outflow(SubDomain):
def inside(self, x, on_boundary):
return x[0] > L - 1e3*DOLFIN_EPS and on_boundary
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(3)
noslip = Noslip()
noslip.mark(sub_domains, 0)
inflow = Inflow()
inflow.mark(sub_domains, 1)
outflow = Outflow()
outflow.mark(sub_domains, 2)
file = File("meshes/twomaterials/subdomains.xml")
file << sub_domains
# defining problem space -----------------
# Taylor-Hood function space
VE = VectorElement("DG", mesh.ufl_cell(), 2) # previously was CG, not DG
QE = FiniteElement("DG", mesh.ufl_cell(), 1)
WE = VE * QE
W = FunctionSpace(mesh, WE)
V = FunctionSpace(mesh, VE)
Q = FunctionSpace(mesh, QE)
# boundary conditions-----------------------
# No-slip boundary condition for the liquid pocket
noslip = Constant((0, 0))
bc0 = DirichletBC(W.sub(0), noslip, sub_domains, 0) # subdomain marked 0 (noslip)
# Inflow boundary condition for velocity
YMAX = H
inflow = Expression(("4*(x[1]*(YMAX-x[1]))/(YMAX*YMAX)", "0."), YMAX=YMAX, element = V.ufl_element())
bc1 = DirichletBC(W.sub(0), inflow, sub_domains, 1) # subdomain marked 1 (inflow)
# Boundary condition for pressure at outflow
# x0 = 0
zero = Constant(0)
bc2 = DirichletBC(W.sub(1), zero, sub_domains, 2) # subdomain marked 2 (outflow)
# Collect boundary conditions
bcs = [bc0, bc1, bc2]
# "Define subdomain markers and integration measure"
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)
# Define an Expression for viscosity
class Viscosity(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index]==1: # liquid
values[0]==1 # viscosity in liquid is 1
else:
values[0]==3e15 # solid viscosity much higher
def value_shape(self): # from https://fenicsproject.discourse.group/t/problems-interpolating-a-userexpression-and-plotting-it/1303/2
#return (1,)
return ()
eta = Viscosity(markers)
# Define variational problem
I = Identity(2)
def epsilon(u): # strain rate tensor
return 0.5*(grad(u) + grad(u).T)
def sigma(u,p): # Cauchy stress tensor
return 2*eta*epsilon(u) - p*I # from https://fenicsproject.org/pub/tutorial/html/._ftut1008.html
# Define trial and test functions
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
f = Constant((0, 0)) # set f, body force, to 0
a = inner(sigma(u,p), grad(v))*dx + q*div(u)*dx
L = inner(f, v)*dx
# Compute solution
w = Function(W)
solve(a == L, w, bcs)
(u, p) = w.split(True)
plot(u)
plt.show()
plot(p)
plt.show()
```

The error I keep getting is shown below. Since there are no values being computed, my plots also come out blank. I have spent many hours trying to figure out where it is coming from but I am not even sure what part of the code is causing this malfunction Any help would be greatly appreciated!

```
/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:715: RuntimeWarning: divide by zero encountered in double_scalars
length = a * (widthu_per_lenu / (self.scale * self.width))
/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:715: RuntimeWarning: invalid value encountered in multiply
length = a * (widthu_per_lenu / (self.scale * self.width))
/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:767: RuntimeWarning: invalid value encountered in less
short = np.repeat(length < minsh, 8, axis=1)
/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:780: RuntimeWarning: invalid value encountered in less
tooshort = length < self.minlength
```