The problem I am tring to simulate involves multiple domains. It consists of 3 layers stacked on one another. One set of equations must be solved over all three layers, and another set of equations only applied to the middle layer. In that middle layer the equations are coupled, and must therefore be solved simultaneously.

In order to figure out how to do this, I have started with a much simpler problem. 2 diffusion equations that are not even coupled:

- domain 1, D2, is unit square
- domain 2, D2, is the left half of D1
- basic diffusion equation of species “a” over D1
- del(a) = 1
- a = a0 @ x=0
- 0 gradient on all other boundaries

- basic diffusion equation of species “b” over
**only**D2- del(b) = 1
- b = b0 @ x=0
- 0 gradient on all other boundaries

I have included my code below. It does not work.

Could someone please let me know what I am doing wrong with my current code, or if this is even the correct approach in the first place.

# Resulting garbage output

```
In[1]: a.vector()[:]
Out[1]: array([nan, nan, nan, ..., 1., 1., 1.])
In[2]: b.vector()[:]
Out[2]: array([nan, nan, nan, ..., 1., 1., 1.])
```

# Code

```
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
# Define Boundary equations
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
# Define mesh
resolution = 64
full_domain = Rectangle(Point(0, 0), Point(1, 1))
b_domain = Rectangle(Point(0, 0), Point(0.5, 1))
full_domain.set_subdomain(1, b_domain)
mesh = generate_mesh(full_domain, resolution)
# Initialize boundary instance(s)
left_boundary = LeftBoundary()
# Define domains
domains = MeshFunction("size_t", mesh, 2, mesh.domains())
# Define boundaries on main mesh
boundaries = MeshFunction("size_t", mesh, 1)
boundaries.set_all(0)
left_boundary.mark(boundaries, 1)
# Plot subdomains to make sure that's working properly
plot(domains)
plt.show()
# Generate dx and ds for all meshes and submeshes
dx = Measure('dx', domain=mesh, subdomain_data=domains)
# Define input Data
a0 = 1
b0 = 1
f_a = Constant(1)
f_b = Constant(1)
# Define function space and basis functions
degree: int = 1
fe = FiniteElement("CG", mesh.ufl_cell(), degree=degree)
me = MixedElement([fe, fe])
V = FunctionSpace(mesh, me)
# Define trial and test functions
a, b = TrialFunction(V)
v_a, v_b = TestFunction(V)
# Define Dirichlet Boundary conditions
bcs = [DirichletBC(V.sub(0), a0, boundaries, 1),
DirichletBC(V.sub(1), b0, boundaries, 1)]
# Define variational form
F = inner(grad(a), grad(v_a))*dx(0) + f_a*v_a*dx(0) + \
inner(grad(b), grad(v_b))*dx(1) + f_b*v_b*dx(1)
# Seperate LHS and RHS of F
A, L = lhs(F), rhs(F)
# Solve the problem
u = Function(V)
solve(A == L, u, bcs)
a, b = u.split()
#Graph the results
plt.figure()
im = plot(a)
ax = plt.gca()
cbar = ax.figure.colorbar(im, ax=ax)
cbar.ax.set_ylabel("", rotation=-90, va="bottom")
plt.title("a")
plt.show()
plt.figure()
im = plot(b)
ax = plt.gca()
cbar = ax.figure.colorbar(im, ax=ax)
cbar.ax.set_ylabel("", rotation=-90, va="bottom")
plt.title("b")
plt.show()
```