dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Sat, Oct 04, 2008 at 11:16:45AM +0200, Martin Sandve Alnæs wrote:
> 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.
> >
> 
> Are cell() and facet() threadsafe? Imagine a function being
> evaluated in two threads iterating over different mesh partisions...

No, but we could fix it if necessary.

-- 
Anders

Attachment: signature.asc
Description: Digital signature