Subject: [UFC-dev] evaluate_dof

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 07 Oct 2009 17:56:09 +0100

In the UFC manual, there is nothing specific for the function

double ufc::finite_element::evaluate_dof(unsigned int i,
                                          const function& f,
                                          const cell& c) const = 0;

on whether or not f is defined on the cell c (as in it's an FE function 
on the mesh that contains c). The FFC generated code assumes that f is 
defined in terms of c.

This is causing some issues on the DOLFIN side when interpolating on 
non-matching meshes. Is it an idea to pass on optional argument to ufc 
functions which take a ufc::function to indicate whether or not f is 
defined on c?

Garth

---

From: Anders Logg <[email-masked]>
Date: Wed, 7 Oct 2009 22:22:12 +0200

> on the mesh that contains c). The FFC generated code assumes that f is

Does it? It just assumes that f can be evaluated at points that within
the cell c as far as I know.

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 07 Oct 2009 21:36:15 +0100

> the cell c as far as I know.

In

   demo/function/nonmatching-interpolation/cpp/P3.h

and inside the function

   double evaluate_dof(unsigned int i,
                               const ufc::function& f,
                               const ufc::cell& c) const

there is a call to ufc::function::evaluate(....),

   f.evaluate(values, coordinates, c);

The cell c is the cell for which the dof functional is being evaluated, 
but not necessarly the cell on which f lives. The UFC manual for
ufc::function::evaluate says that it evaluates the function at the point 
'coordinates' in the cell 'c'. This inconsistency makes the Function 
design tricky for non-matching meshes.

Garth

> 
> --
> Anders
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> UFC-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufc-dev

---

From: Anders Logg <[email-masked]>
Date: Wed, 7 Oct 2009 22:41:49 +0200

> design tricky for non-matching meshes.

I see the problem now. f has no way to know if c is one of its own
cells or a cell in another mesh so the cell information is basically
useless.

We might want to consider adding an extra piece of information to
ufc::cell which is a hash value for the mesh to which it belongs. Then
one can easily check if it's a cell in this mesh or another mesh.

I'll set up a blueprint for it.

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 07 Oct 2009 21:52:26 +0100

> cells or a cell in another mesh so the cell information is basically

Even worse, it will likely lead to errors in application code.

> one can easily check if it's a cell in this mesh or another mesh.

That would work, or ufc::finite_element::evaluate_dof could take bool 
which indicates whether or not the cells match. This would require an 
extra ufc::function member function

   evaluate(double* values, const double* coordinates)

The UFC code would call the right function and dolfin::Function wouldn't 
have to worry about it.

Garth

---

From: Anders Logg <[email-masked]>
Date: Wed, 7 Oct 2009 23:01:58 +0200

> The UFC code would call the right function and dolfin::Function wouldn't

I'm not sure that would work. It would be difficult to set such a flag
and get it right. You set the flag, then call some other function
which calls some other function. In the end, you may or may not end up
in the same mesh and it might be difficult to know when you first make
the call. But you do know the mesh when you create the ufc::cell (in
UFCCell) and you can compare the mesh hashes when it finally reaches
evaluate().

> > I'll set up a blueprint for it.

https://blueprints.launchpad.net/ufc/+spec/extend-ufc-data

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 07 Oct 2009 22:07:14 +0100

Anders Logg wrote:
> the call. 

This is a problem with design too many calls back and forth that ones 
doesn't know what's going on.

But you do know the mesh when you create the ufc::cell (in
> UFCCell) and you can compare the mesh hashes when it finally reaches

I find it misleading that code reaches an eval function with a ufc::cell 
argument which is meaningless.

---

From: Anders Logg <[email-masked]>
Date: Wed, 7 Oct 2009 23:11:21 +0200

> doesn't know what's going on.

There's just one callback in the UFC interface which is function::evaluate.
I don't see any way to avoid it.

But I agree it can be confusing. At least the Function classes are
less confusing now than they used to be. Another reason for the many
redirections is reuse of code. We could remove some of them copying
code which is worse.

> argument which is meaningless.

Yes, I agree. The current ufc::cell argument in function::evaluate is
meaningless, but it won't be if we add a hash/id for the mesh.

