dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Feb 19, 2008 at 10:41:24AM +0100, Kristen Kaasbjerg wrote:
> Anders Logg wrote:
> > On Mon, Feb 18, 2008 at 02:54:38PM +0100, Kristen Kaasbjerg wrote:
> >   
> >> Anders Logg wrote:
> >>     
> >>> On Mon, Feb 18, 2008 at 02:29:23PM +0100, Kristen Kaasbjerg wrote:
> >>>   
> >>>       
> >>>> Anders Logg wrote:
> >>>>     
> >>>>         
> >>>>> On Mon, Feb 18, 2008 at 02:15:47PM +0100, Kristen Kaasbjerg wrote:
> >>>>>   
> >>>>>       
> >>>>>           
> >>>>>> Anders Logg wrote:
> >>>>>>     
> >>>>>>         
> >>>>>>             
> >>>>>>> Nice. The obvious thing would be to implement this in DiscreteFunction
> >>>>>>> and map that to the function call
> >>>>>>>
> >>>>>>>   virtual void Function::eval(real* values, const real* x) const;
> >>>>>>>
> >>>>>>> so that any Function (discrete, constant or user-defined) can be
> >>>>>>> evaluated at an arbitrary point.
> >>>>>>>
> >>>>>>> It should be possible to implement this for any kind of element, and
> >>>>>>> the code will look about the same as the code you have done for simple
> >>>>>>> elements.
> >>>>>>>
> >>>>>>> We might add some kind of caching so that evaluation at multiple
> >>>>>>> points that lie close to each other is efficient. (But maybe GTS is
> >>>>>>> smart and handles this already.)
> >>>>>>>
> >>>>>>>   
> >>>>>>>       
> >>>>>>>           
> >>>>>>>               
> >>>>>> ok guys, I have made a dirty hack in the C++ Function class
> >>>>>> in order to get the desired functionality. Looks very much like
> >>>>>> Dags code. Could I ask you to take a quick look at it (see below)
> >>>>>> to see if I have done anything alarming. So now both the cell searching 
> >>>>>> and the function evaluation can be done from python (and perhaps be 
> >>>>>> condensed into one function call if desired) and it
> >>>>>> seems to work.
> >>>>>> Thanks for your help along the way.
> >>>>>> Kristen
> >>>>>>
> >>>>>> -----------------------------------------------------------------------------------------------------------
> >>>>>> void Function::my_eval(real* values, const real* x,
> >>>>>>                        const ufc::cell& ufc_cell,
> >>>>>>                        const ufc::finite_element& finite_element,
> >>>>>>                        Cell& cell)
> >>>>>> {
> >>>>>>     if (!f)
> >>>>>>         error("Function contains no data.");
> >>>>>>     //step #1: get expansion coefficient on the cell
> >>>>>>     uint n = finite_element.space_dimension();
> >>>>>>     real* coefficients = new real[n];
> >>>>>>     this->interpolate(coefficients,ufc_cell,finite_element,cell);
> >>>>>>
> >>>>>>     //step #2: multiply with basis functions on the cell
> >>>>>>     real* basis_val = new real[finite_element.value_dimension(0)];
> >>>>>>     for(uint i=0; i<n; i++)
> >>>>>>     {
> >>>>>>         finite_element.evaluate_basis(i,basis_val,x,ufc_cell);
> >>>>>>         values[0] += basis_val[0]*coefficients[i];
> >>>>>>     }
> >>>>>> }
> >>>>>>     
> >>>>>>         
> >>>>>>             
> >>>>> Looks about right, but remember to delete the pointers coefficients
> >>>>> and basis_val.
> >>>>>
> >>>>> Extending this to non-simple elements should be fairly simple. Add
> >>>>> something like this:
> >>>>>
> >>>>>   // Compute size of value (number of entries in tensor value)
> >>>>>   uint size = 1;
> >>>>>   for (uint i = 0; i < finite_element->value_rank(); i++)
> >>>>>     size *= finite_element->value_dimension(i);
> >>>>>
> >>>>> Then iterate over the number of values for each basis function (size),
> >>>>> not only the first.
> >>>>>
> >>>>> Then we just need to include finding the element (using
> >>>>> IntersectionDetector) in this function, remove the arguments ufc_cell,
> >>>>> finite_element and cell, and then this can be added to DiscreteFunction.
> >>>>>
> >>>>>   
> >>>>>       
> >>>>>           
> >>>> Is the ufc::finite_element available in the Function class ?
> >>>> Kristen
> >>>>     
> >>>>         
> >>> No, but it's available in DiscreteFunction.
> >>>
> >>>   
> >>>       
> >> Ok, should I try to implement as much of this function as I can ?
> >>     
> >
> > Yes, that would be nice.
> >
> >   
> >> How are the Function and DiscreteFunction classes related and
> >> what type is the FEM solution you get out from dolfin ?
> >>     
> >
> > It's a so-called envelope-letter design (with a twist).
> >
> > Basically, Function acts as the front-end for users, but does
> > everything internally by calls to a pointer to a GenericFunction.
> > This pointer is instantiated to either a DiscreteFunction,
> > UserFunction or ConstantFunction depending on the arguments to the
> > constructor of Function.
> >
> > So when you call u.eval() for a Function, then you call
> > Function::eval(), which in turn calls GenericFunction::eval(), which
> > is overloaded by for example DiscreteFunction::eval() depending on the
> > representation of the function.
> >
> > In Function, you need to do something like
> >
> > void Function::eval(real* values, const real* x) const
> > {
> >   if (!f)
> >     error("Function contains no data.");
> >   f->eval(values, x);
> > }
> >
> > Then add eval() to the GenericFunction interface and implement eval()
> > in DiscreteFunction, UserFunction (should return the same error as in
> > Function now...) and ConstantFunction.
> >
> > See if you can find your way around...
> >
> >   
> Ok, have it implemented for simple elements now.
> Seems to work for all subclasses of GenericFunction.
> Should I bundle what I have and send it to you ?

Would be nice.

> I'm a little uncertain on the following things:
> - point on boundary between two or more elements - is it then enough to 
> do the evaluation in one of the elements ?

Yes, at least for now. We do the same thing now when postprocessing
for example DG solutions.

> - for UserFunction eval (which is only called when not overloaded) calls 
> the eval(x) of class Function (via f->eval(x)) to see if that then has 
> been overloaded, correct ?

Yes, but now when eval() is available for all function types, Function
needs to call f->eval() which should return an error for UserFunction
(since we should never reach there).

> - In DiscreteFunction.eval: a smart way to create the point (needed in 
> gts) depending of the spacial dimension of the finite_element.

ok. We'll figure this out.

-- 
Anders