Verify the divergence theorem

mathematics
1

Hi everyone!
I need some help trying to verify the divergence theorem in fenics. I am using dolfin 2018

The MWE is as follows:

from dolfin import *

# Create a mesh
mesh = UnitSquareMesh(32, 32)

class rectangle(SubDomain):
    def inside(self,x, on_boundary):
        return between(x[0],(0.25,0.75)) and between(x[1],(0.25,0.75))
Rec = rectangle()

class InnerBoundary(SubDomain):
    def inside(self,x, on_boundary):
        incell = Rec.inside(x,True)
        return incell and bool(near(x[0],0.25) or near(x[0],0.75) or near(x[1],0.25) or near(x[1],0.75))
inner_bndr = InnerBoundary()

# Define unit cell facet indicator function
inner_bndrindicator = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
inner_bndrindicator.set_all(0)
inner_bndr.mark(inner_bndrindicator,1)

# Define domain indicator function
domainindicator = MeshFunction('size_t', mesh, mesh.topology().dim())
domainindicator.set_all(1)
Rec.mark(domainindicator,0)

# Define the finite element space
V = FunctionSpace(mesh, "CG", 2)

# Define the expression
class MyExpression(UserExpression):
    def eval(self, value, x):
        value[0] = exp(x[0])
        if x[0] >= 0.25 and x[0] <= 0.75 and x[1] >= 0.25 and x[1] <= 0.75:
            value[0] = 0.0

    def value_shape(self):
        return ()

# Define and interpolate the expression
T = MyExpression(degree=2)
T_orig = interpolate(T, V)

## Test the divergence theorem on the fluid domain
# Define the normal vector on the boundary
n = FacetNormal(mesh)
dS_facet = Measure('dS', domain=mesh, subdomain_data=inner_bndrindicator,subdomain_id=1) # fluid solid interface
ds = Measure('ds', domain=mesh)  # exterior boundary of the unit mesh
facetorientationform = Constant(0)*dx(domain=mesh, subdomain_data=domainindicator)
flux = assemble(inner(n('+'), grad(T_orig)('+'))*dS_facet+facetorientationform+inner(n, grad(T_orig))*ds)

dx = Measure('dx', domain=mesh, subdomain_data=domainindicator,subdomain_id=1)
volume_f_integral = assemble(div(grad(T_orig))*dx) 

print("Volume integral of fluid div(grad(T)):", volume_f_integral)
print("Surface integral of fluid dot(n, grad(T)):",flux)
print("Difference (should be close to zero):", volume_f_integral - flux)

image

In the above picture, the temperature field is on the left. On the right is the domain indicator meshfunction, which it differentiates the fluid domain(yellow) and solid domain(blue), in the center, the blue rectangle is at T=0 and acts like some BC.

Running the above code it gives:

Volume integral of fluid div(grad(T)): -456.8569440160415
Surface integral of fluid dot(n, grad(T)): -321.87865028985516
Difference (should be close to zero): -134.97829372618634

May I know why, the surface integral and the volume integral are not matched? Is it because of the T=0 rectangle placed in the center making it unphysical?

Another question is when we change the mesh resolution from 32 to 64 why the integral also seems to double?

Volume integral of fluid div(grad(T)): -887.7423581283879
Surface integral of fluid dot(n, grad(T)): -645.0589192858026
Difference (should be close to zero): -242.68343884258525

Thanks for any help!!