Does project not conserve integral?


86
views
3
28 days ago by
I observed that using project to project from a finer grid onto a coarser grid does not conserve the integral of a quantity.
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$vVc
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.

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[0])",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)))
​
This gives the output 1.72225749247
                                      1.73416246012


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?
Community: FEniCS Project

1 Answer


6
28 days ago by
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)


Now 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.

Thank you, I didn't find the issue on bitbucket when I was digging..

However the code provided by you in the latest post of the discussion on bitbucket.org did not work for me, since "NodalInterpolant" seems not to be known.
I also could not find any documentation thereof. Could you please elaborate this a little bit for me?
written 28 days ago by Lukas O.  
1
NodalInterpolant does not exist yet. It is a proposed remedy for the problem with hidden interpolation. Remind the Python Zen: explicit is better than implicit.

# Should raise error "not able to integrate; integrand not polynomial"
assemble(u_on_mesh_0*dx(mesh_1))

# This will give the current behaviour with hidden interpolation
# Eventually simplified inteface is possible
el = FiniteElement(...)
interpolant = NodalInterpolant(u_on_mesh_0, ufl.FunctionSpace(el, mesh_1)
assemble(interpolant*dx(mesh_1))​
written 28 days ago by Jan Blechta  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »