problem with discontinuous flux at an internal interface

Hello,
I am working with a problem related to different materials and discontinuous flux at the interface.
Here, I am considering the Laplace equation over omega=[0,1]x[0,1]
-div(D(x,y)*u(x,y))=0
Let u=u1 and D(x,y)=D1 on [0,1]x[0,1/2], u=u2 and D(x,y)=D2 on [0,1]x[1/2,1],
where u1 = x - 2*D2*y+1/2*D1+D2, u2 = x+D1*y
We can compute the flux at the interface y=1/2
D1*dot(grad(u1)*n1)=-2*D1*D2
D2*dot(grad(u2)*n2)=-D1*D2
Now I want to recover this solution by Fenics. To do this, I applied DirichletBC for extenal boundaries. Since the flux is discontinuous, I added kcoeff*dot(grad(u),n)*v*ds(0) to BilinearForm. However, this term did not effect on the results and the results were not correct. Can you please show me the way to fix it?
Thanks a lot,
Best regards,
Dang

My code is following:

ufl file
Poisson.ufl
# Compile this form with FFC: ffc -l dolfin Poisson.ufl.
element = FiniteElement("Lagrange", triangle, 1)
u = TrialFunction(element)
v = TestFunction(element)
kcoeff = Coefficient(element);
n = FacetNormal("triangle");
f = Coefficient(element)
a = kcoeff*inner(grad(u), grad(v))*dx - kcoeff*dot(grad(u),n)*v*ds(0)
L = f*v*dx

FunctionSpace1Dtri.ufl
# Compile this form with FFC: ffc -l dolfin FunctionSpace1Dtri.ufl
element = FiniteElement("DG", "triangle", 0)

main.cpp

#include <dolfin.h>
#include "Poisson.h"
#include "FunctionSpace1Dtri.h"

using namespace dolfin;

double Darray[2] = {3e-3,3e-3};

void ExactSol(Mesh mesh, FunctionSpace V, double*Darray,Function& exactsol)
{
  for (VertexIterator v(mesh);!v.end();++v)
  {
    Point p = v->point();
    double val = 0.0;
    if (p.y()<=1.0/2.0) // u1
    {
      val = p.x()-2.0*Darray[1]*p.y()+1.0/2.0*Darray[0]+Darray[1];
    }
    else // u2
    {
      val = p.x()+Darray[0]*p.y();
    }
    exactsol.vector()->setitem(v->index(),val);
  }
}

class Omega0 : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[1] <= 0.5;
  }
};

class Interface : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return fabs(x[1] - 0.5)<1e-10;
  }
};

int main()
{
  // Create mesh
  UnitSquare mesh(8, 8);

  Poisson::FunctionSpace V(mesh);
  Function u0(V);
  ExactSol(mesh, V, Darray,u0);
  plot(u0);

  class DirichletBoundary : public SubDomain
  {
    bool inside(const Array<double>& x, bool on_boundary) const
    {
      return on_boundary;
    }
  };

  MeshFunction<unsigned int>sub_domains(mesh, mesh.topology().dim(),1);
  MeshFunction<unsigned int> interface_marker(mesh, mesh.topology().dim()-1,1);

  // Mark Omega0 as sub domain 0
  Omega0 omega0;
  omega0.mark(sub_domains, 0);

  Interface interface;
  interface.mark(interface_marker, 0);

  Function interface_fun(V); // for checking the interface
  interface_fun.vector()->zero();
  for (FacetIterator fc(mesh);!fc.end();++fc)
  {
     if (interface_marker.values()[fc->index()]==0)
     {
       for (VertexIterator v(*fc);!v.end();++v)
       {
          interface_fun.vector()->setitem(v->index(),1);
       }
    }
  }
  plot(interface_fun);

  FunctionSpace1Dtri::FunctionSpace V1(mesh);
  Function kcoeff(V1);

  // set diffusion coefficients for two subdomains
  for (CellIterator c(mesh);!c.end();++c)
  {
    kcoeff.vector()->setitem(c->index(),Darray[sub_domains.values()[c->index()]]);

  }
  plot(kcoeff);

  DirichletBoundary dirichlet_boundary;
  DirichletBC bc0(V, u0, dirichlet_boundary);
  // Compute solution
  Function u(V);

  Poisson::BilinearForm a(V, V);
  a.kcoeff = kcoeff;
  //a.ds = interface_marker;
  Constant f(0.0);
  Poisson::LinearForm L(V);
  L.f = f;

  solve(a == L, u, bc0);
  plot(u);

  Function error(V);
  for (VertexIterator v(mesh);!v.end();++v)
  {
    error.vector()->setitem(v->index(),u0.vector()->getitem(v->index())-u.vector()->getitem(v->index()));
  }
  plot(error);

  return 0;
}