ufl team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Monday 06 April 2009 16:47:27 Anders Logg wrote:
> On Mon, Apr 06, 2009 at 04:36:18PM +0200, Martin Sandve Alnæs wrote:
> > On Mon, Apr 6, 2009 at 4:24 PM, Anders Logg <logg@xxxxxxxxx> wrote:
> > > On Mon, Apr 06, 2009 at 03:40:03PM +0200, Martin Sandve Alnæs wrote:
> > >> On Mon, Apr 6, 2009 at 3:07 PM, Anders Logg <logg@xxxxxxxxx> wrote:
> > >> > On Mon, Apr 06, 2009 at 02:42:33PM +0200, Martin Sandve Aln=E6s 
wrote:
> > >> >> 2009/4/6 Johan Hake <hake@xxxxxxxxx>:
> > >> >> > On Monday 06 April 2009 13:44:14 Martin Sandve Aln=E6s wrote:
> > >> >> >> I'm not sure, I was just thinking about this...
> > >> >> >>
> > >> >> >> Say we have:
> > >> >> >>
> > >> >> >> =A0 V =3D VectorElement("CG", triangle, 1)
> > >> >> >> =A0 P =3D FiniteElement("CG", triangle, 1)
> > >> >> >> =A0 TH =3D V+P
> > >> >> >> =A0 up =3D Function(TH)
> > >> >> >>
> > >> >> >> then up.shape() =3D=3D (3,) and thus len(up) =3D=3D 3 and you
> > >> >> >> can do
> > >> >> >>
> > >> >> >> =A0 ux, uy, p =3D up # only considers the value shape
> > >> >> >> =A0 u =3D as_vector((ux, uy))
> > >> >> >>
> > >> >> >> while ufl.split works like this:
> > >> >> >>
> > >> >> >> =A0 u, p =3D split(up) # deals with the mixed element structure
> > >> >> >> =A0 len(u) =3D=3D 2
> > >> >> >>
> > >> >> >> We can only implement __iter__ to do one thing!
> > >> >> >
> > >> >> > A naive solution could be to overload the __len__ and __iter__
> > >> >> > function=
> > >> >
> > >> > s in
> > >> >
> > >> >> > ufl.Functions, so they return the subfunctions as they are used
> > >> >> > in spli=
> > >> >
> > >> > t()?
> > >> >
> > >> >>=20
> > >> >> Yes, that's the general idea, but then iteration over a Function
> > >> >> and another ufl expression won't work the same way. If that's what
> > >> >> we choose to do I will probably remove iteration over general
> > >> >> expressions, because this will be too confusing.
> > >> >
> > >> > If I do
> > >> >
> > >> >  TH = V + P
> > >> >  w = Function(TH)
> > >> >
> > >> > then I expect w to be vector-valued function of size d + 1.
> > >>
> > >> Yes, nobody disagrees with that.
> > >>
> > >> > So I suggest the split functionality to be a special call:
> > >> >
> > >> >  ux, uy, p = w
> > >>
> > >> Ok, just like I just implemented __iter__ in UFL.
> > >> Requires that __iter__ is _not_ overloaded for Function in dolfin.
> > >>
> > >> >  (ux, uy), p = w.split()
> > >>
> > >> This line is not as innocent as it looks...
> > >>
> > >> Consider these two lines:
> > >> A)  u, p = w.split()
> > >> B)  (ux, uy), p = w.split()
> > >>
> > >> We already have (A), and it works fine.

Yes, but you cannot use it to define mixed ufl.forms, as the mixing of 
ufl.split and dolfin.split is non-trivial. 

> > >> B is like A plus this line:
> > >> C)  ux, uy = u
> > >>
> > >> so B will _either_ need to use __iter__ on u for splitting u ->
> > >> (ux,uy), _or_ w.split() must already have returned a
> > >> tuple(tuple(ux,uy), p), in which case A will result in u being a tuple
> > >> (which nobody wants).
> > >
> > > Yes, B will need to use __iter__ for splitting u into ux an px.
> >
> > Then (ux, uy) won't have the same type as p, unless you solve the
> >
> > below problem that is:
> > >> Note that you didn't differentiate between ufl.split(w) and w.split(),
> > >> they are still not the same.
> > >
> > > They need to be merged somehow so that w.split() in Python does
> > > something sensible, that is, returns Functions which can be used as any
> > > other Functions (both in forms and for plotting).
> >
> > Good luck :)
>
> Thanks.

