dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

-- 
Anders