← Back to team overview

dolfin team mailing list archive

Re: SubFunctions

 



Kent Andre wrote:
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?

With a map, at least for the creation of the view. It might be possible to use backend supplied tools for the view behind the scenes. Offsets don't work when the dof map is renumbered.

Garth


Kent





References