dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jan 29, 2008 at 07:31:16PM +0100, Dag Lindbo wrote:
> Anders Logg wrote:
> > 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?
> 
> Discrete vector function (think velocity components in Stokes 
> Taylor-Hood demo).
> 
> > 
> > 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.
> > 
> 
> Does this work for VectorElements? I feel like a newbie here (and DOLFIN 
> has definitely seen a lot of changes since last year).
> /Dag

Yes, just change to

  DG0 = VectorElement("DG", "your shape", 0)

And note that you can do all this in Python now (without calling FFC,
compiling, linking etc).

-- 
Anders