← Back to team overview

dolfin team mailing list archive

Re: [HG DOLFIN] Automatically interpolate user-defined functions on assignment

 



Martin Sandve Alnæs wrote:
On Thu, Mar 12, 2009 at 11:05 PM, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
Anders Logg wrote:
On Thu, Mar 12, 2009 at 05:41:42PM +0100, Martin Sandve Alnæs wrote:
On Thu, Mar 12, 2009 at 5:26 PM, Anders Logg <logg@xxxxxxxxx> wrote:
On Thu, Mar 12, 2009 at 04:59:43PM +0100, Martin Sandve Alnæs wrote:
On Thu, Mar 12, 2009 at 4:33 PM, Anders Logg <logg@xxxxxxxxx> wrote:
On Thu, Mar 12, 2009 at 01:42:55PM +0100, Martin Sandve Alnæs wrote:
On Thu, Mar 12, 2009 at 1:16 PM, Johan Hake <hake@xxxxxxxxx> wrote:
On Thursday 12 March 2009 12:58:33 Garth N. Wells wrote:
Anders Logg wrote:
On Thu, Mar 12, 2009 at 10:46:14AM +0100, Martin Sandve Alnæs wrote:
On Thu, Mar 12, 2009 at 10:32 AM, Anders Logg <logg@xxxxxxxxx> wrote:
On Thu, Mar 12, 2009 at 10:25:36AM +0100, Martin Sandve Alnæs wrote:
I have serious problems with the idea of letting user-defined
functions change nature to discrete functions as a side effect
of other operations, effectively hiding the user-defined
implementation of eval.

(I know this concept wasn't introduced in this discussion, but it's
related).
You mean by calling u.vector()? Yes, I agree it's problematic.

But I don't know how to handle it otherwise. Consider the following
example:

 u = Function(V)
 solve(A, u.vector(), b)

First u is a user-defined function since it doesn't have a
vector. Then it becomes discrete when we ask for the vector.

The problem I think is that there is no way for us to check whether a
user has overloaded eval() without trying to call it.
If we didn't require that user-defined functions had a function space
for various operations, this wouldn't be as large a problem.
Then you could do

class MyFunction(Function)
{
  MyFunction(V): Function(V) {}
  MyFunction(): Function() {}
  void eval(...) { ... }
}

MyFunction f;
f.vector(); // error, no function space!

MyFunction f(V);
f.vector(); // ok, should interpolate and make f discrete

This is already possible, but function spaces get attached
unneccesarily:


1) It's not necessary to attach function spaces to functions for
assembly, since a user-defined function is evaluated through
ufc::finite_element::evaluate_dofs for each cell anyway.
This would remove the need for the side-effect in

MyFunction f;
MyForm a;
a.f = f; // attaches function space to f
Isn't the function space used to deduct the value rank during assemble now? If
we go for this approach we probably need to readd the rank function and make
it virtual?
I don't see how the current approach with rank deducted from
some function space is any safer than deducting it from the
finite element of the form. Which is to say: this is unsafe, the
user should define the rank when defining eval.
We just removed the rank() and dim() arguments from the Function
interface as they are not needed.
Are not needed because a Function always has a FunctionSpace
or for some other reason? If we don't require a FunctionSpace
I think they should be re-added. I expect to be able to use a Function
for other things than assembly, e.g., if I get some arbitrary Function
to some code:

void foo(Function & f)
{
   ... compute something with f
}

I can't know the dimensions of arguments to f.
This means I can't write generic code to do something with Functions.
One of the simplifying assumptions we made a while back when
redesigning the Function classes was that a Function always has a
FunctionSpace. This still seems to be a good assumption and perhaps we
should enforce it even stronger than we do now.

It would also make it easy to check for errors during assembly. The
form knows the correct FunctionSpace and the assembler may compare
that FunctionSpace against the FunctionSpace of the Function.
I'm happy with "a Function always has a FunctionSpace" concept if it is
strictly enforced at construction time. The current situation is not so,
and since I'm allowed to create a user function without a function space
I expect to be able to use it where it makes sense.

C++ is often much easier to work with if objects are completely initialized at
construction time, then there are a lot of mistakes you're just not
allowed to make.
I agree completely. This is the approach we tried to follow in the
recent redesign of Function. That's why one first needs to create a
mesh, then a function space, then a form. I'm not sure why we
decided not to enforce the requirement of having a function space for
coefficients.

It seems it wouldn't be anymore difficult than what we have now. For
example, the Poisson demo would be

  UnitSquare mesh(32, 32);
  PoissonFunctionSpace V(mesh);

  PoissonBilinearForm a(V, V);
  PoissonLinearForm L(V);
  Source f(V);
  L.f = f; // without side effect

That's OK, but if you have >5 coefficients with >5 different functions
spaces, it's not nice.

This feature will without doubt lead to many users writing code where
function spaces are duplicated. If you have a form where more than one
coefficient shares same function space, you shouldn't use this feature.
Even if it is well documented, people don't generally read documentation
until things break.

If we strictly require that functions have function spaces, we can instead
check that the form coefficient function spaces match the given function
in "L.f = f;". Then your problem is reduced to runtime "typechecking".


The problem is the number of lines a user must program, and having intelligible names for the spaces.


Tests are hard, I know... But here I was thinking about simple
things like matching dimensions and vector sizes etc.
Yes, and I agree it's important. The tests are missing simply because we've
either been lazy or been pressed on time (and most often the latter).

I think that this point is being exaggerated. We can and should add more
tests, but there are tests. DOLFIN is not devoid of error checking.

True, I'm probably exaggerating, I just had some bad experiences.

Speaking of which, nobody answered my email about interpolate
three months ago. Many of those problems are still there.
Email is not a good bugtracker :-/

Should I just add unit tests that fail when I encounter some problem?


This is usually the most effective way to get things fixed - a failing unit test or demo ;), so go ahead and add it.

Garth


Martin




Follow ups

References