dolfin team mailing list archive

[Kent Andre <kent-and@xxxxxxxxx>]

On on., 2009-08-19 at 09:23 +0100, Garth N. Wells wrote:
> 
> mspieg wrote:
> > 
> > On Aug 17, 2009, at 8:39 AM, Kent Andre wrote:
> > 
> >> On ma., 2009-08-17 at 13:14 +0100, Garth N. Wells wrote:
> >>>
> >>> Johan Hake wrote:
> >>>> On Monday 17 August 2009 14:03:55 Garth N. Wells wrote:
> >>>>> Johan Hake wrote:
> >>>>>> On Monday 17 August 2009 13:33:21 Garth N. Wells wrote:
> >>>>>>> Please send any comments on how a SubFunction should work as I'm 
> >>>>>>> going
> >>>>>>> to work on it soon (the Cahn-Hilliard demo runs in parallel(!) ;) 
> >>>>>>> with
> >>>>>>> the exception of the extraction of sub-Functions for output). I'm
> >>>>>>> inlined to change a sub function such that for
> >>>>>>>
> >>>>>>>    Function u0 = u[0]
> >>>>>>>
> >>>>>>> the Function u0 will only point the vector belonging to u rather than
> >>>>>>> creating a new vector. The advantages are:
> >>>>>>>
> >>>>>>> - No need to copy a vector
> >>>>>>> - DofMap will not require 'is_view'
> >>>>>>> - Possible to do things like
> >>>>>>>     u0  = 0.0;
> >>>>>>>     u1  += v;
> >>>>>>>     u0.interpolate();
> >>>>>> What would Function::vector() for a SubFunction then return? The 
> >>>>>> original
> >>>>>> full Vector?
> >>>>> Yes.
> >>>>
> >>>> Ok.
> >>>>
> >>>>>> It would be cool to add a view of the original Vector that only
> >>>>>> represents the values of the dofs in the SubFunction, without coping
> >>>>>> data. I fiddled around with this when adding slicing for the 
> >>>>>> PyDOLFIN la
> >>>>>> interface, but realized that it would be too difficult.
> >>>>> This is more or less what I plan to do, although internally. A user
> >>>>> wouldn't see the vector, but operations like interpolate would only
> >>>>> involve part of the vector.
> >>>>
> >>>> I see. Other direct Vector operations would then operate directly on 
> >>>> the 
> >>>> shared Vector, like get, set, aso.
> >>>>
> >>>>> We could add a class like
> >>>>>
> >>>>>    VectorView(GenericVector& x, DofMap& dof),
> >>>>>
> >>>>> which could derive from GenericVector, to provide views. It isn't a
> >>>>> priority for me though.
> >>>>
> >>>> Yes, I also thought along these lines, however I did not think of 
> >>>> doing it 
> >>>> using a DofMap, which really is the natural thing. I will also not have 
> >>>> possibility to priorities it for now.
> >>>>
> >>>
> >>> On second thought,
> >>>
> >>>    VectorView(GenericVector& x, std::vector<uint> map& map)
> >>>
> >>> would be better. The DofMap could produce the map,
> >>>
> >>>    std::vector<uint> map = dofmap.view();
> >>>
> >>> Garth
> >>>
> >>
> >> This would be great! This is sort of the reason why I implemented block
> >> matrices in the first place, but this is far more general and there
> >> is no need to mess with the block structure during assembly. 
> > 
> > Hi All, 
> >    just adding my me-too this suggestion (the Index Set (IS) in PETSc, 
> > provides similar functionality which is very useful)
> > 
> > In particular, I find that the new version of SubFunction  a bit 
> > misleading when it returns the entire vector rather than the sub vector 
> > corresponding to just the subfunction.
> > (I'm whining because it's  breaking some of my previous code, but that's 
> > less important).
> > 
> > Consider the following from a Stokes problem
> > 
> > Function w;  //full velocity-pressure solution
> > Function u = w[0]; 
> > Function p = w[1];
> > Function ux = u[0];
> > Function uy = u[1];
> > 
> > if you plot or write to a file  u,p,ux, or uy you get what you expect 
> > just the portion of the solution that is the subfunction,
> > 
> > but if you request something like
> > ux.vector().max() 
> > you'll get the unexpected result of the maximum of the entire vector i.e.
> > 
> > ux.vector().max() == uy.vector().max() == p.vector().max()
> > 
> > which in this case will be the maximum pressure.  This really threw me 
> > (particularly since I was trying to calculate the maximum vector 
> > magnitude and getting very strange solutions)
> > 
> 
> I've added a line for now to print a warning when extracting the vector 
> of a sub-Function. This should help stop anyone being caught out 
> expectantly by the change.
> 
> > I very much like the idea of not having all these copies floating 
> > around, but it seems like a subfunction should still retain its original 
> > consistent behavior of being restricted to just a portion of the 
> > original function.
> > 
> 
> A view should take care of this, but I'm not sure if we should have
> 
> a) A separate class VectorView with some limited functionality (get, 
> set, size, etc) which holds a pointer to the underlying vector object.
> 
> or
> 
> b) Provide an optional std::map to the constructor of the vector classes.
> 
> Option (a) is easy an not intrusive, but not easy to provide a wide 
> range of functionality. Computing norms would be tricky. Option (b) is 
> more work but would allow us to take advantage of functionality provided 
> by the linear algebra backend, like index sets in the case of PETSc.
> 
> Garth


How will you implement the VectorView ? With an offset, size and Vector
pointer? 

Kent