dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Mon, Jul 07, 2008 at 11:28:40AM +0200, Martin Sandve Alnæs wrote:
> 2008/7/7 Garth N. Wells <gnw20@xxxxxxxxx>:
> >
> >
> > Anders Logg wrote:
> >> On Fri, Jul 04, 2008 at 11:25:52AM +0200, Martin Sandve Alnæs wrote:
> >>> 2008/7/3 Martin Sandve Alnæs <martinal@xxxxxxxxx>:
> >>>> 2008/7/3 Martin Sandve Alnæs <martinal@xxxxxxxxx>:
> >>>>> I've imple
> >>>>>
> >>>>> 2008/7/3 Dag Lindbo <dag@xxxxxxxxxx>:
> >>>>>> Martin Sandve Alnæs wrote:
> >>>>>>> 2008/7/2 Dag Lindbo <dag@xxxxxxxxxx>:
> >>>>>>>> Martin Sandve Alnæs wrote:
> >>>>>>>>> I just added checks for valid Functions to the assembler, but this
> >>>>>>>>> turned out to break some demos so I changed it to warnings.
> >>>>>>>> This seems like helpful checks to make!
> >>>>>>>>
> >>>>>>>>> The reason is that Function(mesh, 0.0), which creates a ConstantFunction,
> >>>>>>>>> is used for zero vector functions as well, which leaves the value shape
> >>>>>>>>> of such a function undefined. This makes error checking impossible,
> >>>>>>>>> which is not very nice.
> >>>>>>>>>
> >>>>>>>>> In my opinion, the correct next step is to keep these warnings,
> >>>>>>>>> and turn them into errors after "a while". We should either add
> >>>>>>>>> support for ConstantFunctions with tensor shape (would be nice),
> >>>>>>>>> or make the demos and apps use user-define functions.
> >>>>>>>> Adding tensor shape to the existing concept of ConstantFunction seems
> >>>>>>>> like the way to go. IMHO, the convenience if this over user defined
> >>>>>>>> functions is too great to discard.
> >>>>>>>>
> >>>>>>>> Would it not just amount to having ConstantFunction::rank() return the
> >>>>>>>> proper value as determined by the mutable member size, instead of 0? (on
> >>>>>>>> line 28 of ConstantFunction.cpp)
> >>>>>>> No, because this size isn't known at construction time, only when
> >>>>>>> interpolate is called with a finite_element object. A suitable (set
> >>>>>>> of) constructor(s) needs to be defined, mirrored in Function, and
> >>>>>>> handled in the swig interface for python.
> >>>>>> Yes, I see that this is needed if the rank and shape of the function is
> >>>>>> to be set constructor. Clearly, this makes for a nice user interface! My
> >>>>>> suggestion was more along the lines of a quick fix to get rid of the
> >>>>>> warnings: rank() returns zero if size == 1, returns 1 if size == 2 or 3 etc.
> >>>>>>
> >>>>>> /Dag
> >>>>>>
> >>>>>>> // Scalar constructor
> >>>>>>> Function f(mesh, 0.0);
> >>>>>>>
> >>>>>>> // General tensor constructor:
> >>>>>>> Function f(Mesh& mesh, uint rank, uint* shape, uint* values);
> >>>>>>>
> >>>>>>> // vector
> >>>>>>> uint rank = 1;
> >>>>>>> uint shape[1] = { 2 };
> >>>>>>> real values[2] = { 1., 2. };
> >>>>>>> Function f(mesh, rank, shape, values);
> >>>>>>> // or
> >>>>>>> uint size = 2;
> >>>>>>> real values[2] = { 1., 2. };
> >>>>>>> Function f(mesh, 1, &size, values);
> >>>>>>>
> >>>>>>> // 3x3 tensor:
> >>>>>>> uint rank = 2;
> >>>>>>> uint shape[2] = { 3, 3 };
> >>>>>>> real values[9] = {
> >>>>>>>    1., 2., 3.,
> >>>>>>>    4., 5., 6.,
> >>>>>>>    7., 8., 9. };
> >>>>>>> Function f(mesh, rank, shape, values);
> >>>>>>>
> >>>>>>
> >>>>> I just implemented this, but can't reproduce the old behaviour at the
> >>>>> same time, so the demos must be updated and I guess many user
> >>>>> applications as well. I'd like to have consensus to do this change
> >>>>> before I continue with the demos and check in.
> >>>>>
> >>>>> The typical example is a zero function intended to be a vector:
> >>>>>  Function f(mesh, 0.0);
> >>>>>
> >>>>> Which now would become:
> >>>>>  uint fsize = 3;
> >>>>>  real fvalues[3] = { 0.0, 0.0, 0.0 };
> >>>>>  Function f(mesh, 1, &fsize, fvalues);
> >>>>>
> >>>>> I could evt. add a vector constructor:
> >>>>>  real fvalues[3] = { 0.0, 0.0, 0.0 };
> >>>>>  Function f(mesh, 3, fvalues);
> >>>> And:
> >>>>
> >>>>  Function f(mesh, 3, 0.0);
> >>>>
> >>>> which does the same as the above.
> >>>>
> >>> I'd rather not wait a month before I merge, so I'll push these changes.
> >>> Somebody who knows the logic should look at dolfin_function_pre.i
> >>
> >> I think Ola knows these typemaps best.
> >>
> >>> and fix the typemaps there, until then the python syntax for constant
> >>> vector functions will be:
> >>>
> >>> f = Function(mesh, size, array)
> >>>
> >>> E.g.
> >>> f = Function(mesh, 3, numpy.array((1.,2.,3.)))
> >>>
> >>> Or simply:
> >>> f = Function(mesh, 3, 0.0)
> >>
> >> I think that it would be nicer to create a vector-valued function in
> >> C++ by
> >>
> >>   Array values(0.0, 0.0, 0.0);
> >>   Function f(mesh, values);
> >>
> >> and in Python by
> >>
> >>   f = Function(mesh, [0.0, 0.0, 0.0])
> >>
> >
> > Looks nice to me.
> >
> > Garth
> 
> I've done this and added
> 
> Array<uint> shape(2,2);
> Array<real> values(1.0, 2.0,
>                                 3.0, 4.0);
> Function f(mesh, shape, values);
> 
> as well.
> 
> I kept this notation:
>   Function f(mesh, 3, 0.0);
> since it's shorter than defining the Array first, and the size (3) can
> be a variable.

Very nice, but there are already 20 different ways or so one may
initialize a Function so at some point maybe we need to do a cleanup
or at least document the different options very carefully.

-- 
Anders



> In Python, the typemaps for Array doesn't seem to work,
> so at the moment the syntax there is now:
> 
> f = Function(mesh, 3, 0.0)
> f = Function(mesh, ArrayDouble(1.0, 2.0, 3.0))
> f = Function(mesh, ArrayUInt(2,2), ArrayDouble(1.0, 2.0, 3.0, 4.0))
> 
> Note that ArrayUInt(2) will give a length-two array with zeros,
> in contrast with how I used it in the commit message...
> This isn't a problem, since one of the first two lines above
> should be fine for vector functions.
>

Attachment: signature.asc
Description: Digital signature