← Back to team overview

dolfin team mailing list archive

  1. dolfin team
  2. Mailing list archive
  3. Message #06321

Re: Evaluating the FEM solution at an arbitrary point

  Thread Previous • Date Previous • • Thread Next 
Kristen Kaasbjerg wrote:

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 ?

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 ? - 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 ? - In DiscreteFunction.eval: a smart way to create the point (needed in gts) depending of the spacial dimension of the finite_element. 

Kristen


_______________________________________________
DOLFIN-dev mailing list
DOLFIN-dev@xxxxxxxxxx
http://www.fenics.org/mailman/listinfo/dolfin-dev
Hhhhmmm, something seems to have changed between version 0.7.1 and 0.7.2. When calling the eval function of class DiscreteFunction, the call to finite_element->evaluate_basis returns strange small numbers like 5.81e-268. Has anything changed between these two versions ? 
My implementation:
-------------------------------------------------------------------------------------------------------- 
void DiscreteFunction::eval(real* values, const real* x)
{
   uint n = finite_element->space_dimension();
   //step #1: locate the cell that contains x
   // check also if x has the correct dimension
   Point p(x[0],x[1]); // generalize to arbitrary dimension
   IntersectionDetector id(mesh);
   Array<uint> cells;
   id.overlap(p,cells);
   uint cell_index = cells[0];
   Cell cell(mesh,cell_index);
   UFCCell ufc_cell(cell);

   //step #2: get expansion coefficient on the cell
   real* coefficients = new real[n];
   this->interpolate(coefficients,ufc_cell,*finite_element);

   //step #3: multiply with basis functions on the cell
   real* basis_val = new real[finite_element->value_dimension(0)];
   real value(0.);
   for(uint i=0; i<n; i++)
   {
       finite_element->evaluate_basis(i,basis_val,x,ufc_cell);
       value += basis_val[0]*coefficients[i];
       cout << basis_val[0] << endl;
   }
   values[0] = value;
   delete [] coefficients;
   delete [] basis_val;
}
------------------------------------------------------------------------------------------------------------ 
I think I'll need some expertise to sort out this strange behavior.
I think it might be related to the python wrappings. At least, when
I make a C++ version of my demo, the correct values are returned
from the evaluate_basis function call. Whereas in python undefined
behavior is observed.

Kristen

References

  Thread Previous • Date Previous • • Thread Next