How to handle Implicit Euler method when replacing domain and mesh

I tried to implement the implicit Euler method for a time dependent problem.
[tex]\partial_t c+\nabla\cdot(cu-k\nabla c)=f[/tex].
Please see the attached file.

The code run well with the unit square domain.

mesh = Mesh(unit_square.GenerateMesh(maxh=0.5) )

When I replaced the unit square on line 32 (above) by a new domain (a unit square with two subdomains) by line 33.

mesh = Mesh(geo.GenerateMesh(maxh=0.5))

and line 83

a += IfPos(u*n,1,0)*B_ds_out*ds(skeleton=True,definedon=mesh.Boundaries("right|top|left|bottom"))
by line 84

a += IfPos(u*n,1,0)*B_ds_out*ds(skeleton=True,definedon=mesh.Boundaries("GammaS|GammaD")) I think this would give a similar results. But the new one doesn’t make any sense.

Could you please tell me what the mistakes I made? Do I need to change anything in the bilinear and linear forms? Any help would be appreciated.
Thank you so much.

Attachment: neumann.ipynb