← Back to team overview

dolfin team mailing list archive

Re: Function and DofMap

 



Anders Logg wrote:
On Tue, Sep 09, 2008 at 04:31:41PM +0200, Martin Sandve Alnæs wrote:
2008/9/9 Garth N. Wells <gnw20@xxxxxxxxx>:

Anders Logg wrote:
On Tue, Sep 09, 2008 at 11:30:26AM +0100, Garth N. Wells wrote:
Anders Logg wrote:
On Tue, Sep 09, 2008 at 09:06:00AM +0200, Anders Logg wrote:
On Mon, Sep 08, 2008 at 11:53:09PM +0200, Johan Hake wrote:
On Monday 08 September 2008 21:45:27 Martin Sandve Alnæs wrote:
2008/9/8 Anders Logg <logg@xxxxxxxxx>:
On Mon, Sep 08, 2008 at 11:12:14AM +0200, Martin Sandve Alnæs wrote:
2008/9/8 Johan Hoffman <jhoffman@xxxxxxxxxx>:
2008/9/8 Dag Lindbo <dag@xxxxxxxxxx>:
Anders Logg wrote:
There seems to be a problem (among many) with the current design of
the Function classes (see thread "evaluating higher order mesh
function").

In particular, the finite element is missing in DiscreteFunction.
My suggestion would be to just add it and let a DiscreteFunction
consist of the following four items which are always available:

  mesh, x, dof_map, finite_element

Is this enough, and what other issues to we need to fix?
I'm not sure I agree that the dof map and finite element should be
owned by the discrete function. There was a great suggestion from
Martin, in a thread "Abstraction idea" from 06/05/2008, to create a
class FunctionSpace where the mesh, element and dof_map(s) are
aggregated. Citing Martin:
U = FunctionSpace(mesh, dofmapset, form, 0) # or something similar
u = Function(U)
v = Function(U)

This seems a solid approach to me since it would provide a way of
encapsulating the mathematical formulation of the problem, which is
more or less const and likely to be reused by many discrete
functions in a solver.

It seems to me that there is an obvious risk that a lot of redundant
initialization would occur if all discrete functions should own
their own elements and dof maps. There seems to be consensus that
the mesh should be "global" for efficiency reasons, so why not treat
the function space the same way?

Is there a problem with an approach where the funciton _always_ owns
the vector and _never_ owns the function space (and mesh)? A very
strict design would avoid shared/smart pointers, provide a
comprehensible user interface and probably help the parallellization
effort.

/Dag
If the Function always owns the vector, there are cases you'll have
to make unneccessary copies of a vector, in particular such scenarios
may occur when trying to combine dolfin with something else.

If the Function never owns the function space, it must always be
constructed explicitly by the user. This may not be a bad thing.
However, if the Function is loaded from a file, nobody owns the
FunctionSpace.
Conceptually, I agree with Dag (and Martin?) that it is natural to
have global function spaces. And if the explicit construction of such
spaces can be made simple, it may not be a bad thing but a natural
part in setting up the mathematical problem. And I do not really like
that functions should be initialized from a form, which defines an
equation.

I think one idea was to not force less mathematically oriented users
to worry about function spaces. I guess there are (at least) 2 types
of functions: (i) functions part of the form, and (ii) functions not
part of the form, but used in pre/postprocessing etc.

For (i) it may be natural to construct the function space from the
form, and for (ii) it may be convenient in some cases, but it is not
really obvious that this is the best solution.

Maybe an explicit construction of a function space can come with a
default, such as a nodal basis of piecewise linears?

/Johan
So:
  FunctionSpace V(mesh);
  Function f(V);
gives a function f on piecewise linears?
That's ok with me.


About ownership, I think the only both robust and intuitive solution
is that an object never should store a reference (or regular pointer)
to another object. But as long as we are aware of the cost of doing
this and state it clearly in the documentation, I'm ok with keeping
e.g. the Mesh references like we do.

Any time object A stores a reference to object B, the user must
take care that the lifetime of B exceeds the lifetime of A. There
are no exceptions to this. This puts some real limitations on the
way the user must structure his program, e.g. he must sometimes
(often?) keep objects around longer than they're explicitly needed.

This may be a good thing, since it forces the user to think about
dependencies and object lifetimes, and the objects in question
use some memory.
I think this is ok. There are many ways to create a segfault in C++.
If you program in C++, you will have to think about memory.

But if we use references instead of shared_ptr,
we should never have default values:
- A Function has a reference to a Mesh, which is ok since
  it's always created outside.
- If a DiscreteFunction is to have a reference to a Vector, or a
  Function is to have a reference to a FunctionSpace, it cannot
  create its own without adding memory management code.

Every place we accept these limitations and requirements of how
the user structures his programs, we can use references and be
done with it. But don't think that the pretty syntax means the user
doesn't have to think about memory management, since all the
responsibility for memory management (object destruction order)
is in fact placed on the user, and errors from the users side will
lead to invalid references we cannot detect and segfaults.
I want the syntax to be simple and pretty, but I don't necessarily
want to hide the user from problems that are part of the design of
C++. It isn't Python or Java. You should be expected to know what
you are doing. :-)
It's not only about knowing what you're doing. It forces very
hard restrictions on the design/flow of your program, which can

