Hi, I'm very new to Fenics so this question is mostly about trying to understand how the automatic integration works, because I'm seeing some strange behaviour that I don't quite understand.
For example code, lets just say I'm solving the Poisson problem
mesh = UnitSquareMesh(N, N) # N = 4,8,16,32
finemesh = UnitSquareMesh(128, 128)
exact = Expression('sin(pi*x[0])*sin(pi*x[1])',domain=mesh);
f = Expression('2*pow(pi,2.0)*sin(pi*x[0])*sin(pi*x[1])',domain=mesh);
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
F = f*v*dx # f is source term
sol = Function(V)
bc = DirichletBC(V, g, Dir_boundary); #g is boundary conditions
solve(a == F, sol, bc)
Now I want to show that sol satisfies the Galerkin orthogonality condition, so I compute the left and right hand side of the condition:
#Compute left hand side \int grad(u) \dot grad(v)
R = inner(grad(sol), grad(v))*dx(finemesh);
left[ii] = assemble( R )
#Compute right hand side \int f*v
right[ii] = assemble(v*f*dx(finemesh));
where I use a finer mesh to avoid any error using quadrature. Now I should be able plug in any H^1_0 function for the test function v, and as N gets bigger the ratio ratio = left/right should go to 1.
If I use
v = Expression('x[0]*(x[0]-1.0)*x[1]*(x[1]-1.0)', domain=finemesh, degree=4)
in the above code, the ratio goes to 1.
If I use
v = interpolate(Expression('x[0]*(x[0]-1.0)*x[1]*(x[1]-1.0)', domain=finemesh, degree=4), V)
the ratio also goes to one, which makes sense since this v is in the test space.
However if I use
v1 = Expression('x[0]*(x[0]-1.0)*x[1]*(x[1]-1.0)', domain=finemesh, degree=4)
v = v1 - interpolate(v1, V)
the ratio goes to something like -.5. This last one is important because it is similar to an error, and one of the things I'm trying to code up is an a posteriori error representation formula. But mostly something is happening in the integration code that I don't understand, and I'd like to. Can anyone tell me what is happening with this last integration?
Thanks in advance for any help.
