Solving a PDE with Neumann BC between subdomains

Hello,

I am trying to model a pressure-temperature PDE system, where the pressure variable is defined only on part of the domain. In order to understand how to solve it with FEniCS, I wrote a simplified test problem that should be similar to my system. The domains are:

$$\Omega_2 = \left[\tfrac{1}{4},\tfrac{3}{4}\right]\times \left[\tfrac{1}{4},\tfrac{3}{4}\right], \Omega_1=[0,1]\times[0,1]-\Omega_2$$

Domain and boundaries: http://imgur.com/DXf6qNM

The equations:

$$\Omega_1:\left\{\begin{array}{l}
\Delta P + P - T = f(x)\\
\Delta T + T + P = g(x)
\end{array}\right.$$

$$\Omega_2:\{\Delta T + 2T = 0$$

with Dirichlet boundary conditions on the edges of the unit square, and Neumann boundary condition for $P$ between $\Omega_1$ and $\Omega_2$ ($\nabla P\cdot n=0$ on $\Gamma_2$, $\nabla P\cdot n=1$ on $\Gamma_3$, and $\nabla P\cdot n=-1$ on $\Gamma_4$). The solutions should be

$$P(x,y)=\left(x-\tfrac{1}{4}\right)^2\left(x-\tfrac{3}{4}\right)^2+y$$

$$T(x,y)=\sin(x+y)$$

This is my attempt at the solution:

from dolfin import *
import numpy as np
import scipy.io as IO

mesh = RectangleMesh(0, 0, 1, 1, 20, 20, "crossed")

# Boundaries
class OuterBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

# East and West boundaries
class Internal_WE_Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return  between(x[1], (0.25,0.75)) and \
            (near(x[0], 0.75) or near(x[0], 0.25))

# North boundary
class Internal_N_Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.25,0.75)) and near(x[1], 0.75)

# South boundary
class Internal_S_Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.25,0.75)) and near(x[1], 0.25)

# Create subdomains
class Omega2Domain(SubDomain):
    def inside(self, x, on_boundary):
        return  between(x[0],(0.25,0.75)) and \
            between(x[1],(0.25,0.75))

# Mark subdomains
subdomains = CellFunction("size_t", mesh)
subdomains.set_all(0)       # Omega_0
Omega2Domain().mark(subdomains, 1)  # Omega_1

# Mark boundaries
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
OuterBoundary().mark(boundaries, 1)
Internal_WE_Boundary().mark(boundaries, 2)
Internal_N_Boundary().mark(boundaries, 3)
Internal_S_Boundary().mark(boundaries, 4)

# Initialize spaces & functions
V = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V, V])
[T, P] = TrialFunction(W)
[v1, v2] = TestFunction(W)

P_Exact = Expression("pow((x[0]-0.25)*(x[0]-0.75),2) + x[1]")
P_Lapl_Exact = Expression("12*x[0]*x[0] - 12*x[0]+2.75")
T_Exact = Expression("sin(x[0]+x[1])")
T_Lapl_Exact = Expression("-2*sin(x[0]+x[1])")

bc_T = DirichletBC(W.sub(0), T_Exact, boundaries, 1)
bc_P = DirichletBC(W.sub(1), P_Exact, boundaries, 1)

dx = Measure("dx")[subdomains]
ds = Measure("ds")[boundaries]

a_form = (-inner(grad(P),grad(v2)) + inner(P,v2) - inner(T,v2)) * dx(0) \
    +(-inner(grad(T),grad(v1)) + inner(T,v1) + inner(P,v1)) * dx(0) \
    +(-inner(grad(T),grad(v1)) + 2*inner(T,v1)) * dx(1)

b_form = (P_Lapl_Exact + P_Exact - T_Exact) * v2 * dx(0) \
    + (T_Lapl_Exact + T_Exact + P_Exact) * v1 * dx(0)
    + Constant(-1.0) * v2 * ds(3) + Constant(1.0) * v2 * ds(4)

u = Function(W)
solve(a_form == b_form, u, [bc_T, bc_P])

(T, P) = u.split(deepcopy = True)
plot(P)
interactive()

However, the line (in definition of b_form)

    + Constant(-1.0) * v2 * ds(3) + Constant(1.0) * v2 * ds(4)

seems to have no effect on the solution, and thus I am unable to impose any Neumann BC's between $\Omega_1$ and $\Omega_2$. I would appreciate any help in dealing with this issue.

Thanks in advance,

safinenko

EDIT: a short version of the question: Let's say I have mesh defined over union of two subdomains, $\Omega=\Omega_1\cup\Omega_2$. Let $\Gamma$ be the boundary between the two subdomains. If I want to solve a PDE only on $\Omega_1$, is there a way to impose a Neumann BC over $\Gamma$?