Finite element method for elliptic interface problem

Why we can't use the codes to deal with the elliptic interface problem?

The model is described as :
-\nabla u1 = f in \Omega1
-\nabla u2 = f in \Omega2
u1=u2 on interface
\partial{u1}/\partial{n1}+\partial{u2}/\partial{n2}=0 [attachment=undefined]numerical solution.png[/attachment] [attachment=undefined]true solution.png[/attachment]

from ngsolve import *
from netgen.geom2d import SplineGeometry
from numpy import pi
import netgen.gui

def twoDomain():
geometry = SplineGeometry()

# point coordinates ...
pnts = [ (0,0), (0.5,0), (0.5,1), (0,1), (1,0), (1,1)]
pnums = [geometry.AppendPoint(*p) for p in pnts]
# start-point, end-point, boundary-condition, domain on left side, domain on right side:
lines = [ (pnums[0],pnums[1],1,1,0), (pnums[1],pnums[2],2,1,2),
(pnums[2],pnums[3],3,1,0), (pnums[3],pnums[0],4,1,0), (pnums[1], pnums[4],5,2,0),
(pnums[4],pnums[5],6,2,0), (pnums[5],pnums[2],7,2,0)]
for p1,p2,bc,left,right in lines:
geometry.Append( ["line", p1, p2], bc=bc, leftdomain=left, rightdomain=right)
return geometry

mesh = Mesh(twoDomain().GenerateMesh (maxh=0.1)) #interface 也是边界边，注意一下
#Draw(mesh)
#print ("boundary labels: ", mesh.GetBoundaries()) #interface 也是边界边，注意一下

order = 1

uex = sin(pi*x)*sin(pi*y)
source = 2*pi*pi*sin(pi*x)*sin(pi*y)


V1 = H1(mesh, order=order, dirichlet=[1,3,4], definedon=[1])
V2 = H1(mesh, order=order, dirichlet=[5,6,7], definedon=[2])
V = FESpace([V1,V2])

u1,u2 = V.TrialFunction()
v1,v2 = V.TestFunction()

n = specialcf.normal(2)
h = specialcf.mesh_size

a = BilinearForm(V)
a += SymbolicBFI (grad(u1)*grad(v1),definedon=[1])
a += SymbolicBFI (grad(u2)*grad(v2),definedon=[2])

# nitsche term

cf1 = 0.5*(grad(u1)[0] + grad(u2)[0])*(v1-v2)

cf2 = 0.5*(grad(v1)[0]+ grad(v2)[0])*(u1-u2)

cf3 = 20*order*(order+1)*(u1-u2)*(v1-v2)/h


a += SymbolicBFI (-cf1-cf2+cf3,BND, skeleton=True, definedon=[2])

a.Assemble()

f = LinearForm(V)
f += SymbolicLFI ( source * v1, definedon = [1])
f += SymbolicLFI ( source * v2, definedon = [2])
f.Assemble()

u = GridFunction(V)
u.vec.data = a.mat.Inverse(V.FreeDofs()) * f.vec

"""
u.components[0].Set(uex,BND)
u.components[1].Set(uex,BND)
res = f.vec.CreateVector()
res.data = f.vec - a.mat*u.vec
u.vec.data += a.mat.Inverse(V.FreeDofs()) * res
"""

"""
L2 projection
"""

Vpost = L2(mesh, order = order)
#Vpost = H1(mesh, order = order)
up = Vpost.TrialFunction()
vp = Vpost.TestFunction()
ap = BilinearForm(Vpost)
ap += SymbolicBFI(up*vp)
fp = LinearForm(Vpost)
fp += SymbolicLFI(u.components[0]*vp,definedon=[1])
fp += SymbolicLFI(u.components[1]*vp,definedon=[2])

ap.Assemble()
fp.Assemble()
up = GridFunction(Vpost)
up.vec.data = ap.mat.Inverse()*fp.vec

err = Integrate((up-uex)*(up-uex),mesh)
print(err)
Draw(up)

Draw(uex, mesh, "exact") # exact solution