Kristian Oelgaard wrote:
>
>
> On 1 February 2010 15:40, Garth N. Wells <gn...@cam.ac.uk> 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@lists.launchpad.net
>> Unsubscribe : https://launchpad.net/~ffc
>> More help : https://help.launchpad.net/ListHelp
>>
>
_______________________________________________
Mailing list: https://launchpad.net/~ffc
Post to : ffc@lists.launchpad.net
Unsubscribe : https://launchpad.net/~ffc
More help : https://help.launchpad.net/ListHelp
