Proposals for the DUNE meeting 2010 in Münster

AD: I tried to summarize some of the proposals from the discussion page, please feel free to modify the proposal if my summary seems wrong.

AD: Add below any of the points were you a bit more explanation is required (please being as specific as possible). But please do not start a discussion here!

CE: Please describe the problems. Often it is not clear why a certain change is proposed. Thus it is not possible to discuss the suggested solution and/or alternative approaches.

1: Performance issues

1.a: EntityKey

• Proposer: Robert and Martin
• Problem description: For “storage” of Entities in array-like structures we introduce the concept EntityKey (appropriate name to be found) which only contains the minimal information needed to create the corresponding Entity from this. An appropriate method will be available on the Grid class.
• Proposal: An example implementation is available in ALUGrid (3d). The EntityKey only has to store the minimal information needed to create an Entity(Pointer). This will allow to efficicently create subsets of grids which at the moment is not really possible in an good way. A genernal EntityPointer implementation could be based on EntityKeys and thus reducing implementation effort for grid implementers.
• Impact on Users: None, the new concept is not a must-use, but a can use.

1.b: Entity (pointer)

• Proposer: Robert and Martin

• Problem description: There are a number of problems with the EntityPointer.

  • It does not behave like a pointer and in fact it can not do so, because Entities can be on-the-fly objects.\
  • Dereferencing an EntityPointer returns a reference, therefore the EntityPointer basically has to contain the Entity so that two concepts are rolled into one.\
  • Storing Entities through the EntityPointer is not really feasible because EntityPointers contain the Entity. To have allow light weight storage the concept of an EntityKey is suggested. For many grids this need not contain more than a pointer or a few indices. The grid then provides methods to convert an EntityKey to an Entity.\
  • Having fewer classes to code reduces the work for implementing a new grid.\
  • There should be some performance gain for example with MetaGrids since they do not have to store there own EntityPointer the HostEntityPointer…

• Proposal:

  Removal of EntityPointer: The transition should be not to difficult since the typedef for EntityPointer will point to Entity and on the Entity class we introduce the method

  const Entity* operator -> () const { return this; }
  

  After this the interface will be much better to understand better to use. In fact, all methods returning a reference to an object should be revised.

  1. Storage of entities: For “storage” of Entities in array-like structures we use the concept EntityKey, see 1a).

1.c: Copy vs. Reference

• Proposer: Robert and Martin

• Problem description: As in the previous topic this is about a class having to store all kinds of objects just to be able to return references. This leads to rather large objects (Intersections stores a Geometry, two EntityPointer, which again each store an Entity and a Geometry…) For the method on the Geometries the object returned should not be restricted to a FieldMatrix/FieldVector anymore. Many grids (Yasp/Prism for example) could return a light weight object on which all methods could be implemented much faster. We need an interface or a concept for at least a constant matrix. Here the methods on DenseMatrix could be used.

• Proposal:

  1. Have methods in grid interface return objects instead of references. In general we should remove most references. This will eliminate problems related to lifetime of a cache object etc. When ever copies are expensive we should us use lightweight interface classes and move the storage to a shared container class.
     Examples are the geometry methods, inside/outside, index set methods.
  2. Change return value of jacobianInverseTransposed: not only should the methods on the geometry return objects instead of references also the mapping methods should noy return a FieldMatrix but some class satisfying a matrix like interface.

1.d: Semantics of method contains()

• Proposer: ???

• Problem description: ???

• Proposal:

  Fix behaviour of the method contains() if the entity that’s being tested is not from the same grid as the grid view to be undefined.
  AD: is this the proposal?

Interface cleanup

• Proposer: ???

• Problem description: ???

• Proposal:

  Remove the following methods from the GenericReferenceElement:
  - checkInside (const FieldVector< ctype, dim-codim > &local, int i) const
  - global
  - initializeTopology()
  - initialize()

  AD: There seems to be no other finished proposal for this topic yet?

2: Interface extensions

subIndex interface

