Does project not conserve integral?
Say I want to project a function $u$u from the fine space $V_f$Vƒ onto the coarser space $V_c$Vc , where the fine space is created by refining the coarse space.
Then the following problem is solved:
$inner(Pu,v)*dx=inner(u,v)*dx$inner(Pu,v)*dx=inner(u,v)*dx for all $v\in V_c$v∈Vc
Now, the piecewise constant 1 function should be in both spaces. So I would expect that the integral of the projection over the coarse space should be (up to some small error) the same as the integral of $u$u over the fine space.
However, as seen in the following example this is not the case.
This gives the output
from __future__ import print_function from dolfin import * # Form compiler options parameters["form_compiler"]["optimize"] = True parameters["form_compiler"]["cpp_optimize"] = True mesh = UnitIntervalMesh(3) # mesh = UnitSquareMesh(3,3) refined_mesh = refine(mesh) func = Expression("exp(x)",degree=5) V = FunctionSpace(mesh, "DG",1) V1 = FunctionSpace(refined_mesh, "DG",1) dx = Measure("dx",domain = mesh) dx1 = Measure("dx",domain = refined_mesh) u_coarse = Function(V) u_fine = Function(V1) u_fine.interpolate(func) print(str(assemble(u_fine*dx1))) u_coarse = project(u_fine,V) print(str(assemble(u_coarse*dx)))
I was wondering if there are too less points used in order to solve the linear problem, so I tried to increase the quadrature degree but this seems not to have any effect. Am I missing something obvious?
u_coarse = project(u_fine,V)
basically translates to
u_coarse = Function(V) solve(TrialFunction(V)*TestFunction(V)*dx == u_fine*TestFuntion(V)*dx, u_coarse)
u_fine*TestFunction(V)*dx is evaluated by taking quadrature over the cells of the coarse mesh. Observe that integrand is not polynomial on the coarse cells, it is piece-wise polynomial. This is a limitation of the code generation and the assembly chain how it currently exists. The bad thing is that the approximate interpolation of the nonpolynomial integrand happens behind the scene without user's consent. This was already discussed on a related thread, https://bitbucket.org/fenics-project/uflacs-deprecated/issues/49/representation-fails-with-coefficient-from, and suggested remedy was to require explicit UFL request for interpolation on non-matching mesh, e.g. using
NodalInterpolant(f_fine), otherwise fail with an error message.
The limitation that one cannot integrate on coarse cells iterating over fine cells goes deep to the current design of the assembly chain. So don't expect that there will be an easy way around this soon. In a specific cases one can sometimes find a workaround, for example assembling a 1-form with DP0 test functions (i.e., one degree of freedom per cell) can be done on fine mesh and then summing up the contributions from the fine cells to owning coarse cells.