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