Call for Change

Subject: Indices, attaching external data to grid entities


0) INTRO
========

We discussed at several occassions the syntax and semantics of the
index methods in the entities and were still not settled. The
requirements are:

- Mapping of entities to entries in a contiguous array
- Fast, i.e. O(1) (with SMALL constant), access from entity to data
- Fast access from codim 0 entity to all subentities
- Levelwise and leafwise organization of data arrays
- Globally unique id's of entities for parallel algorithms
- Efficient reorganization of data arrays in case of adaptive refinement
- Possibility of specialized implementations for the simpler cases

The solution that I devised requires the following changes:

- Define four index methods in every entity
- Define additional four index methods in codim 0 entity (subindex)
- New size method in grid for size of compressed leaf index set
- Suggestion that entities also deliver old compressed index after
  mesh adaptation
- Suggestion for the interface of a mapper class
- Finally the mappers should be placed in dune/grid/common since they
  can be used to attach any type of data with a grid and not just
  degrees of freedom


1) IN THE ENTITIES
==================

Four different indices are required from each entity:

- global index

This index unique for every entity in the whole grid. It is also
persistent over adaptation of the grid.

Method:

  Grid::PersistentIndexType persistentIndex () const;

- compressed level index

This index is unique and consecutive for all entities of a given codim
on a given level in a given process.

Method:

  int compressedIndex () const;

- global leaf index

This index is globally unique and persistent for all leaf entities of
the grid. The index is inherited by all copies of an
entity (as defined in the "Grid White Paper"). It is accessibly at
least on all copies that are part of a codim 0 leaf entity. This is
formally defined in the "Grid theory" document.

Method:

  Grid::PersistentIndexType persistentLeafIndex () const;

- compressed leaf index 

This index is unique and consecutive for all leaf entities of a given codim
in a given process. The index is inherited by all copies of an
entity (as defined in the "Grid White Paper"). It is accessibly at
least on all copies that are part of a codim 0 leaf entity.


Method:

  int compressedLeafIndex () const;


- sub index

The indices of subentities of a codim 0 entity are accessible directly
from the codim 0 entity. This gives another four methods only on codim
0 entities:

Suggested methods:

  PersistentIndexType subPersistentIndex () const;
  int subCompressedIndex () const;

  PersistentIndexType subPersistentLeafIndex () const;
  int subCompressedLeafIndex () const;

Note: These methods are member templates parametrized by the codimension.

- New type Grid::PersistentIndexType

This type is exported by each Grid to denote the data type used for
the persistent indices. In parallel implementations on 32-bit
processors these cannot be assumed to fit into an int. I wrote a class
template 

		 bigunsignedint<k> in dunu/common/bigunsignedint.hh

that support arbitrily large unsigned ints where the number of bits to
use is a template parameter. The PersistentIndexType must have the
operators ==, !=, <, >, <=, >= and an output function should be
supplied for debugging.


RATIONALE:
----------

The compressed indices are used by the subsequently described mapper
classes to map an entity to an entry in an array (of some user defined
type). This is needed for the leaf entities (e.g. assembly of the global
system to be solved) as well as for the entities on a given level
(e.g. in the geometric multigrid method).

The persistent indices are used only during grid adaptation. The
mapper could use them to build up a map of "old compressed indices". Then,
after adaptation, the new compressed index can be related to the old
compressed index for each entity.


2) IN THE GRID
==============

- Size methods

The new method

  int size (int codim) const;

should return the number of leaf entities of the given codim on that
process, i.e. the size of the compressed leaf index set.

The old method

  int size (int level, int codim) const;

still returns the number of entities of the given codim on the given
level, i.e. the size of the compressed level index.


Note: These methods should be implemented as efficiently as
possible. A mapper (discussed below) might use them frequently.

The following methods are useful when allocating stiffness matrices
  (calculation of nonzeroes), so I would suggest to add them:

//! number of entities per level, codim and geometry type in this process
int size (int level, int codim, GeometryType type) const;


//! number of leaf entities per codim and geometry type in this process
int size (int codim, GeometryType type) const;

