dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jan 29, 2008 at 08:33:01AM -0500, Jake Ostien wrote:
> Dag Lindbo wrote:
> > Hello all,
> >
> > I saw that there is a convenient way to evaluate a DiscreteFunction at 
> > all vertices of its mesh, i.e. u.evaluate( val ). Is there a convenient 
> > way to evaluate a function on each cell (e.g. on the mid-point of each 
> > triangle)?
> >
> > If not, I would really appreciate some suggestions how to go about. Is 
> > the more general function
> >
> > interpolate(real* coefficients,const ufc::cell& cell,
> > 	    const ufc::finite_element& finite_element)
> >
> > the way to go? How do I evaluate the sum of basis functions in a point 
> > if I know the coefficients?
> >   
> I have been playing with a similar functionality in my own project, 
> where I need to evaluate some field on a cell.  What I ended up doing is 
> adding two methods to the Function class, specifically for 
> DiscreteFunctions, that return the local data, uint* dofs (which are the 
> coefficients) and the DofMap* dof_map.  Then I can iterate over cells 
> and evaluate the coefficients on each cell however I want.  So my loop 
> looks something like:
>  
>   (say field is a DiscreteFunction) 
> 
>   for (; !cell.end(); ++cell) {
>     ufc_cell.update(*cell);
>     dof_map->tabulate_dofs(dofs, ufc_cell);
>    
>     field.vector().get(array, dof_map->local_dimension(), dofs);
> 
>     ...do something with array...
>   }
> 
> There is probably a better way to do this, and I'd like to hear it, 
> too.  Otherwise I can submit my own changes.
> 
> Jake

What type of function is it that you want to evaluate on cells?

Is it a finite element function? Defined by what element? Or is it a
user-defined function, something like sin(x)?

In any case, the best way to compute the values on cells would be to
just define a projection to piecewise constants. Then all values will
be in the Vector() that holds the degrees of freedom (if you should
need them).

Here's an example:

  element = FiniteElement(your space here)
  DG0 = FiniteElement("DG", "your shape", 0)
  v = TestFunction(DG0)
  u = TrialFunction(DG0)
  f = Function(element)

  a = v*u*dx
  L = v*f*dx

Solve the system and u will be the projection of f onto constants.

-- 
Anders