He, he who is going to implement this?

We have two discussions here. Should we use __iter__ for splitting Functions. 
We have not landed on this, but 

The cpp.Function is correctly splitted, and it gets instantiated with the 
correct FunctionSpace through the c++ interface. We can use this information 
to instantiate a _new_ ufl.Element and use this to instantiate the 
correspondning ufl.Function.

I imagine that this is not good enough, as the splitting is not done through 
the ufl.split, so it cannot be used to define mixed forms.

One could maybe modularise the split function so one can call it inside the 
constructors of the new Functions, maybe implement it as a generator 
expression? But my knowledge of how ufl treat this is zero. 

I see that if we just use the ufl.split on the result from stokes we get:

  type(u): <class 'ufl.tensors.ListTensor'>
  type(p): <class 'ufl.indexed.Indexed'>
  type(u).__bases__: (<class 'ufl.expr.WrapperType'>,)
  type(p).__bases__: (<class 'ufl.expr.WrapperType'>,)
  u: [[ w_4[0], w_4[1] ]]
  p: w_4[2]

The same result printed for u, v = w.split()

  type(u): <class 'dolfin.ufl.function.DiscreteSubFunction'>
  type(p): <class 'dolfin.ufl.function.DiscreteSubFunction'>
  type(u).__bases__: (<class 'ufl.function.Function'>,
                      <class 'dolfin.cpp.DiscreteFunction'>,
                      <class 'dolfin.ufl.function.Function'>)
  type(p).__bases__: (<class 'ufl.function.Function'>,
                      <class 'dolfin.cpp.DiscreteFunction'>,
                      <class 'dolfin.ufl.function.Function'>)
  u: <Function in <Function space of dimension 22272 (
        <Mixed element: (<CG2 on a <triangle of degree 1>>,
                         <CG2 on a <triangle of degree 1>>)>)>>
  p: <Function in <Function space of dimension 2868 (
        <CG1 on a <triangle of degree 1>>)>>

So ufl.split does not return two ufl.Functions, right? Only some sort of view? 
Is that the WrapperType? Are these classes inheritable and could they be 
instantiated inside an __init__ function? I also see (or not see!) that the 
information about the original function is hidden somewhere.

