Hi

I have a simple question in a time-dependent diffusion problem.

Why flux integration on the interior domain has different result from that defined in the From equation as a constant flux? This difference is such that in through time steps, computed flux reaches a constant (~0.17 in this setting) but different from Neumann boundary condition (0.2) that specified in the right hand side Form. Also, computed flux trend is affected by diffusion constant.

The complete code is

```
from fenics import *
from mshr import *
import matplotlib.pyplot as plt
circle1= Circle(Point(30, 30), 30)
circle2= Circle(Point(30, 30), 15)
Domain= circle1 - circle2
mesh= generate_mesh(Domain, 15)
# Neumann b.c. will defined in inner boundary
B= CompiledSubDomain('pow(x[0]-30, 2) + pow(x[1]-30, 2) < pow(15+0.5, 2)')
B_Marker= MeshFunction('size_t', mesh, 1)
B.mark(B_Marker, 1)
ds= Measure('ds', domain= mesh, subdomain_data= B_Marker)
V= FunctionSpace(mesh, 'P', 1)
u0= interpolate(Constant(1), V)
u= TrialFunction(V)
v= TestFunction(V)
D= Constant(1)
g= Constant(2)
dt= Constant(0.1)
a= u*v*dx + dt*D*dot(grad(u), grad(v))*dx
L= u0*v*dx + dt*g*v*ds(1)
n= FacetNormal(mesh)
flux= dt*D*dot(grad(u0), n)*ds(1)
DS= assemble(1*ds(1))
A= assemble(a)
u= Function(V)
flux_array= []
for i in range(200):
b= assemble(L)
solve(A, u.vector(), b)
u0.assign(u)
flux_array.append(assemble(flux)/DS)
plt.plot(flux_array)
```