Divergence theorem on anisotropic diffusion problem

8 months ago by
Hi All,
The following problem is bugling my mind for a while now. I would appreciate if someone could shed light on it.
File attached: Schematic (15.4 KB)

Our domain is a bi-unit square plate with all Dirichlet BCs. A subdomain is marked in the center of our square ($\Omega_{\mathrm{subdomain}}$) and the interface facets of
subdomain are also marked separately (see figure attached).
After solving anisotropic diffusion equation with diffusivity tensor $D$, we want 
to calculate the following $M-tensor$ on nodal solutions:
\mathbf{M} = \mathrm{grad}[c] \otimes \mathrm{grad}[c]
and finally obtain the x-component of the following vector:
\int_{\Omega_{\mathrm{subdomain}}} \mathrm{div}[\mathbf{M}] \; \mathrm{d} \Omega
However, from divergence theorem, we know that:
\int_{\Omega_{\mathrm{subdomain}}} \mathrm{div}[\mathbf{M}] \; \mathrm{d} \Omega =
\int_{\Gamma^{+}_{\mathrm{subdomain}}} \mathbf{M}\widehat{\mathbf{n}} \; \mathrm{d} \Gamma
But I could not get the same results when I use RHS( denoted as method I in code) and  LHS( denoted as method II in code) of the above equation. I was able to get similar values for RHS and LHS for isotropic
diffusion but not for anisotropic diffusion (the values get closer when anisotropy decreases).
I have two questions:
1-what is the result for this discrepancy. Could our divergence theorem be violated in (Galerkin numeral setting) or the problem is with my implementation ?
2- Since I need to perform mesh refinement study on this problem, which result is reliable?
method I or method II
Thank you very much!
from dolfin import *
# Create mesh and define function space ;
mesh = UnitSquareMesh(30,30)
W = FunctionSpace(mesh,"Lagrange",2)
# Define volumetric source ;
f = Constant(0.0)
# Anisotropic diffusivity ;
theta = pi/6.0
c = cos(theta)
s = sin(theta)
d1 = 1000.0
d2 = 1.0
RMat = as_matrix([[c, s],[-s, c]])
D0 = as_matrix([[d1, 0],[0, d2]])
D = RMat * (D0 * RMat.T)
# Define boubdary conditions ;
u_B = Expression('sin(pi * x[0])', degree=5)
def bottom_boundary(x,on_boundary):
tol = 1E-14
return on_boundary and abs(x[1]) < tol
Gamma_0 = DirichletBC(W,u_B,bottom_boundary)
u_T = Constant(0)
def other_boundary(x,on_boundary):
tol = 1E-14
return on_boundary and abs(x[1]) > tol
Gamma_1 = DirichletBC(W,u_T,other_boundary)
bcs = [Gamma_0, Gamma_1]
# Define variational problem ;
u = TrialFunction(W)
v = TestFunction(W)
a = inner(D*nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx
# Compute solution ;
u = Function(W)
solve(a == L, u, bcs)
# Compute flux integrals ;
# Marking inner-square subdomain
class Omega6(SubDomain):
def inside(self, x, on_boundary):
return ( x[1] >= 1.0/3.0 and x[1] <= 2.0/3.0 ) and ( x[0] >= 1.0/3.0 and x[0] <= 2.0/3.0 )
materials = CellFunction('size_t', mesh)
subdomain6 = Omega6()
subdomain6.mark(materials, 6)
plot(materials, title = "Marked dubdomain")
dx = Measure("dx")[materials]
# Marking each facets of inner-square separately
loc = 1.0/3.0
leftSegment = AutoSubDomain(lambda x: near(x[0],loc) and x[1] >= loc and x[1] <= 1.0-loc)
bottomSegment = AutoSubDomain(lambda x: near(x[1],loc) and x[0] >= loc and x[0] <= 1.0-loc)
rightSegment = AutoSubDomain(lambda x: near(x[0],1.0-loc) and x[1] >= loc and x[1] <= 1.0-loc)
topSegment = AutoSubDomain(lambda x: near(x[1],1.0-loc) and x[0] >= loc and x[0] <= 1.0-loc)
markers = FacetFunction("size_t",mesh)
plot(markers,title = "marked facets")
dS = dS[markers]
# Functional M-tensor ;
M_tensor = outer(as_vector(grad(u)),as_vector(D*grad(u)))
# method I (sum int_{Gamma_i} Jn d\Gamma)
# Useful quantities
ex = Constant((1.0,0.0))
ey = Constant((0.0,1.0))
# Left surface
nVec = -ex # (nVec is unit normal vector)
M_Left = assemble( dot(ey,M_tensor('-')*nVec)*dS(2) )
# Bottom surface
nVec = -ey
M_Bottom = assemble( dot(ey,M_tensor('-')*nVec)*dS(3) )
# Right surface
nVec = ex
M_Right = assemble( dot(ey,M_tensor('-')*nVec)*dS(4) )
# Top surface
nVec = ey
M_Top = assemble( dot(ey,M_tensor('-')*nVec)*dS(5) )
print "method I = ", M_Left + M_Bottom +M_Top + M_Right
# method II \int_Omega1 div[J] d \Omega
div_M = (div(M_tensor)[1])*dx(6)
output_div_M = assemble(div_M)
print "method II = ", output_div_M


Community: FEniCS Project
This is strange, the divergence theorem should work for any tensor. I guess there might be a problem with the unclear identification of the '+' and '-' sides, see the discussion here: https://fenicsproject.org/qa/1553/and-in-interior-facets
written 8 months ago by Adam Janecka  
Do you see how this would be in addition to the fact that he is solving a non-well-posed-discrete system?
written 8 months ago by pf4d  
It wouldn't. It is just strange that the divergence theorem doesn't hold. It should be satisfied for any (no matter how constructed) tensor.
written 8 months ago by Adam Janecka  
... given that tensor field is differentiable, that is.
written 8 months ago by pf4d  

1 Answer

8 months ago by
I haven't looked at all you code, but I can see from your variational problem that you are solving a non-symmetric system; you can find a very good explanation of this problem and solution in Chapter 11 of Zienkiwicz and Taylor (2005).

Or you can check how others have solved this problem with FEniCS here:


and may help you with your mesh refinements later on too.
.....     Did this help?
written 8 months ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »