Comment 8 for bug 901069

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
>