Hello. I am wondering how FEniCS comptues interior facet integrals. In this example, I am computing the integral of $f(x,y)=1$ under the interior facets for A, and $f(x,y)=x$ for B.
from dolfin import *
mesh = RectangleMesh(0.0, 0.0, 1.0, 1.0, 2, 2)
V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
f = project(Expression("x[0]"), V)
A = assemble(avg(v) * dS)
B = assemble(f * avg(v) * dS)
print A.array()
print sum(A.array())
print B.array()
print sum(B.array())
This results in the output:
[ 0. 0.60355339 0.60355339 0.35355339 1.70710678 0.35355339
0.60355339 0.60355339 0. ]
4.82842712475
[ 0. 0.10059223 0.24285113 0.05892557 0.85355339 0.29462783
0.36070226 0.50296116 0. ]
2.41421356237
When I compute A, it seems like FEniCS computes the integral under each facet (which is just a rectangle in this case) then half of that value is allotted to one vertex connected to that facet and the other half is allotted to the other vertex. For example, the value associated with the bottom left vertex is $ \frac{ \sqrt{ \frac{1}{2} }} {2} $, which is half the area of a rectangle from the bottom left vertex to the center vertex with height 1.
However, I'm not sure what's going on in B when I take the facet integrals under $f(x,y)=x$. The sum of B.array() is what I expect. I assume that the value at each vertex is just a weighted sum of the connected interior facet integrals, but I don't understand how these weights are assigned.
Thanks,
Jake
