dolfin team mailing list archive

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

Johan Hake wrote:
> On Sunday 24 January 2010 00:03:41 Garth N. Wells wrote:
>> Johan Hake wrote:
>>> On Saturday 23 January 2010 14:55:08 Garth N. Wells wrote:
>>>> Johan Hake wrote:
>>>>> On Saturday 23 January 2010 08:42:14 Garth N. Wells wrote:
>>>>>> Is it correct that behind the scenes that
>>>>>>
>>>>>>   U0 = Function(V)
>>>>>>   U  = Function(V)
>>>>>>   U0.vector()[:] = U.vector()[:]
>>>>>>
>>>>>> involves a GenericVector::get(..) call and a GenericVector::set(..)
>>>>>> call? If so, it isn't ideal since it introduces unnecessary new/delete
>>>>>> operations and unnecessary copying of data.
>>>>> None of GenericVector::get(..) or GenericVector::set(..) are invoked,
>>>>> see __getslice__ and __setslice__ in la_post.i.
>>>>>
>>>>>    U0.vector()[:]
>>>>>
>>>>> involves
>>>>>
>>>>>   GenericVector::operator =(..)
>>>>>
>>>>> and
>>>>>
>>>>>   U.vector()[:]
>>>>>
>>>>> involves
>>>>>
>>>>>   GenericVector::copy()
>>>>>
>>>>> However the latter is unnecessary as you instead can do:
>>>>>
>>>>>   U0.vector()[:] = U.vector()
>>>>>
>>>>> invoking the assignment operator of U0's vector with U's vector.
>>>> What happens if I do
>>>>
>>>>   x = U.vector()[:]
>>> It just triggers the copy method of GenericVector, which is the same
>>> behavior as for other itterable Python types.
>>>
>>>> ? Is x a numpy array?
>>> No you need to call array() to accomplish that.
>> OK. What I'm trying to do is
>>
>>     # Get vectors
>>     u_vec  = u.vector()[:]
>>     u0_vec = u0.vector()[:]
>>     v0_vec = v0.vector()[:]
>>     a0_vec = a0.vector()[:]
> 
> You should not need to make a copy of the vectors here.
>

How can I avoid it?

What I have now (which is working) is

    # Get vectors (copies of)
    u_vec  = u.vector().array()
    u0_vec = u0.vector().array()
    v0_vec = v0.vector().array()
    a0_vec = a0.vector().array()

    # Update acceleration and velocity
    a_vec = (1.0/(2.0*beta))*( (u_vec - u0_vec - v0_vec*dt)/(0.5*dt*dt)\
          - (1.0-2.0*beta)*a0_vec )
    v_vec = dt*((1.0-gamma)*a0_vec + gamma*a_vec) + v0_vec

Garth

>>     a_vec = 1.0/(2.0*beta)*((u_vec-u0_vec - v0_vec*dt)/(0.5*dt*dt) \
>>          - (1.0-2.0*beta)*a0_vec)
>>     v_vec = dt*((1.0-gamma)*a0_vec + gamma*a_vec) + v0_vec
> 
> The expression cause a lot of intermediate vectors, but should work, if beta, 
> gamma and dt are all scalars.
> 
>>     # Update (U0 <- U0)
>>     v0.vector()[:] = v_vec
>>     a0.vector()[:] = a_vec
> 
> This should work.
> 
> Johan
> 
>> I guess I need to call array() to make this work?
>>
>> Garth
>>
>>> Johan
>>>
>>>> Garth
>>>>
>>>>> Johan