problems with plotting subdomains

I created a divided domain with the stokes-equation in the first subdomain and the mixed-poisson-equation (darcy) in the second subdomain.
But now i get the following error when i want to plot the solution:

UMFPACK problem related to call to solve
Warning: UMFPACK reports that the matrix being solved is singular
...
at the end:
assert vmax>=vmin, "empty range, please specify vmin and/or vmax"
Assertion error: empty range, please specify vmin and/or vmax

Can anyone help?
Thanks!

Here is the code:

#-*- coding: utf-8 -*-
 from dolfin import *
 import numpy as np

 # Define mesh
 mesh = UnitSquare(32,32)

 #Subdomain 1
 # Gitter übergeben
 subdomains = CellFunction("uint", mesh)

 # Klasse des Teilgebiets
 class Domain_1(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 0.5)) # Koordinatenangabe des Teilgebiets

 # Objekt der Klasse erstellen
 sub_domain1 = Domain_1()
 sub_domain1.mark(subdomains,0)

 # Definition Funktionenräume
 U = FunctionSpace(mesh, "CG", 2)
 V = FunctionSpace(mesh, "CG", 1)
 W = U*V

 # Definition Trial- und Testfunktion
 (u, p) = TrialFunctions(W)
 (v, q) = TestFunctions(W)

 # Randbedingungen
 p_in = 1
 p_out = 0
 noslip = DirichletBC(W.sub(0), (0),
                       "on_boundary && \
                       (x[1] <= DOLFIN_EPS | x[1] >= 0.5-DOLFIN_EPS)")

 inflow = DirichletBC(W.sub(1), p_in, "x[0] <= 0.0 + DOLFIN_EPS*1000")
 outflow = DirichletBC(W.sub(1), p_out, "x[0] >= 0.5 - DOLFIN_EPS*1000")

 bcp = [noslip,inflow, outflow]

 # Definition f
 f = Expression("0")

 # Variationsformulierung
 a = inner(grad(u), grad(v))*dx + div(v)*p*dx(0) + q*div(u)*dx(0)
 L = inner(f,v)*dx(0)

 # Lösung berechnen
 w = Function(W)
 solve(a == L, w, bcp)
 (u, p) = w.split()

 # Subdomain 2
 # Gitter übergeben
 subdomains = CellFunction("uint", mesh)

 # Klasse des Teilgebiets
 class Domain_2(SubDomain):
    def inside(self,x,on_boundary):
        return between(x[0], (0.5,1.0)) # Koordinatenangabe des Teilgebiets

 # Objekt der Klasse erstellen
 sub_domain2 = Domain_2()
 sub_domain2.mark(subdomains,1)

 # Define function spaces and mixed (product) space
 BDM = FunctionSpace(mesh, "BDM", 1)
 DG = FunctionSpace(mesh, "DG", 0)
 CG = FunctionSpace(mesh, "CG", 1)
 W = MixedFunctionSpace([BDM, DG, CG])

 # Define trial and test functions
 (sigma, u, p) = TrialFunctions(W)
 (tau, v, q) = TestFunctions(W)

 #Define pressure boundary condition
 p_in = 1
 p_out = 0
 noslip = DirichletBC(W.sub(1), (0),
                       "on_boundary && \
                       (x[1] <= 0.5 + DOLFIN_EPS | x[1] >= 1.0-DOLFIN_EPS)")

 inflow = DirichletBC(W.sub(2), p_in, "x[0] <= 0.5 + DOLFIN_EPS*1000")
 outflow = DirichletBC(W.sub(2), p_out, "x[0] >= 1.0 - DOLFIN_EPS*1000")
 bcp = [noslip,inflow, outflow]

 # Define f

 #f = Expression("0")

 f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")

 # Define variational form
 a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(1) + inner(p,q)*dx(1) + u*q*dx(1)
 L = f*v*dx(1)

 # Compute solution
 w = Function(W)
 solve(a == L, w, bcp)
 (sigma, u, p) = w.split()

 # plot
 plot(u, axes = True, interactive=True, title = "u")
 plot(p, axes = True, interactive=True, title = "p")
 plot(sigma, axes = True, interactive=True, title = "sigma")