dolfin team mailing list archive

["Matthew Knepley" <knepley@xxxxxxxxx>]

On Jan 29, 2008 7:39 AM, Anders Logg <logg@xxxxxxxxx> 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?
>
> 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.

Just so you are aware, you can probably get your original result
by reducing the DG0 space to 1 point quadrature.

   Matt

> --
> Anders
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which
their experiments lead.
-- Norbert Wiener