ffc team mailing list archive

[Kristian Oelgaard <k.b.oelgaard@xxxxxxxxx>]

On 1 February 2010 15:40, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
> Marie Rognes wrote:
> > Garth N. Wells wrote:
> > > Marie Rognes wrote:
> > >
> > > > Garth N. Wells wrote:
> > > >
> > > > > There seems to be a nasty bug in the quadrature optimisation. I can
> > > > > compute a result, but it's wrong. I can try to look into it in more
> > > > > detail later, but I'm using a combination of BDM and P0 functions which
> > > > > might be enough of a hint for someone.
> > > > Could you give a test script?
> > >
> > >
> > > Below is a script which reproduces a bug when the quadrature
> > > optimisation is turned on. The norm of vectors assembled with and
> > > without using optimisation dffer.
> >
> > Same problem arises with CG.
> >
> > This looks like Harish's old bug where:
> >
> >       un = dot(u0, n)/2.0
> >
> > does not work, but
> >
> >    un   = 0.5*dot(u0, n)
> >
> > does.
>
> It's a strange one. I thought that Kristian had fixed it.
>
> Interesting is that
>
>  un   = dot(u0, n)/2.0
>  F    = 1.0/(s0**2 + 1.0)
>
> doesn't work, but
>
>  un   = dot(u0, n)/2.0
>  F    = s0**2
>
> does work, as does
>
>  un   = 0.5*dot(u0, n)
>  F    = 1.0/(s0**2 + 1.0)

Strange, I would have guess that the  'F' expression in the first example was the problem, but that it works in the last example is a mystery. I'll have a look at it.

Kristian

> Garth
>
> > --
> > Marie
> >
> > > Garth
> > >
> > >
> > >
> > > from dolfin import *
> > >
> > > mesh = UnitSquare(3, 3)
> > > n = FacetNormal(mesh)
> > >
> > > BDM = FunctionSpace(mesh, "Brezzi-Douglas-Marini", 1)
> > > P0 = FunctionSpace(mesh, "Discontinuous Lagrange", 0)
> > > mixed_space = MixedFunctionSpace([BDM, P0])
> > >
> > > # Function spaces and functions
> > > r  = TestFunction(P0)
> > > U0 = Function(mixed_space)
> > > u0, s0 = split(U0)
> > >
> > > U0.vector()[:] = 1.0
> > >
> > > un   = dot(u0, n)/2.0
> > > F    = 1.0/(s0**2 + 1.0)
> > > L = r*F*un*ds
> > >
> > > param0 = {"representation": "quadrature", "optimize": False}
> > > param1 = {"representation": "quadrature", "optimize": True}
> > > print "Vec norm (no opt):   ", assemble(L,
> > >                              form_compiler_parameters=param0).norm("l2")
> > > print "Vec norm (with opt): ", assemble(L,
> > >                              form_compiler_parameters=param1).norm("l2")
> > >
> > >
> > >
> > > > --
> > > > Marie
>
> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : ffc@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp

Attachment: signature.asc
Description: OpenPGP digital signature