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 _______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
