Cannot do L2 projection between functions on different meshes

[Lizao (Larry) Li]

I am trying to write a geometric multigrid code and I cannot do L2 projection from a fine mesh to a coarse (the restriction operator).
Here is an simplest code to demostrate this problem:

from dolfin import *
mesh_coarse = UnitSquare(2, 2)
V_coarse = FunctionSpace(mesh_coarse, "CG", 1)
mesh_fine = UnitSquare(4, 4)
V_fine = FunctionSpace(mesh_fine, "CG", 1)
f = Expression("sin(x[0])")
f_fine = interpolate(f, V_fine)
u = TrialFunction(V_coarse)
v = TestFunction(V_coarse)
lhs = u*v*dx
A = assemble(lhs)
rhs = f*v*dx
b = assemble(rhs, mesh = mesh_fine) ## This line gives an error because v is an "Argument"
f_coarse = Function(V_coarse)
solve(A, f_coarse.vector(), b)

If I do the assemble without specifying the mesh, f_fine will first be interpolated into V_coarse and then product with v, the end result is just the interpolant but not the L2 projection.
If v=interpolate(f, V_coarse), then assemble(f*v*dx, mesh = V_fine) works and gives the correct answer. However when v is a test function, dolfin complains it cannot be interpolated (which is understandable, because v is not given a concrete value yet).

It seems that the only operation which can go between function spaces defined on different meshes is interpolate so far. It would be nice to implement the above operation which is very valuable in the implementation of geometric multigrids. In terms of coding, it seems to be a matter of passing mesh=mesh_fine into each of the form assembly loops, when v is specified.

[Johan Hake (johan-hake)]

Try:

  project(f_fine, V_coarse)

Which essentially does what you describe, but using

  b = assemble(rhs)

AFAIK this is the L2 projection of f unto V_coarse.

Also the mesh argument to assemble is only used when the form does not contain any mesh.

[Lizao (Larry) Li (creatorlarryli)]

Thank you for the quick reply. However both of the approaches you mentioned don't give the right answer.
The following is a demonstration:

from dolfin import *
mesh_coarse = UnitInterval(2)
V_coarse = FunctionSpace(mesh_coarse, "CG", 1)
mesh_fine = UnitInterval(16)
V_fine = FunctionSpace(mesh_fine, "CG", 1)
# The following function vanishes on the coarse grid but has a non-zero L2 projection onto the piecewise linear space on the coarse grid
f = Expression("sin(2.0*pi*x[0])")
f_fine = interpolate(f, V_fine)
f_coarse = interpolate(f, V_coarse)

# Approach #1 you suggested
u = TrialFunction(V_coarse)
v = TestFunction(V_coarse)
lhs = u*v*dx
A = assemble(lhs)
rhs = f_fine*v*dx
b = assemble(rhs) # Without specify the mesh, f_fine is interpolated into V_coarse first, giving the wrong answer
f_false_l2 = Function(V_coarse)
solve(A, f_false_l2.vector(), b) #The interpolant is L2-projected, nothing happens
print(sqrt(assemble(f_false_l2**2*dx)))

# Approach #2 you suggested
f_project = project(f_fine, V_coarse) #Again, the interpolnt is projected. The result is no better than the interpolant.
print(sqrt(assemble(f_project**2*dx)))

# The only correct way I know to get the right answer, very inefficient
# b_correct = assemble(f_fine*v*dx, mesh = mesh_fine), calculated by manually loop assembly
b_correct = Vector(3)
for i in range(3):
    basis = Function(V_coarse)
    basis.vector()[i] = 1.0
    b_correct[i] = assemble(f_fine*basis*dx, mesh = mesh_fine)
f_correct = Function(V_coarse)
solve(A, f_correct.vector(), b_correct)
print(sqrt(assemble(f_correct**2*dx))) # This gives the right non-zero answer.

That's why I think this is a bug. The code does not do behave as we intend it to.

[Doug Arnold (dnarnold)]

I have to agree with Larry. The project function does not work to
project from a finite element on a fine mesh to a finite element on
a coarser mesh. Here is a simple example. Take the hat function
on the 1D mesh (0, 1/4, 1/2, 3/4, 1) which is 1 at 1/4 are 0 at
the other nodes, and project it onto the continuous piecewise
linear functions on the coarser mesh (0, 1/2, 1). Simple use of
the project function claims that the projection is identically
zero, but that is obviously not correct. Is there an easy way
to correctly compute this in FEniCS?

