ufl team mailing list archive

[Martin Sandve Alnæs <martinal@xxxxxxxxx>]

Sorry, I don't have time for this. Sorry for giving the impression
that I could participate in a discussion like this.

Martin


On Mon, Apr 6, 2009 at 9:37 PM, Johan Hake <hake@xxxxxxxxx> wrote:
> 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
>