ufl team mailing list archive

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

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.
>> 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 :)


>> >> > 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)
>> >>=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.

>> >> 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.

Indeed there may be, but implementing any such should only
be done after carefully considering the implications in UFL.
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.

Martin