from dolfin import *
mesh = UnitInterval(4)
V = FunctionSpace(mesh, 'CG', 1)
u = Function(V)
u.vector()[1] = 0.0 # nodal values of u are 0, 1, 0, 0, 0
mesh0 = UnitInterval(2)
V0 = FunctionSpace(mesh0, 'CG', 1)
u0 = project(u, V0)
print u0.vector().array()

The output is [0., 0., 0.], while the correct answer is [.625, .25, -.125].

[Anders Logg (logg)]

On Thu, Dec 08, 2011 at 08:23:03PM -0000, Doug Arnold wrote:
> I have to agree with Larry. The project function does not work to
> project from a finite element on a fine mesh to a finite element on
> a coarser mesh. Here is a simple example. Take the hat function
> on the 1D mesh (0, 1/4, 1/2, 3/4, 1) which is 1 at 1/4 are 0 at
> the other nodes, and project it onto the continuous piecewise
> linear functions on the coarser mesh (0, 1/2, 1). Simple use of
> the project function claims that the projection is identically
> zero, but that is obviously not correct. Is there an easy way
> to correctly compute this in FEniCS?
>
> from dolfin import *
> mesh = UnitInterval(4)
> V = FunctionSpace(mesh, 'CG', 1)
> u = Function(V)
> u.vector()[1] = 0.0 # nodal values of u are 0, 1, 0, 0, 0

Should be 1.0 here, but it doesn't matter.

> mesh0 = UnitInterval(2)
> V0 = FunctionSpace(mesh0, 'CG', 1)
> u0 = project(u, V0)
> print u0.vector().array()
>
> The output is [0., 0., 0.], while the correct answer is [.625, .25,
> -.125].

The reason is that when a function defined on another mesh, it is
interpolated locally to the mesh of the function space used to define
the form (which in this case is the function space into which the
function is being projected). So 0, 1, 0, 0, 0 will be interpolated
into 0, 0, 0 and then projected.

It might be possible to get around this by experimenting with
QuadratureElement and some tweaking of the eval-chain of Expressions
in DOLFIN, but I'm not sure.

--
Anders

[Jan Blechta (blechta)]

> The reason is that when a function defined on another mesh, it is
> interpolated locally to the mesh of the function space used to define
> the form (which in this case is the function space into which the
> function is being projected). So 0, 1, 0, 0, 0 will be interpolated
> into 0, 0, 0 and then projected.

I checked code generated by project(.., form_compiler_parameters={'representation': 'quadrature'}). Here's relevant piece of compiled rhs of projection to CG1:
_____________________________________________________________________________
    // Loop quadrature points for integral.
    // Number of operations to compute element tensor for following IP loop = 24
    for (unsigned int ip = 0; ip < 2; ip++)
    {

      // Coefficient declarations.
      double F0 = 0.0;

      // Total number of operations to compute function values = 4
      for (unsigned int r = 0; r < 2; r++)
      {
        F0 += FE0[ip][r]*w[0][r];
      }// end loop over 'r'

      // Number of operations for primary indices: 8
      for (unsigned int j = 0; j < 2; j++)
      {
        // Number of operations to compute entry: 4
        A[j] += FE0[ip][j]*F0*W2[ip]*det;
      }// end loop over 'j'
    }// end loop over 'ip'
____________________________________________________________
FE0 are basis (of new element) values at quadrature points. W2 are quadrature weights.

It is obvious that form requires dof values w[][] of coefficient. But these are dofs with respect to new element. So they must be computed by interpolation in DOLFIN.

So it seems that projection is implemted with pre-interpolation in DOLFIN because UFC requires it. Sounds like a major flaw in UFC.

[Martin Sandve Alnæs (martinal)]

It should rather be interpolated to a quadrature element. That is dolfin
functionality, ufc is fine.

Martin
Den 12. apr. 2013 17:06 skrev "Jan Blechta" <email address hidden>
følgende:

> > The reason is that when a function defined on another mesh, it is
> > interpolated locally to the mesh of the function space used to define
> > the form (which in this case is the function space into which the
> > function is being projected). So 0, 1, 0, 0, 0 will be interpolated
> > into 0, 0, 0 and then projected.
>
> I checked code generated by project(..,
> form_compiler_parameters={'representation': 'quadrature'}). Here's relevant
> piece of compiled rhs of projection to CG1:
>
> _____________________________________________________________________________
> // Loop quadrature points for integral.
> // Number of operations to compute element tensor for following IP
> loop = 24
> for (unsigned int ip = 0; ip < 2; ip++)
> {
>
> // Coefficient declarations.
> double F0 = 0.0;
>
> // Total number of operations to compute function values = 4
> for (unsigned int r = 0; r < 2; r++)
> {
> F0 += FE0[ip][r]*w[0][r];
> }// end loop over 'r'
>
> // Number of operations for primary indices: 8
> for (unsigned int j = 0; j < 2; j++)
> {
> // Number of operations to compute entry: 4
> A[j] += FE0[ip][j]*F0*W2[ip]*det;
> }// end loop over 'j'
> }// end loop over 'ip'
> ____________________________________________________________
> FE0 are basis (of new element) values at quadrature points. W2 are
> quadrature weights.
>
> It is obvious that form requires dof values w[][] of coefficient. But
> these are dofs with respect to new element. So they must be computed by
> interpolation in DOLFIN.
>
> So it seems that projection is implemted with pre-interpolation in
> DOLFIN because UFC requires it. Sounds like a major flaw in UFC.
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/901069
>
> Title:
> Cannot do L2 projection between functions on different meshes
>
> Status in DOLFIN:
> Confirmed
>
> Bug description:
> I am trying to write a geometric multigrid code and I cannot do L2
> projection from a fine mesh to a coarse (the restriction operator).
> Here is an simplest code to demostrate this problem:
>
> from dolfin import *
> mesh_coarse = UnitSquare(2, 2)
> V_coarse = FunctionSpace(mesh_coarse, "CG", 1)
> mesh_fine = UnitSquare(4, 4)
> V_fine = FunctionSpace(mesh_fine, "CG", 1)
> f = Expression("sin(x[0])")
> f_fine = interpolate(f, V_fine)
> u = TrialFunction(V_coarse)
> v = TestFunction(V_coarse)
> lhs = u*v*dx
> A = assemble(lhs)
> rhs = f*v*dx
> b = assemble(rhs, mesh = mesh_fine) ## This line gives an error because
> v is an "Argument"
> f_coarse = Function(V_coarse)
> solve(A, f_coarse.vector(), b)
>
> If I do the assemble without specifying the mesh, f_fine will first be
> interpolated into V_coarse and then product with v, the end result is just
> the interpolant but not the L2 projection.
> If v=interpo...

[Jan Blechta (blechta)]

On Fri, 12 Apr 2013 15:17:32 -0000
Martin Sandve Alnæs <email address hidden> wrote:
> It should rather be interpolated to a quadrature element. That is
> dolfin functionality, ufc is fine.

I see. I guess that logic should be that it is first interpolated to
quadrature element of degree corresponding to new space if it has none
or different space.

Do you know what quadrature element degree means? There have been
changes I think.

Generally speaking, shouldn't same operation be performed (i.e.
interpolation to quadrature element and tabulating tensor through
quadrature element version if element of coefficient is non-matching)
while assembling arbitrary form with coefficients? It would possibly
help with convergence to many users which are unaware of this subtle
issue.

Jan