• Proposer: Robert and Martin
• Problem description: The GeometryGrid can use vectors (based on (host) vertex indices) to store the position of the vertices. To create the geometry of an entity, it needs to acquire the (host) indices of the vertices for the entity. For 0 < codim < dim, these indices cannot currently be obtained, because the method subIndex( entity, i, entity.mydimension ) does not exist. Access to an element having this entity as a subentity is not available either.
• Proposal: Index Sets should be able to return subIndices for all codimension entities.
• Impact on Developers: Only grids implementing entities for 0 < codim < dim are really affected by this change (for elements the method already exists and for vertices its trivial). For dune-grid, this is only AlbertaGrid, ALUGrid, and SGrid (GeometryGrid is trivial; it simply forwards the method).
• Impact on Users: None, only additional method variants are introduced.

Persistent container

• Proposer: Robert and Martin
• Problem description: Currently when a data is stored that should be persistent also in an adaptive computation one would have to use LocalIdSet. These ids have to be attached to the grid in some way (e.g., ALUGrid and AlbertaGrid have a “hierarchic index” to attach data to the grid). Though technique is very grid depend, it might be the most efficient. Can we provide access to the technique to the user?
• Proposal: Add a utility class persistent Index container. AlbertaGrid or ALUGrid export their hierarchic indices through the very grid specific HierarchicIndexSet (the implementation is even suboptimal for AlbertaGrid). Together with these, DUNE-Fem implements a concept for fast adaptation of DoF vectors, without copying to a (hash)map and back. To support this concept by a wider variety of grids without such an undefined thing as a HierarchicIndexSet, we suggest to introduce a utility class that, by default, uses the current IdSet to implement a fast, persistent storage using a hash map. For some grids, like AlbertaGrid or ALUGrid, it can then be specialized to do whatever the grid developer thinks most efficient (true believers can still use the LocalIdSet). Even more, for ALUGrid and AlbertaGrid, we can implement the LocalIdSet based on this container (which is nothing but the current implementation, thus reducing code duplication).
• Impact on Developers: Implementation of the Default Persistent container (already done); specialization for grids like ALUGrid (done) and AlbertaGrid. Adding a default implementation for a LocalIdSet based on the persistent container (only if the persistent container is specialized).
• Impact on Users: None, the concept of a LocalIdSet shall remain untouched. Only a utility class is introduced.
• Discussion:
  CE: The problem description only covers half the truth:
  a) hash-based containers offer O (1) access, even with the LocalIdSet
  b) the persistent container dictates the user a trade-off between speed and memory
  RK: hash-based containers offer C O (1) access, where C >= 1. This is not the same as with some efficient implementation (e.g. in ALUGrid).
  MN: ALUGrid implements the IdSet through its “HierarchcIndexSet”. In that sense, you already have the waste of memory + the performance loss. Moreover, it is an external class and you don’t have to use it. It just reflects what the grid implementor thinks optimal. The idea of DUNE is to have slim interfaces, enforcing the LocalIdSet might be quite heavyweight. By the way: The case lies similar with Alberta: The PersistentContainer (though it is not called that) is already provided by the grid, DUNE simply uses it to store the LocalIdSet. Why should we enforce such waste of memory within our “slim” interface?
  CE: OK, then it is still not clear to me how your concepts work. I though you have this sparse index, called HierarchicIndex and you want to store data in a vector with holes. In this case, sotrying data in vectors would imply a waste of memory, not the sparse index itself. If you store data in a hash structure, you have the freedom to use larger hash and reduce the danger of collisions to nearly 0, or to use a small hash and waste CPU cycles for either collision detection, or the search in a bucket. Btw. your last argument is rather ridiculous, as you are suggesting to solve, what you describe as a design error, by repeating the design.
  AD: As with many other parts in dune the decision on how to implement this will be different for each grid implementation. The technique of using a hash map based on the local id can be used for any grid but for some grids alternative approaches are possible (see ALU) and we want to make them available to the user. If the waste of memory caused by ALU is not acceptable for a user then he/she can use an id set based version or switch the grid implementation.
  RK: The PersistentContainer does not waste memory since it allows explicitly to store data for “each” entity in the hierarchical grid. At the moment this can only be done with some “cuke” implementation.
  If I used this concept to store leaf data, then indeed, this would be a slight waste of memory. However, this is not the intent of this class.

