DG VectorFunctionSpace

I might be missing something here, but it looks like you compare
different things: the full jump and the jump in only the normal
component.

--
Anders

On Sat, Dec 11, 2010 at 08:51:26PM -0000, JG wrote:
> New question #137463 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/137463
>
> Hello -- I seem to get different answers when I project a discontinuous vector function component by component or using the vector DG function space. I expected the difference to be machine precision. What could I be doing wrongly? Thanks, --Jay
>
> from dolfin import *
>
> degree=1
> n=10
>
> mesh = UnitSquare(n,n)
> Q = Expression(("(x[1]>=0.5? 10.0 : 0.0)","(x[0]>=0.5? 100.0 : 0.0)"))
>
> VDG = VectorFunctionSpace(mesh, "DG", degree)
> DG = FunctionSpace(mesh, "DG", degree)
>
> q = project(Q, VDG, degree)
>
> qx = project(Q[0], DG, degree)
> qy = project(Q[1], DG, degree)
>
> print 'q and (qx,qy) differs by: %e, %e' % \
> ( sqrt(abs(assemble( dot(q[0]-qx,q[0]-qx)*dx, mesh=mesh))), \
> sqrt(abs(assemble( dot(q[1]-qy,q[1]-qy)*dx, mesh=mesh))) )
>
> n = FacetNormal(mesh)
> print 'jump(q)=%e, jump(qx)=%e, jump(qy)=%e' %\
> ( sqrt(abs(assemble( jump(q,n)*jump(q,n)*dS, mesh=mesh))), \
> sqrt(abs(assemble( jump(qx)*jump(qx)*dS, mesh=mesh))), \
> sqrt(abs(assemble( jump(qy)*jump(qy)*dS, mesh=mesh))) )
>
> Output:
>
> q and (qx,qy) differs by: 5.215406e-08, 3.885747e-07
> jump(q)=3.516078e-14, jump(qx)=1.931411e-07, jump(qy)=2.697398e-06
>
>
>
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp

I expected q and (qx,qy) to be the same. That is, when I projected Q into VDG space, I thought the result would equal the projection of Q[0] into the scalar DG space together with the same projection of Q[1]. So I was surprised when qx came out not equal to Q[0]. (Sorry for confusing you with the jump outputs, please ignore that.) --Jay

I get machine precision using the quadrature representation for the code generation. Just add

    parameters["form_compiler"]["representation"] = "quadrature"

somewhere near the top of the file.

The tensor contraction representation must be introducing round-off errors (which are not helped using parameters["form_compiler"]["precision"] ).

Thanks, that works! (Although I don't understand why it worked... If there is documentation on "quadrature representation for the code generation", please point me to it.) --Jay

The 'quadrature representation' is just the standard quadrature loop implementation. Details of the quadrature code generation and the quadrature optimisations that are activated by

    parameters["form_compiler"]["optimize"] = True

are described in

    http://dx.doi.org/10.1145/1644001.1644009

Thanks.

Sorry to continue with this again. But after I use your suggested quadrature representation parameter, I notice that jumps of both components of the projected functions are zero to machine precision. How can the projection of a discontinuous function into a DG space have zero jumps?

It's because you haven't specified a function space for Q. If a space is not provided, FFC will interpolate the function in a continuous Lagrange space, with a guess at the polynomial order to preserve the convergence order. If you want to interpolate Q in a discontinuous space, try

    VDG = VectorFunctionSpace(mesh, "DG", degree)
    Q = Expression(("(x[1]>=0.5? 10.0 : 0.0)","(x[0]>=0.5? 100.0 : 0.0)"), element = VDG.ufl_element())

Actually, this may still be problematic since if the discontinuity is along a cell edge, the interpolation is not "well defined" when machine precision is taken into account. The code will be attempting to evaluate the function at the discontinuity.

Another approach is to used what we have coined a 'quadrature element'. This is 'space' in which functions can only be evaluated at quadrature points, so Q would not be interpolated, but evaluated at quadrature points.