>
> Martin
> Den 12. apr. 2013 17:06 skrev "Jan Blechta"
> <email address hidden> følgende:
>
> > > The reason is that when a function defined on another mesh, it is
> > > interpolated locally to the mesh of the function space used to
> > > define the form (which in this case is the function space into
> > > which the function is being projected). So 0, 1, 0, 0, 0 will be
> > > interpolated into 0, 0, 0 and then projected.
> >
> > I checked code generated by project(..,
> > form_compiler_parameters={'representation': 'quadrature'}). Here's
> > relevant piece of compiled rhs of projection to CG1:
> >
> > _____________________________________________________________________________
> > // Loop quadrature points for integral.
> > // Number of operations to compute element tensor for following
> > IP loop = 24
> > for (unsigned int ip = 0; ip < 2; ip++)
> > {
> >
> > // Coefficient declarations.
> > double F0 = 0.0;
> >
> > // Total number of operations to compute function values = 4
> > for (unsigned int r = 0; r < 2; r++)
> > {
> > F0 += FE0[ip][r]*w[0][r];
> > }// end loop over 'r'
> >
> > // Number of operations for primary indices: 8
> > for (unsigned int j = 0; j < 2; j++)
> > {
> > // Number of operations to compute entry: 4
> > A[j] += FE0[ip][j]*F0*W2[ip]*det;
> > }// end loop over 'j'
> > }// end loop over 'ip'
> > ____________________________________________________________
> > FE0 are basis (of new element) values at quadrature points. W2 are
> > quadrature weights.
> >
> > It is obvious that form requires dof values w[][] of coefficient.
> > But these are dofs with respect to new element. So they must be
> > computed by interpolation in DOLFIN.
> >
> > So it seems that projection is implemted with pre-interpolation in
> > DOLFIN because UFC requires it. Sounds like a major flaw in UFC.
> >
> > --
> > You received this bug notification because you are a member of
> > DOLFIN Core Team, which is subscribed to DOLFIN.
> > https://bugs.launchpad.net/bugs/901069
> >
> > Title:
> > Cannot do L2 projection between functions on different meshes
> >
> > Status in DOLFIN:
> > Confirmed
> >
> > Bug description:
> > I am trying to write a geometric multigrid code and I cannot do L2
> > projection from a fine mesh to a coarse (the restriction operator).
>...

[Kristian B. Ølgaard (k.b.oelgaard)]

On 12 April 2013 17:50, Jan Blechta <email address hidden> wrote:

> On Fri, 12 Apr 2013 15:17:32 -0000
> Martin Sandve Alnæs <email address hidden> wrote:
> > It should rather be interpolated to a quadrature element. That is
> > dolfin functionality, ufc is fine.
>
> I see. I guess that logic should be that it is first interpolated to
> quadrature element of degree corresponding to new space if it has none
> or different space.
>
> Do you know what quadrature element degree means? There have been
> changes I think.
>

See Section 26.5 (26.5.2) in the FEniCS book.

> Generally speaking, shouldn't same operation be performed (i.e.
> interpolation to quadrature element and tabulating tensor through
> quadrature element version if element of coefficient is non-matching)
> while assembling arbitrary form with coefficients? It would possibly
> help with convergence to many users which are unaware of this subtle
> issue.
>
> Jan
>
> >
> > Martin
> > Den 12. apr. 2013 17:06 skrev "Jan Blechta"
> > <email address hidden> følgende:
> >
> > > > The reason is that when a function defined on another mesh, it is
> > > > interpolated locally to the mesh of the function space used to
> > > > define the form (which in this case is the function space into
> > > > which the function is being projected). So 0, 1, 0, 0, 0 will be
> > > > interpolated into 0, 0, 0 and then projected.
> > >
> > > I checked code generated by project(..,
> > > form_compiler_parameters={'representation': 'quadrature'}). Here's
> > > relevant piece of compiled rhs of projection to CG1:
> > >
> > >
> _____________________________________________________________________________
> > > // Loop quadrature points for integral.
> > > // Number of operations to compute element tensor for following
> > > IP loop = 24
> > > for (unsigned int ip = 0; ip < 2; ip++)
> > > {
> > >
> > > // Coefficient declarations.
> > > double F0 = 0.0;
> > >
> > > // Total number of operations to compute function values = 4
> > > for (unsigned int r = 0; r < 2; r++)
> > > {
> > > F0 += FE0[ip][r]*w[0][r];
> > > }// end loop over 'r'
> > >
> > > // Number of operations for primary indices: 8
> > > for (unsigned int j = 0; j < 2; j++)
> > > {
> > > // Number of operations to compute entry: 4
> > > A[j] += FE0[ip][j]*F0*W2[ip]*det;
> > > }// end loop over 'j'
> > > }// end loop over 'ip'
> > > ____________________________________________________________
> > > FE0 are basis (of new element) values at quadrature points. W2 are
> > > quadrature weights.
> > >
> > > It is obvious that form requires dof values w[][] of coefficient.
> > > But these are dofs with respect to new element. So they must be
> > > computed by interpolation in DOLFIN.
> > >
> > > So it seems that projection is implemted with pre-interpolation in
> > > DOLFIN because UFC requires it. Sounds like a major flaw in UFC.
> > >
> > > --
> > > You received this bug notification because you are a member of
> > > DOLFIN Core Team, which is subscribed to DOLFIN.
> > > https://bugs.launchpad.net/bugs/901069
> > >
> > > Title:
> > > Can...