Comment 6 for bug 901069

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=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.
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/901069/+subscriptions
>