ffc team mailing list archive

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

Kristian Oelgaard wrote:
> 
> 
> 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.
> 

Thanks!

Garth

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