dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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).

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.

-- 
Anders

Attachment: signature.asc
Description: Digital signature