dolfin team mailing list archive

[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

Attachment: signature.asc
Description: Digital signature