Hello,
I have coded the finite element method by hand. In my code I use one point Gaussian quadrature for a right hand side. I evaluate the load function at each triangles centroid, then multiply by the triangle area and 1/3 (representing the hat function evaluated at the triangle centroid), and accumulate into a vector.
In FEniCS, I have gathered that the using
metadata={'quadrature_degree': 1}
as an argument to my measure is supposedly doing this. I was made aware of the caveat regarding Expressions in the following Q&A
http://fenicsproject.org/qa/9096/cant-reproduce-integration-function-using-finite-elements
This leads me to believe I should be able to recover my assembled load vector by supplying an Expression in FEniCS with degree=n, and then increase n to achieve some desired error tolerance between the two vectors. Of course, this turns out to not be true...
Does this mean FEniCS does something else besides one-point Gaussian quadrature for this quadrature degree for elements in FunctionSpace(mesh, 'CG', 1)? Or have I missed something else here?
Here is the FEniCS code I use to assemble the load vector:
nx = 4
ny = 4
mesh = UnitSquareMesh(nx-1, ny-1)
V = FunctionSpace(mesh, 'Lagrange', 1)
d2v = dof_to_vertex_map(V)
u = TrialFunction(V)
v = TestFunction(V)
f = Expression('(-8.*pi*pi) * ( sin((2.0*pi*x[0]) - (0.75*pi)) * sin((2.0*pi*x[1]) - (0.75*pi)) )', degree=6)
L = inner(f,v)*dx(metadata={'quadrature_degree': 1})
b = assemble(L)
print assemble(L).array()[d2v]
[-0.63313542 -2.49406667 1.88016875 1.24703334 0.14618852 0.61389792
-0.12695102 -0.63313542 0.63313542 -0.12695102 1.88016875 -2.38635315
1.24703334 0.61389792 -2.49406667 0.63313542]
The output from my own FEM code using one-point Gaussian quadrature is:
-1.45105684e+00 -1.81659774e+00 2.53657278e+00 7.31081807e-01
-1.81659774e+00 2.15992510e+00 6.66133815e-16 -3.43327353e-01
2.53657278e+00 6.10622664e-16 -2.15992510e+00 -3.76647679e-01
7.31081807e-01 -3.43327353e-01 -3.76647679e-01 -1.11067754e-02
