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 > Just my 2 cents...you guys do all the hard work > cheers > marc > > > > I'm running into some bugs with the new SubFunction design returning the > full vector in a PETSc Vector, > > I have code >> >> Kent >> >> _______________________________________________ >> DOLFIN-dev mailing list >> [email protected] <mailto:[email protected]> >> http://www.fenics.org/mailman/listinfo/dolfin-dev > > ---------------------------------------------------- > Marc Spiegelman > Lamont-Doherty Earth Observatory > Dept. of Applied Physics/Applied Math > Columbia University > http://www.ldeo.columbia.edu/~mspieg > tel: 845 704 2323 (SkypeIn) > ---------------------------------------------------- > > _______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
