dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Tuesday 04 November 2008 14:11:55 Anders Logg wrote:
> On Tue, Nov 04, 2008 at 01:59:33PM +0100, Johan Hake wrote:
> > On Tuesday 04 November 2008 13:15:51 Anders Logg wrote:
> > > On Mon, Nov 03, 2008 at 10:38:40PM +0100, Anders Logg wrote:
> > > > On Mon, Nov 03, 2008 at 06:25:18PM +0100, Martin Sandve Alnæs wrote:
> > > > > 2008/11/3 Anders Logg <logg@xxxxxxxxx>:
> > > > > > On Mon, Nov 03, 2008 at 03:55:31PM +0000, Garth N. Wells wrote:
> > > > > >> Anders Logg wrote:
> > > > > >> > On Mon, Nov 03, 2008 at 02:38:36PM +0000, Garth N. Wells wrote:
> > > > > >> >> Anders Logg wrote:
> > > > > >> >>> On Mon, Nov 03, 2008 at 11:22:26AM +0000, Garth N. Wells 
wrote:
> > > > > >> >>>> Anders Logg wrote:
> > > > > >> >>>>> On Sun, Nov 02, 2008 at 06:29:25PM +0000, Garth N. Wells 
wrote:
> > > > > >> >>>>>> Anders Logg wrote:
> > > > > >> >>>>>>> On Sun, Nov 02, 2008 at 05:52:21PM +0000, Garth N. Wells
> >
> > wrote:
> > > > > >> >>>>>>>> Do we want to insist that Dirichlet bc functions that
> > > > > >> >>>>>>>> do not appear inside a form are constructed with a
> > > > > >> >>>>>>>> FunctionSpace? DirichletBC is supplied with a
> > > > > >> >>>>>>>> FunctionSpace, so if the bc Function does not have a
> > > > > >> >>>>>>>> FunctionSpace, we could attach one automatically.
> > > > > >> >>>>>>>>
> > > > > >> >>>>>>>> Garth
> > > > > >> >>>>>>>> _______________________________________________
> > > > > >> >>>>>>>> DOLFIN-dev mailing list
> > > > > >> >>>>>>>> DOLFIN-dev@xxxxxxxxxx
> > > > > >> >>>>>>>> http://www.fenics.org/mailman/listinfo/dolfin-dev
> > > > > >> >>>>>>>
> > > > > >> >>>>>>> I think this is already handled. Look in the Poisson
> > > > > >> >>>>>>> demo. It uses a Constant to set the BC and it does not
> > > > > >> >>>>>>> have a FunctionSpace attached to it. The DirichletBC
> > > > > >> >>>>>>> class now uses its own FunctionSpace rather than the one
> > > > > >> >>>>>>> that the Function has (if any). There is a check (in
> > > > > >> >>>>>>> DirichletBC::check()) that checks that the FunctionSpace
> > > > > >> >>>>>>> for the Function is the same as the one in the
> > > > > >> >>>>>>> DirichletBC.
> > > > > >> >>>>>>
> > > > > >> >>>>>> It works for Constant, but not for Functions. I was
> > > > > >> >>>>>> getting an error when Function::interpolate is called.
> > > > > >> >>>>>> Function::interpolate leads to eval being called, in
> > > > > >> >>>>>> which case there is a test for the FunctionSpace which
> > > > > >> >>>>>> fails. Constant provides its own eval and therefore
> > > > > >> >>>>>> doesn't have a problem.
> > > > > >> >>>>>>
> > > > > >> >>>>>> For now, I've added a test in DirichletBC for the
> > > > > >> >>>>>> FunctionSpace. What we can add is an attach function if
> > > > > >> >>>>>> there is no FunctionSpace associated.
> > > > > >> >>>>>>
> > > > > >> >>>>>> Garth
> > > > > >> >>>>>
> > > > > >> >>>>> In which demo does this show up? Is there a simple way I
> > > > > >> >>>>> can comment something out to reproduce the error so I
> > > > > >> >>>>> understand what goes wrong?
> > > > > >> >>>>
> > > > > >> >>>> Look at /demo/nls/nonlinearpoisson/cpp.
> > > > > >> >>>>
> > > > > >> >>>> If you change
> > > > > >> >>>>
> > > > > >> >>>>    DirichletBoundaryCondition g(V, t);
> > > > > >> >>>>
> > > > > >> >>>> to
> > > > > >> >>>>
> > > > > >> >>>>    DirichletBoundaryCondition g(t);
> > > > > >> >>>>
> > > > > >> >>>> it will break down.
> > > > > >> >>>>
> > > > > >> >>>> Garth
> > > > > >> >>>
> > > > > >> >>> ok I see the problem now.
> > > > > >> >>>
> > > > > >> >>> The problem is a user may choose to either overload a scalar
> > > > > >> >>> eval function or a tensor eval function and we need to
> > > > > >> >>> decide which one after the callback from
> > > > > >> >>> ufc::function::evaluate(). If the FunctionSpace is not
> > > > > >> >>> known, we can't decide which one to pick.
> > > > > >> >>>
> > > > > >> >>> If we insist that one should be able to pass a Function
> > > > > >> >>> without a FunctionSpace to a DirichletBC, then we must
> > > > > >> >>> remove the scalar eval function.
> > > > > >> >>
> > > > > >> >> Fine with me. I think that it makes things simpler because
> > > > > >> >> the eval interface remains the same for all user-defined
> > > > > >> >> functions.
> > > > > >> >>
> > > > > >> >> Garth
> > > > > >> >
> > > > > >> > ok. It will also look the same as in Python.
> > > > > >>
> > > > > >> We discussed recently passing an object to eval() which contains
> > > > > >> some data. It would be useful the object also carried
> > > > > >> information on the rank and dimension of the function to allow
> > > > > >> checks and switching between 1D/2D/3D problems.
> > > > > >>
> > > > > >> Garth
> > > > > >
> > > > > > Yes, this would be nice, but it won't work as long as we don't
> > > > > > require that a Function always has a FunctionSpace. We could add
> > > > > > some new classes to handle error checking and data for
> > > > > > user-defined Functions:
> > > > > >
> > > > > >  void eval(Values& values, Data& data)
> > > > > >  {
> > > > > >    values[0] = sin(data.x[0]);
> > > > > >  }
> > > > > >
> > > > > > The class Values could check that data is not assigned to any
> > > > > > illegal indices and it could also check that all values have been
> > > > > > assigned etc, but if the Function does not know its
> > > > > > FunctionSpace, this can't be done.
> > > > > >
> > > > > > Another complication related to this but also to the
> > > > > > thread-safety of cell() and facet() is that UFC gets in the
> > > > > > middle of the call sequence:
> > > > > >
> > > > > >  assemble()
> > > > > >
> > > > > >  |--> Function::interpolate()
> > > > > >  |
> > > > > >       |--> FiniteElement::evaluate_dof()
> > > > > >       |
> > > > > >            |--> ufc::finite_element::evaluate_dof
> > > > > >            |
> > > > > >                 |--> ufc::function::evaluate()
> > > > > >                 |
> > > > > >                      |--> Function::eval()
> > > > > >
> > > > > > Since eval() is called from the generated UFC code, any arguments
> > > > > > like cell and facet passed from the assembler will be lost on the
> > > > > > way.
> > > > > >
> > > > > > Should we extend the UFC interface to allow sending a void* to
> > > > > > evaluate_dof which it will propagate to evaluate()?
> > > > >
> > > > > That's possible, but not very safe.
> > > > >
> > > > > A typesafe alternative is that dolfin::Function doesn't inherit
> > > > > ufc::function, but a ufc::function subclass is created which has a
> > > > > dolfin::Function which it calls and a dolfin::Data & which it sends
> > > > > in the call. Then we avoid ufc complications in dolfin::Function,
> > > > > and one such wrapper function can be created for each thread.
> > > > > It doesn't really add any more function calls, since it replaces
> > > > > the existing ufc::function::evaluate in the call stack you wrote
> > > > > above.
> > > >
> > > > Sounds good. I'll try something like this.
> > >
> > > This is now implemented. Take a look in Function::interpolate(),
> >
> > Should it be called FunctionData instead?
>
> Data is shorter and means less to write in the eval() functions (so we
> can make these look as simple as possible).

  "Make everything as simple as possible, but not simpler."

I think writing Function in addition to Data wont hurt.

> We could possibly make it a nested class of Function (Function::Data)
> but that might require that one writes Function::Data in the
> overloaded eval functions and it might be troublesome for SWIG.
>
> > As Data is implemented now it cannot be changed after creation, i.e., you
> > cannot set _facet or _cell. Would it be benefitial to cache the Data
> > object and the dofs, which both are created for each call of interpolate?
>
> Creating it each time is exactly what makes interpolate() thread-safe
> now which it wasn't before.

Ahh...

Johan