← Back to team overview

dolfin team mailing list archive

Re: Evaluate discrete function on cells

 

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?

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.

I have just one comment about getting the degrees of freedom from the Vector() associated with a DiscreteFunction. Let's say I am dealing with stress. Then, if I want to work with the full stress tensor my projection will look like in 3D

element = VectorElement(something)
DG = VectorElement("DG", "tet", 0, 9)

V = TestFunction(DG)
U = TrialFunction(DG)
f = Function(element)


a = dot(V,U)*dx
L = dot(V,f)*dx

Suppose I solve this projection for U. Then when I look at U.vector(), the degrees of freedom associated with a cell are not adjacent entries in the vector. Specifically for a UnitCube(1,1,1),

real* arr = new real[U.vector().size()];
U.vector().get(arr);

The stress components associated with a cell will look like arr(cell_nuber + offset*num_stress_components), where in this simple mesh offset = mesh.numCells(). For more general meshes I assume this would depend on the mesh, and the associated DofMap. I added hooks to grad the dof_map and dofs in the DiscreteFunction so I wouldn't have to think so hard. :)

Jake



References