dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

-- 
Anders