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:
> > >   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
> > >
> >
>
> --
> 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
>