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