> > >> >> > I have __no__ clue what so ever to what other things this
> > >> >> > potentially c=
> > >> >
> > >> > ould
> > >> >
> > >> >> > break ;)
> > >> >>
> > >> >>=20
> > >> >> For some reason I had explicitly implemented __iter__ using
> > >> >> NotImplemented to signal it shouldn't be used, and as I said
> > >> >> I don't remember why that was:
> > >> >>=20
> > >> >>
> > >> >> >> It might be for ambiguities similar to this that I didn't
> > >> >> >> implement this earlier, but I don't remember any
> > >> >> >> more compelling reasons...
> > >> >>
> > >> >>=20
> > >> >>
> > >> >> >> In PyDOLFIN, Function.split() is something else?
> > >> >> >>
> > >> >> >> It creates SubFunctions, which can't be used in forms.
> > >> >> >> Right?
> > >> >> >
> > >> >> > As it is now you are actually able to use the returned
> > >> >> > subfunction to d=
> > >> >
> > >> > efine
> > >> >
> > >> >> > forms. It's a bug which I will fix :P.
> > >> >> >
> > >> >> > It turns out that the subfunctions has two version of the
> > >> >> > FunctionSpace=
> > >> >
> > >> > =2E One
> > >> >
> > >> >> > in python, which is just the original FunctionSpace that was used
> > >> >> > to de=
> > >> >
> > >> > fine
> > >> >
> > >> >> > the MixedFunctionSpace and one which is the correct
> > >> >> > FunctionSpace, incl=
> > >> >
> > >> > uding
> > >> >
> > >> >> > the correct offset, which is stored by the c++ version. Nice!
> > >> >> > Not...
> > >> >> >
> > >> >> > I will fix this so the versions correponds.
> > >> >>
> > >> >>=20
> > >> >> I won't even try to wrap my head around that paragraph,
> > >> >> this is why I'm continuing to use C++ while getting UFL done :-/
> > >> >> (Don't waste your time trying to explain it better,
> > >> >> I just don't want to know at this point)

That's ok, but then I have to understand some ufl internals, see question 
about ufl.split above.

> > >> >>=20
> > >> >>=20
> > >> >>
> > >> >> >> I think we should keep that explicit.
> > >> >> >
> > >> >> > While I see the ambiguity I implemented this in PyDOLFIN.
> > >> >> >
> > >> >> > Now in the stokes demo we can do:
> > >> >> >
> > >> >> > =A0u, p =3D problem.solve()
> > >> >> >
> > >> >> > This is what you thought of Anders? You desides what get pushed
> > >> >> > :)
> > >> >
> > >> > It looks good, but perhaps it should not be implemented using
> > >> > __iter__ but instead the solve() function checks if there it is
> > >> > solving a mixed system and in that case calls split() before
> > >> > returning. So that it sometimes returns a single solution and
> > >> > sometimes a tuple. Or might that also be confusing?
> > >>
> > >> We can't do that, because the user might want the solution function
> > >> in the mixed space for further use. There's no "un-split".
> > >
> > > Yes we can. We can do as in the assembler, that the assemble()
> > > function returns an object (which is a matrix, vector or number) but
> > > also have the possibility of providing an extra argument in which to
> > > put the solution:
> > >
> > >  problem.solve(solution=w)
> > >
> > > It's also the case that someone wanting more control and doing
> > > something "advanced" will assemble and solve manually anyway and have
> > > full control over when to split.
> >
> > Ok.

So the solve, function will return a the result of w.split(). If the user want 
to the whole solution it will use a kwarg?

> > >> >> But what does this mean?
> > >> >>   u, p =3D ufl.split(up)
> > >> >> or
> > >> >>   u, p =3D up.split()
> > >> >> ?
> > >> >>=20
> > >> >> Unless there is a solution to make these the same (which seems
> > >> >> complicated), I still think this should be done explicitly.
> > >> >>=20
> > >> >> Also, I've said before: PyDOLFIN.Function should _not_ implement
> > >> >> any special functions (__foo__) since they may interfere with UFL.
> > >> >> The "u, p =3D up" feature breaks with this.
> > >> >
> > >> > I'm not sure that's an absolute rule. There may be exceptions.
> > >> > One such is __repr__.
> > >>
> > >> Repr has a very specific meaning and should not be
> > >> implemented for dolfin.Function, because a proper
> > >> implementation of it would require a full string representation
> > >> of the mesh.
> > >>
> > >> Repr is used for comparing objects many places in
> > >> UFL and as such is particularly important to _not_
> > >> overload. I said this several times before, and nobody
> > >> argued then. Feel free to overload __str__.
> > >
> > > ok, good, __str__ then. So there are exceptions. There may be others.

Now we explicitly say that 

  dolfin.Function.__repr__ == ufl.Function.__repr__

and 
  
  dolfin.Function.__str__ returns "<Function in %s>" % str(self._V)

see output above.

> > Indeed there may be, but implementing any such should only
> > be done after carefully considering the implications in UFL.

Agree.

> > And that is not the way this is being done, because:
> > >> If you disagree with UFL design, please say so at ufl-dev.
> > >> Nobody has done so, and I can't read minds.
> > >
> > > This is about DOLFIN design: which functions to put in the DOLFIN
> > > Function class and of course we need to be careful when overloading
> > > UFL special functions.
> >
> > That's just the problem: you're treating this as just dolfin design
> > but it's not just that when it puts restrictions on the design of ufl.
> > That's the result of the multiple inheritance pattern used for Function.
>
> It's natural that the requirements in DOLFIN may influence some
> decision in UFL. But so far it's only __repr__ (which should be
> renamed __str__) that's overloaded 

Is this correct, refer to what I wrote above.

> and we'll be careful in choosing what to overload next.

Johan