--
Anders

---

From: [email-masked]
Date: Thu, 8 Oct 2009 08:55:06 +0200

> The UFC code would call the right function and dolfin::Function wouldn't

I think adding this function is a good solution.
Evaluating the dof on a cell on which the finite element is not easy
and should be avoided.

However, evaluating the basis functions in points should always be defined
(at least with sufficient continuity). And this is typically what is needed.
There should be no need to use evaluate_dofs when evaluating a function.

Kent

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 08 Oct 2009 08:39:14 +0100

> There should be no need to use evaluate_dofs when evaluating a function.

We do actually want to use evaluate_dofs because dofs are what we're 
after on the given cell. It's just the function f from which we compute 
the dofs is not defined on the given cell. What we want to do is 
interpolate f on a different mesh, for which we needs the dof 
functionals on a given cell.

Garth

---

From: Kent Andre <[email-masked]>
Date: Thu, 08 Oct 2009 09:58:09 +0200

> functionals on a given cell.

I thought you had a 

f = \sum_i f_i N_i 

(f_i given or unknown)

Then 

f(x_j) = \sum f_i N_i(x_j) 

Here it is not important to know how the fem dofs are defined, ie. if 
it is a normal component on a facet or something. Hence, evaluate_dof 
should not be called in this application. 

It is only important to know the fem-dofs when we assign Dirichlet
conditions. And I guess you will assign them on the right mesh ? 

Kent

---

From: Anders Logg <[email-masked]>
Date: Thu, 8 Oct 2009 13:47:39 +0200

> conditions. And I guess you will assign them on the right mesh ?

We need to evaluate the dofs when we want to compute the interpolation
of f into another space:

  \pi_h : V \rightarrow W
  f \mapsto \pi_h f

Then we need to evaluate the dofs of W for f. This results in a call
to evaluate_dof, which in turn will result in one or more calls to
evaluation of f, and that evaluation is at points in a cell of the
mesh for W, not for V.

If V and W have the same mesh, then the extra ufc::cell argument to
ufc::function::evaluate can be used to speed up the evaluation, but if
not we first need to locate the cell to which the point belongs. There
is currently no way to see the difference between the two cases.

--
Anders

---

From: Kent Andre <[email-masked]>
Date: Thu, 08 Oct 2009 14:33:19 +0200

> is currently no way to see the difference between the two cases.

ok, I think I see your point. There are two ways to 
evaluate f in W. 

Either we can evaluate v \in V as
f(x_j) = \sum f_i N_i(x_j)
in a set of points x_j defined by W. 

Points defined by e.g. one element in W may be in several different
elements in V. Hence, we need a search, possibly for each point, 
but there is no need to evaluate the dof. Instead a look-up is
needed in the coefficient vector.  
 
The other approach, as you describe here, is to evaluate the dofs on W. 

In both cases I guess an additional cell argument may speed-up the
computation. 

The pi_h mapping takes V to the nodal points representation in W 
while the other mapping takes f \in V to 'f evaluated at quadrature
points in W'.  

Kent

---

From: Anders Logg <[email-masked]>
Date: Thu, 8 Oct 2009 14:40:49 +0200

> The other approach, as you describe here, is to evaluate the dofs on W.

Evaluating the dofs of W for f will lead to point evaluation of f as
in your example above.

> while the other mapping takes f \in V to 'f evaluated at quadrature

I don't see a different between these two. In both cases we need to
evaluate f at an arbitrary point.

The difference is between the following three cases:

1. Interpolation / restriction to a cell in the same space. In this
case we can just pick the dofs from the vector.

2. Interpolation / restriction to a cell in the same mesh but with a
different function space (say interpolation from P2 to P1). In this
case we need to evaluate the linear combination but we know the cell
so we can just grab the dofs from the vector and evaluate the basis
functions.

3. Interpolateion / restriction to a cell of a different mesh. In this
case we need to first locate the cell for the given point before we
can proceed as in (2).

The Function class currently handles this by first differentiating
between (1) and (2) in Function::restrict(). The problem is
differentiating between (2) and (3) which is currently not possible
since we can't check if a ufc::cell comes from the same mesh.

--
Anders