1) Be a major source of bugs in nontrivial apps (see Garths email), which
are not locally visible because they depend on the global program flow.

2) Make it impossible to initialize e.g. Function in e.g a file reader,
since the caller of the file reader would need to get the objects
Function depends on. This is not limited to file readers, but is
a recurring pattern in nontrivial apps.

If we want to use dolfin or want dolfin to be used in apps that
are more complicated than the traditional "read input, compute
something, output something" app, these restrictions become
a larger problem.

Anyway, I like the idea about having a FunctionSpace class which
several Functions may share. The problem we need to solve is
reading from file:

 FunctionSpace V(mesh);
 Function u(V);
 file >> u;

The last line just fills out the data in both u and V.

This will lead to side effects as V might be changed when doing

 FunctionSpace V(mesh);
 Function u(V);
 Function v(V);
 file >> u;

V will be changed, both for u and v. In fact, the mesh will also be
changed.

The best thing would be if we could do

 file >> (mesh, V, u);
This is _exactly_ the kind of issue that smart pointers solve.

Btw, I tried to search the swig documentation for shared_ptr, and
found nothing...
I don't know what exactly they mean by "shared_ptr support".
It seems to be a set of typemaps that should kick in at the "right places".
They are defined in

   <Lib/python/boost_shared_ptr.i>

and used very rudimentary in

  <Examples/test_suite/li_boost_shared_ptr.i>

in the source tree of the 1.3.36 release. It seems that it is not specific for
boost::share_ptr but should also figure out tr1::shared_ptr

I think:

  %include <boost_shared_ptr.i>

at the appropriate place should do the trick.

Johan
I was about to start sketching on a FunctionSpace class when I
realized that this class might look different for different types of
Functions.

In particular, some functions (DiscreteFunctions) need to have a
dof map and a finite element, while this is not needed for other
functions.

This means that we might need to duplicate the hierarchy of different
function classes in a number of different function space classes.

I see two other options:

1. Require that a function space consists of

  mesh, finite_element, dof_map

and pick a suitable (trivial) dof map for example constant functions.

2. Implement a number of different function space classes but have a
single Function class which holds data but where each function call
is passed on to the specific type of function call in the
corresponding function space class.
Here are some more thoughts about the shared_ptr issue.

As I understand it, using shared_ptr to store the data in
DiscreteFunction means that users will need to use shared_ptr.
In particular, a Function must be initialized with a shared_ptr to a
Mesh, not a reference to a Mesh. Then in all demos, we would have to
change from

  Mesh mesh("mesh.xml");

to

  shared_ptr<Mesh> mesh = new Mesh("mesh.xml");

This doesn't look very pretty, and it would make the interface
inconsistent. One would need to know that some classes likes Mesh
should "always" be stored in a shared_ptr (since it will most
certainly be used to initialize a Function), while other classes like
Matrix, File, KrylovSolver, can be created as usual.

