dbc

    from dolfin import *

    parameters["linear_algebra_backend"] = "uBLAS"

    Lx = Ly = 1.0
    R = Constant(2500)
    P = Constant(2.026)
    t = 0.0
    T = 20.0
    dt = 0.01

    mesh = UnitSquareMesh(64, 64)


    class onBase(SubDomain):

        def inside(self, x, on_boundary):
            return bool(x[1] < DOLFIN_EPS and on_boundary)

    class onLid(SubDomain):

        def inside(self, x, on_boundary):
            return bool(x[1] > (Ly - DOLFIN_EPS) and on_boundary)

    class onSides(SubDomain):

        def inside(self, x, on_boundary):
            return bool(not(near(x[1], 0) or near(x[1], Ly)) and on_boundary)

    class onFaces(SubDomain):

        def inside(self, x, on_boundary):
            return bool(on_boundary)

    base = onBase()
    lid = onLid()
    sides = onSides()
    faces = onFaces()

    # Create the function spaces for u, p, c
    V = VectorFunctionSpace(mesh, "CG", 1, dim=2)
    W = FunctionSpace(mesh, "CG", 1)
    X = FunctionSpace(mesh, "CG", 1)
    VWX = MixedFunctionSpace([V, W, X])  # mixed functions space

    # Define trail and test functions
    du_ = TrialFunction(VWX)
    v_u, v_p, v_c = TestFunctions(VWX)

    # define functions
    u_ = Function(VWX)  # current solution
    u0_ = Function(VWX)  # solution from previous step

    # Split mixed functions
    du, dp, dc = split(du_)
    u, p, c = split(u_)
    u0, p0, c0 = split(u0_)
    I = Identity(2)

    #define the boundary conditions
    no_slip = [DirichletBC(VWX.sub(0), Constant((0.0, 0.0)), lid),
               DirichletBC(VWX.sub(0), Constant((0.0, 0.0)), base),
               DirichletBC(VWX.sub(0), Constant((0.0, 0.0)), sides)]

    temp = [DirichletBC(VWX.sub(2), Constant(0), lid),
            DirichletBC(VWX.sub(2), Constant(0), base)]

    bcs = no_slip + temp

    # Class representing the initial conditions
    class InitialConditions(Expression):
        def __init__(self):
            pass

        def eval(self, values, x):
            values[0] = 0
            values[1] = 0
            values[2] = 0
            values[3] = 0.01

        def value_shape(self):
            return (4,)

    # Create initial conditions and interpolate
    u_init_ = InitialConditions()
    u_.interpolate(u_init_)
    u0_.interpolate(u_init_)

    # Crank-Nicholson time stepping
    U = 0.5*(u + u0)
    C = 0.5*(c + c0)

    L = (1. /dt)*inner(u - u0, v_u)*dx - \
                     inner(p, div(v_u))*dx - \
                     R*P*inner(c*I[:, 1], v_u)*dx +  \
                     P*inner(grad(U), grad(v_u))*dx + \
                     inner(div(U), v_p)*dx + (1./dt)*inner(c - c0, v_c)*dx - \
                     inner(dot(U,I[:, 1]), v_c)*dx + \
                     inner(grad(C), grad(v_c))*dx + \
                     inner(0.5*(dot(u0, grad(u0)) + dot(u, grad(u))), v_u)*dx + \
                     inner(0.5*(dot(u0, grad(c0)) + dot(u, grad(c))), v_c)*dx

    J = derivative(L, u_, du_)
    problem = NonlinearVariationalProblem(L, u_, bcs, J)
    solver = NonlinearVariationalSolver(problem)

    while (t < T):
        t += dt
        u0_.vector()[:] = u_.vector()
        solver.solve()