dolfin team mailing list archive

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

2008/10/4 Anders Logg <logg@xxxxxxxxx>:
> On Fri, Oct 03, 2008 at 05:39:28PM -0500, Matthew Knepley wrote:
>> On Fri, Oct 3, 2008 at 5:22 PM, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
>> >
>> >
>> > Anders Logg wrote:
>> >> On Fri, Oct 03, 2008 at 11:06:46PM +0100, Garth N. Wells wrote:
>> >>>
>> >>> Anders Logg wrote:
>> >>>> User-defined functions are currently defined by overloading eval().
>> >>>> The eval() function is available in two different versions, for vectors
>> >>>> (tensors) and scalars:
>> >>>>
>> >>>>   void eval(double* values, const double* x) const;
>> >>>>
>> >>>>   double eval(const double* x) const;
>> >>>>
>> >>>> This might be ok, but functions may also depend on the current cell,
>> >>>> the current facet, and perhaps on the current time and other data.
>> >>>> These are accessible by this->cell, this->facet.
>> >>>>
>> >>>> So, a function depends on x and other things and it's not obvious what
>> >>>> additional data a function does depend on. Are there any suggestions
>> >>>> for what the interface should look like for user-defined functions?
>> >>>>
>> >>>> One option could be to let functions be defined by expressions
>> >>>> (function pointers) given to the constructor of Function. That way,
>> >>>> one could define a function by
>> >>>>
>> >>>>   void source(const double* x) { sin(x[0]); }
>> >>>>
>> >>>>   f = Function(mesh, source);
>> >>>>
>> >>>> This way, there could be many different types of function pointers,
>> >>>> all taking a different number of arguments, and we could differentiate
>> >>>> between the different types.
>> >>>>
>> >>> I'm not a fan of function pointers, and I don't see how the above
>> >>> solution would help if a Function depends on other data.
>> >>
>> >> It would help to do things like
>> >>
>> >>   double foo(const double* x) { return sin(x[0]); }
>> >>   double bar(const double* x, real t) { return t*sin(x[0]); }
>> >>
>> >>   f = Function(foo);
>> >>   g = Function(bar);
>> >>
>> >
>> >
>> > I don't follow how this would work - Function will call foo, but it
>> > doesn't know what t is? PETSc uses pointers to function in SNES, but as
>> > far as I recall the function must have a prescribed interface.
>>
>> It has a context argument for user data.
>>
>>   Matt
>
> I remember now that's exactly why we don't use function pointers.
> We use functors (subclassing and overloading a callback function)
> instead. That allows users to add any data they want to their class
> and it will always be accessible to the user-overloaded function.
>
> Our problem here is different. Basically a user-defined function in
> DOLFIN may depend on a number of things (all provided by DOLFIN), so
> in general the function interface might look like
>
>  void eval(double* values,
>            const double* x,
>            const Cell& cell,
>            uint local_facet_number,
>            double t,
>            ... maybe other stuff);
>
> But often a user-defined function is scalar-valued and depends just on
> the point x, and then it's more convenient to just provide
>
>  double eval(const double* x) { return sin(x); }
>
> without having to list all the unused arguments.
>
> My suggestion is that we keep the cell() and facet() functions as now,
> which means these remain "hidden" function arguments. (They are mostly
> used by special functions such as MeshSize, FacetNormal etc anyway.)
>
> But then we provide 4 different versions of eval(), and the user may
> choose which one fits best:
>
>    double eval(const double* x);
>    double eval(const double* x, double t);
>
> --> void eval(double* values, const double* x);
>    void eval(double* values, const double* x, double t);
>
> All calls during assembly will go to
>
>    void eval(double* values, const double* x);
>
> which directs calls to the other 3.
>
> I think I can make this work. I'll put this into NewFunction later
> today and you may have a look.
>
> --
> Anders

Are cell() and facet() threadsafe? Imagine a function being
evaluated in two threads iterating over different mesh partisions...

--
Martin