There seem to be two options:

1. Don't use shared_ptr anywhere (at least not in public interfaces)

If we choose this option, we will need to decide who owns what. For
example, we could say that a Function always owns the Vector, but
never the Mesh.

This works out also for reading from files:

  Mesh mesh;
  Function u(mesh);
  file >> u;

This would lead to side effects (the mesh will be changed), but that
might be ok.

2. Use shared_ptr everywhere

If we choose this option, we need to redefine Mesh so that it is not a
mesh class, but a shared_ptr for a mesh class (which we rename to
something else). Sundance does something like this and it gives a nice
and pretty syntax:

  Mesh mesh = new Mesh();             // no need to delete
  Function u = new Function(mesh);    // no need to delete

Essentially, it would look like Java. See here:

  http://software.sandia.gov/sundance/basics.html

in the section "Handles and Memory Management".

I think we need to pick one of these options and not a mix since it
would make the interface inconsistent.

Opinions? 1 or 2?

I don't think that mixing is inconsistent. At the moment we use a
mixture of references and pointers, and there is usually a logical
choice. Most users never need to work with the pointers. If we use smart
pointers, they will be used primarily internally in place of where we
currently use plain pointers.

For the specific case of a Function, part of the issue will be resolved
by requiring that a Function own it's own vector (forgetting about
sub-functions for the time being). Ownership of a finite element, or
possibly a FunctionSpace in the future, remains an issue because it's
natural that these are shared.

If you look at uBLASVector which now uses a smart pointer, not a single
piece of code outside uBLASVector had to be changed, but an advanced
user can now create a uBLASVector from a uBLAS object and not worry
about it going out of scope,

Garth
I'm all for it if it can be hidden (or visible only by choice), but
isn't it the case that the mesh (or function space) will always have
to be a shared_ptr?

In all applications, one will need to write

  shared_ptr<Mesh> mesh = new Mesh(...);

I think that it's logical for a Function to have a reference to a Mesh.
I'm not sure exactly the nicest way to deal with a FunctionSpace. If we
want pretty syntax + plus smart pointer flexibility, a Function could
have a pointer to the FunctionSpace, in which case it *does not* own the
FunctionSpace, and a smart pointer in which case it may or may not  own
the FunctionSpace.

The class Function would look like

  class Function
  {
    public:

      Function(FunctionSpace& space) : V(space) {}

      Function(shared_prt<FunctionSpace> space) : V(0), shared_V(space)
        {}

    private:
      FunctionSpace* V;
      shared_ptr<FunctionSpace>  shared_V;
  }

and a Function could be initialised by

  FunctionSpace V;
  Function u(V);

in which case it does not own the FunctionSpace, or by

  shared_ptr<FunctionSpace> V(new FunctionSpace);
  Function u(V);

in which case it may or may not own the FunctionSpace. The former would
be used in simple programs, and the latter in, for example,
LinearPDE::solve.

Garth

?
Sounds a bit dangerous to me... It does enable sharing objects
if you want to, but keeps the possibility to shoot of your foot.

Anyway, a little tip: shared_ptr can take a deleter policy as input
to its constructor, so you don't need the extra FunctionSpace* to
implement this approach. Something like this should do the trick:

      class NoDeleter {
        public:
          void operator() (FunctionSpace *p) {}
      };


I don't follow this.

...

      Function(FunctionSpace& space) : shared_V(&space, NoDeleter()) {}

(not tested, but modified from some working code).


And by the way, a thing to look out for is when using shared_ptr
for array pointers. Then you need to use a deleter policy like above,
except with function body { delete [] p; } instead.

This sounds like the trick I was asking for earlier: a way to provide
two different ways to initialize objects, either passing a reference
or a shared_ptr (something like decrease_count).

Do we want to do it this way?


It's not ideal in terms of safety, but it does provide what we want - simplicity and flexibility.

Garth



------------------------------------------------------------------------

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




Follow ups

References