dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

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 get any error messages?

Johan

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