Note: These methods would also allows us to determine at run-time
whether a grid is all simplex or all cubes.


3) MAPPER
=========

Data is stored in consecutive arrays external to the grid. In order to
access the data of an entity the entity must be mapped to an entry in
an array. More generally a fixed number of data items n_c might be
associated with an entity of codim c. Note that n_c here is fixed per codim
but might be different for each codim. 

A mapper class provides such a mapping. Its interface might look like:


  /*! This implementation manages a fixed number of entries for a single codimension.
   * Template parameters are:
   *    G:     a Dune grid type
   *    c:     a valid codim
   *    nc:    number of entries to associate with entities of codim c
   */
  template <typename G, int c, int nc>
  class SimpleMapper {
  public:

	//! associate mapper with a grid
	SimpleMapper(G& grid);

	//! get total number of entries associated with leaf entities
	int size () const;

	//! get total number of entries associated with a given level
	int size (int level) const;
 
	//! map j'th item of entity to entry
	template<int cc> // this is necessary for multiple codim mappers
	int leafmap (const typename G::Traits::template codim<cc>::Entity& e, int j) const; 
  
	//! map j'th item of entity to entry
	template<int cc> // this is necessary for multiple codim mappers
	int levelmap (const typename G::Traits::template codim<cc>::Entity& e, int j) const;  
  
	//! map j'th item of subentity of codim 0 entity to entry
	template<int cc> // this is now the subentity's codim
	int subleafmap (const typename G::Traits::template codim<0>::Entity& e, int i, int j) const; 
  
	//! map j'th item of subentity of codim 0 entity to entry
	template<int cc> // this is now the subentity's codim
	int sublevelmap (const typename G::Traits::template codim<0>::Entity& e, int i, int j) const; 
  

	// methods for adaptation phase
	// WHY don't we require from the mesh to provide old compressed indices?

	//! start adaptation phase, backup old map
	void preadapt();

	//! terminate adaptation phase
	void postadapt();

	//! map j'th item of entity to entry
	template<int cc> // this is necessary for multiple codim mappers
	int oldleafmap (const typename G::Traits::template codim<cc>::Entity& e, int j) const; 
  
	//! map j'th item of entity to entry
	template<int cc> // this is necessary for multiple codim mappers
	int oldlevelmap (const typename G::Traits::template codim<cc>::Entity& e, int j) const;  
  
	//! map j'th item of subentity of codim 0 entity to entry
	template<int cc> // this is now the subentity's codim
	int oldsubleafmap (const typename G::Traits::template codim<0>::Entity& e, int i, int j) const; 
  
	//! map j'th item of subentity of codim 0 entity to entry
	template<int cc> // this is now the subentity's codim
	int oldsublevelmap (const typename G::Traits::template codim<0>::Entity& e, int i, int j) const; 

  private:
	G& g;
  };
 
This mapper here is a special version that allows n_c data items for
entities of codim c and n_i=0 for i!=c. 
It provides mapping to levelwise arrays (levelmap, sublevelmap) and to
arrays containing data associated with leaf entities only (leafmap,
subleafmap). These mapping methods can be implemented efficiently using the
compressed indices of the entities.

In the case of mesh adaptation the mapping of entities to data array
entries changes (in general) and the data arrays have to be
reorganized (data has to be copied from old arrays to new arrays, a
solution with partially filled arrays is not useful when solving
linear systems). This requires that after adaptation the old and the
new mapping for each entity are available. There are two possibilities
to accomplish this:

- The mapper uses the persistent indices to build up a map (in the
  sense of the STL) for the old mapping. The index methods at the
  entity always deliver the new index. This requires at least O(NlogN)
  time and O(N) space where N is the number of entities.

- The entities provide also the old mapping through additional member
  functions. This can be implemented in O(N) time and space (in UG,
  how about Alberta and ALU3d?) and is much easier

Therefore I suggest to move this functionality to the grid.

These methods would look like:

  int oldCompressedIndex () const;
  int oldCompressedIndex () const;
  template<int cc> int oldSubCompressedIndex (int i) const;
  template<int cc> int oldSubCompressedLeafIndex (int i) const;