dolfin team mailing list archive

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

2008/7/7 Anders Logg <logg@xxxxxxxxx>:
> 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);

Ok! I wasn't sure which array class we're supposed to use in C++ dolfin...

> and in Python by
>
>  f = Function(mesh, [0.0, 0.0, 0.0])

Yes:

>> Also, it would be nice if the numpy typemaps used the right numpy
>> macros which converts any python sequence to the wanted format. Then
>> we could use tuples instead of array above.

We might get this for free with Array, since it's a std::vector and I
think it already has such typemaps.

--
Martin