dolfin team mailing list archive

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

Johan Hake wrote:
> On Sunday 24 January 2010 16:12:37 Garth N. Wells wrote:
>> 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?
> 
> a_vec and v_vec are new vectors. None of the four vectors below get modified 
> by the a_vec and v_vec expressions so no need of copying, and the v0 and a0 
> assignment should work with GenericVectors too.
> 

Do you mean that just

  a_vec = 1.0/(2.0*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) \
           - (1.0-2.0*beta)*a0 )

where u and u0 are GenericVectors should work?

> Do you get any error messages?
>

No errors. What I have now seems to work fine.

Garth

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