ffc team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

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)

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