I/O

• Proposer: Robert and Martin

• Problem description: For large-scale computation, being able to save a checkpoint and resume the computation from such a checkpoint is vital (you cannot afford to lose all data near the end of the computation). A key problem here is writing out and reading back the (adaptive) hierarchical grid.

• Proposal: Methods for Backup and Restore (reading and writing of grids) on the grid class: We propose the following interface for a factory-like structure (see also FS#80):

  template struct BackupRestoreFactory { void backup( const Grid& grid, const std::string& path, const std::string& fileprefix );

  Grid* restore( const std::string& path, const std::string& fileprefix );
  

  };

Facilitate implementation/testing of new features

• Proposer: Robert and Martin
• Problem description: Adding a new feature like the EntityKey to dune-grid is very obscure. Just adding a method like key() to the entity implementation only solves half the problem. How can the user call this method? How can he know that the grid implements an extension like EntityKeys (otherwise I might want to resort to EntityPointers)?
• Proposal: Add a public method impl() to each interface class in dune-grid (e.g., Entity, IndexSet, …). This allows easy access to non-standard features. Code using these features can easily be recognized from the use of impl(). Additionally, we introduce a namespace Extensions, which constains “capabilities” about such extensions (e.g., hasEntityKey).
• Impact on Developers: Add the method impl() to all interface classes. Add a file “extensions.hh” containing the extension capabilities (maybe in dune/grid/utility?)
• Impact on Users: No direct impact. Will allow easier access to non-standard features (no workaround to access getRealImplementation needed).
• Discussion:
  CE what exactly is the proposed solution?
  AD Is more description still required?

3: GridFactory interface

3.a: Problems with boundary projection and load balancing/checkpointing

• Proposer: Robert and Martin

• Problem description: During I/O operations (like load balancing or backup / restore), grids currently have to store and read back boundary segments. Due to the virtual interface, this needs very complicated (and currently undeveloped) mechanisms. Moreover, boundary segments are not part of the grid hierarchy, but form a kind of callback during adaptation. This contradicts our paradigm of strict separation of grid and user data (e.g., callbacks).

• Proposal: Make the boundary projections “user data”. A “boundary segment handle” could be given to adapt(), which is called on insertion of a new vertex. The callback interface could look as follows:

  template< class Grid > struct BoundarySegmentHandle { FieldVector< typename Grid::ctype, Grid::dimensionworld > project( const typename Grid::LeafIntersection &boundary, const FieldVector< typename Grid::ctype, Grid::dimension-1 > &x ) const; };

• Discussion:
  CE what exactly is the proposed solution?
  CE So basically you are pointing out a problem, but not yet proposing a solution. Right?! So this is a point for open discussion.
  MN: The proposal is to remove a complicated construction and replace it with something more DUNE-style. Implementing the callback strategy in ALUGrid would be quite easy, since the described interface is nearly what ALU uses internally. The proposal is – again – to remove code duplication by moving certain problems to the user. We can also provide a utility class, that provides the current behavior. The user will only have to pass the “BoundarySegmentHandle” of this utility class to the adapt method of the grid.

3.b: ParameterTree interface for the GridFactory

• Proposer: Jö Fahlke

• Problem description: Often the program will read a configuration file via the ConfigParser, making that configuration the natural place to let the program user specify something like heapSize.

• Proposal:

I propose the give the GridFactory an additional constructor signature of the form

  GridFactory(const ParameterTree &);

This constructor is meant for passing options from a (user written) parameter file to the grid creation process. An empty ParameterTree here is equivalent to calling the default constructor. Parameters unknown to the particular grid manager are ignored.

1.  Use 1. to support passing UG's `heapSize` parameter as a
    "heap\_size" `ParameterTree` option to the `GridFactory`.

• Open Questions: How does this interact with FS#699: “A general dynamic parameter interface for grids”?

3.c: ParameterTree interface for the StructuredGridFactory

• Proposer: Jö Fahlke

• Description: The StructuredGridFactory is useful to create structured coarse grids for Dune’s unstructured grid managers. Often the parameters for these grids should be specified via a configuration file. It would be useful to let the StructuredGridFactory do the interpretation of these options. This avoids having to rewrite the interpretation time and time again, and allows that the options have standardized defaults and names.

• Solution Proposal:

  1. On the StructuredGridFactory add the signatures

     static shared_ptr<GridType>
     createSimplexGrid(const FieldVector<ctype,dimworld>& lowerLeft,
     const FieldVector<ctype,dimworld>& upperRight,
     const array<unsigned int,dim>& elements,
     const ParameterTree &params);
     static shared_ptr<GridType>
     createCubeGrid(const FieldVector<ctype,dimworld>& lowerLeft,
     const FieldVector<ctype,dimworld>& upperRight,
     const array<unsigned int,dim>& elements,
     const ParameterTree &params);
     

     The additional ParameterTree parameter is simply passed to the underlying GridFactory and not otherwise acted upon by the StructuredGridFactory itself. This is meaningful only in connection with the proposal “ParameterTree interface for the GridFactory”.

1. On the StructuredGridFactory add the signatures

   static shared_ptr<GridType>
   createSimplexGrid(const ParameterTree &params);
   static shared_ptr<GridType>
   createCubeGrid(const ParameterTree &params);
   

   These extract the properties of the grid to create from the ParameterTree and then call the functions introduced in 1. The following options are recognized:

   • “bbox.lower” -> lowerLeft (default: vector of 0)
   • “bbox.upper” -> upperRight (default: vector of 1)
   • “elements” -> elements (default: array of 1)

2. On the StructuredGridFactory add the following ParameterTree option to the to the functions introduced in 2.

   • “global_refines” -> number of globalRefines() to do after the grid has been created according to the other options. (default: 0)

• Open Questions: Why can’t we use an abstract interface for the parameter tree?

3.d: Parallel grid creation

AD: There does not seem to be a proposal, just different ways things are done (UG vs. ALU).

4: dune-localfunctions

• Proposer: Jö Fahlke
• Description: dune-localfunctions currently provides an interface and a lot of implementations for local-valued finite elements. Local-valued means that the values returned are given in the coordinate system of the reference element, which makes a difference for vector-valued finite elements. For vector-valued finite elements the local values can usually be transformed into global ones, but the transformation is generally different for different kind of finite elements. Some vector-valued finite elements might even be easier to implement such that they return the global values directly.
  Another issue making sure that the finite elements are matched correctly between neighboring grid elements. This is currently only generically possible for scalar-valued elements with at most one DoF per sub-entity. In general doing this correctly requires both transformations of the global shape function values (usually multiplication by -1) as well as permutations of shape functions.
• Solution Proposal: Main differences to local-valued interface (Details in dune-localfunctions/doc/dune-localfunctions-manual.tex (PDF: /pdf/meetings/2010-11-devmeeting/dune-localfunctions-manual.pdf), sample implementation in dune-pdelab):
  • BasisTraits now needs to distinguish between local and global domain.
  • Basis (and thus the FiniteElement) stores the information necessary to do all the transformations ==> Finite element is constructed anew for each grid element.
  • FiniteElementFactory class with a semi-standard create() method. FiniteElementFactory is responsible for pre-computing information. The create() method cannot be standardized completely since the information required by different kinds of finite elements is simply too, well, different. There are however certain pieces of information which reoccur often: the geometry of the grid element, the ordering of its vertices and the orientation of faces. The encoding of this information should be standardized:
    • Standard interface for Geometry information (Dune::Geometry/BasicGeometry).
    • Standard interface for global ordering information of the vertices (new). One VertexOrdering object provides ordering information for sub-entities of a certain subset of the codimensions. There are different implementations of VertexOrdering objects planned, differing in the subset of codimensions they support.
      • A general VertexOrdering class using the information obtained from global ids exists. It provides information for all codimensions.
  • Wrapper classes for simple local-valued finite elements exist.

• Open Questions:
  • What about dynamic polymorphism?
  • What about adaptation?
  • The VertexOrdering stuff certainly needs more input.

• Discussion:
  CG: Creating the FE anew for each element might be expensive due to the construction of the local keys from the VertexOrdering (see e.g. Pk3d). One solution would be to cache FEs depending on the vertex ordering (at least if they do not depend on the geometry). Another solution would be to cache/precompute the LocalKeys/LocalFEs for all possible orderings in the factory and use them on creation. Both variants require access to some index/id/hash for the VertexOrdering that can be computed and compared/hashed as fast as possible.
  MN+RK: dune-localfunctions was intended for (local) shape functions and should not depend on any grid-specific stuff.
  • The basis should not contain more information than the (local) basis functions themselves. This allows using the same Basis in, e.g., a DG context.
  • If an additional transformation (e.g., Piola transformation), then we can add such a class and put a typedef into the local finite element.
  • Multiplying basis functions with minus 1 (how can I write the number in this Wiki?) can be done by negating the DoF. So at most the criterion is needed in dune-localfunctions (which should be a general description – like the LocalKeys – and not a fixed implementation).

4.a: dune-grid interface to obtain topologic information for global functions

• Proposer: Christian Engwer
• Problem description: dune-localfunctions require additional information to handle the mapping from the reference element to the grid correctly. As Jö pointed out it is not sufficient to specify the local-global-mapping, but certain topological and orientation details are required as well. Up to now there is no general interface for these information. There are e.g. the so-called twist used in dune-fem or the variant used in the Pk implementation of Sven. Neither of these concepts is sufficient for the general case and it might expensive to obtain the necessary information through the current grid interface. In order to write generally usable localfunctions all localfuncitons implementations should use the same interface to the topology of a local element. dune-grid should provide an additional interface to obtain the missing information needed for dune-localfunctions. This interface should be grid agnostic, so that we on’t introduce additional dependencies between dune-grid and dune-localfunctions.
• Proposed solution: Currently there are two approaches:
  a) On the localfunctions meeting in Heidelberg Andreas and Martin suggested to add something he called twist-provider to the grid interface. It is an enumeration of all possible “twists” (or permutations of corner Ids) that could possibly exist in a grid (i.e., just the “identity” for YaspGrid). Then, each such number (0, …, n) must be assigned semantics (e.g., the permutation). Instead of a permutation of corner Ids, we suggested a geometric interpretation because it is easier to understand for us. But the key idea is easy access to the twist (no ordering of vertices required). In this sense, this approach is a superset of b).
  b) Sven and myself suggested to use a permutation instead of a coordinate mapping. Discussing this frther, we realized, that the permutation isn’t enough, as it doesn’t provide enough information to detect a flip-operation in the mapping of the element. In general we have to provide the same information the geometry for geometric aspects for the topological aspects. Discussing this further and testing certain finite elements, we think the an mapping for the indices, and an integer for the orientation (the sign of the jacobian) should be enough.
  We should decide on such an interface on the meeting to make different implementations of local functions compatible.
  AD: I added the suggestion from the Heidelberg meeting (one pdf and the interface class twists.pdf, twistprovider.hh)

5: Flyspray tasks

• 841: Remove dependency to WML
• 830: Use svn tags for bugfix releases
• 827: rbegin of FieldVector and ISTL classes has surprising semantics
• 819: Remove enum FieldVector::size
• 818: Topological embedding and indices of subentities are missing
• 801: Int2Type vs. integral_constant
• 791: Unify interface class LevelIterator and LeafIterator
• 766: Separate GenericGeometries and the SmallObjectPool
• 746 / 661: Extract ‘function from expression’ functionality from dgf parser
• 743: Add a way to determine whether a grid is hybrid
• 699: A general dynamic parameter interface for grids
• 670: Remove obscure “DUNE developer mode”
• 532: Entity<0>::count should have a dynamic codim argument
• 515: Move method ‘contains’ from IndexSet to GridView
• GridFamilyTraits to restrictive

6: Organisation

6.a: branches

6.b: developer status

6.c: Make discretization modules core modules

AD: What is the suggestion here and please describe what for you a core module is. I always thought that core modules contained interfaces which all developers had agreed on and where changes to the interfaces were voted on. To me therefore neither dune-fem nor dune-pdelab seem to meet the requirements of a core module. We could try to agree on some more discretization interfaces which we can put into dune-localfunction (or into some new module), which would be a good thing.