dolfin team mailing list archive

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

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.

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.

--
Martin