Creating an Interface Adaption to Your Own Grid Data Structure

Creating an Interface Adaption to Your Own Grid Data Structure

If you want to use GrAL, you typically want to use its generic components, but not necessarily its concrete grid data structures - at least that's how it is intended. Therefore, you will need to create an interface layer satisfying the syntax specified in the Concepts section. Actual working code along the lines of this section can be found in the triang2d module.

For sake of simplicity, let us assume that your grid data structure is a plain good old C (or Fortran) one: Cell-vertex incidences are stored in a plain array cells, such that cells[3*c +v] gives the index of the v'th vertex of cell c. Likewise, an array geom holds the coordinates of vertex v at indices 2*v and 2*v+1.

Let's call our grid wrapper class Triang2D:

class Triang2D {
  int * cells;

  // to be continued ...
};

We are going to implement as much of the GrAL interface for Triang2D as possible. From the incidence information present in the array cells, we can define the following types more or less directly: Vertex, VertexIterator, Cell, CellIterator, VertexOnCellIterator. We will implement them as nested classes of Triang2D (note: the actual code used separate classes). Also, in order to save some work, we will implement the Vertex and VertexIterator by one and the same class Triang2D::Vertex, and do the same with Triang2D::Cell and Cell / CellIterator.

Moreover, as show below, we can also implement edges and related iterators. So, in sum, we can support the following concepts, without adding any new data structure (only the EdgeIterator will have an internal table to keep track of visited edges):

Element handles could be simply integers in this case; however, we chose to make them distinguishable types, which makes it easier to catch misuse errors, e.g. taking vertex handles for cell handles mistakingly.

class Triang2D {

  // typedef int vertex_handle; // simplest possibility
  typedef vertex_handle_int<Triang2D> vertex_handle;
  typedef cell_handle_int  <Triang2D> cell_handle;

  class Vertex {
     // to be continued ...
  };
  typedef Vertex VertexIterator;

  class Cell {
     // to be continued ...
  };
  typedef Cell   CellIterator;

  class VertexOnCellIterator {
     // to be continued ...
  };
};

We first show the definition of Triang2D::Cell, and then Triang2D::VertexOnCellIterator.

class Triang2D {
   class Cell {
     typedef Cell     self;
     typedef Triang2D grid_type;
   private:
      cell_handle      c;
      grid_type const* g;
   public:
      Cell() : g(0) {}
      Cell(grid_type const& gg, cell_handle cc = 0) : g(&gg), c(cc) {}  
   
      self      & operator++() { ++c; return *this;}
      self const& operator* () const { return *this;}
      bool IsDone() const { return c >= TheGrid.NumOfCells();}
  
      cell_handle handle() const { return c;}     
      grid_type   const& TheGrid() const { return *g; } 

      unsigned NumOfVertices() const { return 3; }
      VertexOnCellIterator FirstVertex() const 
       { return VertexOnCellIterator(*this);}
   };

};

(In practice, we would have to define FirstVertex() outside of Triang2D::Cell.) The class Triang2D::Vertex is of course very similar and not shown here.

Now for the incidence iterator Triang2D::VertexOnCellIterator which models VertexOnCellIterator:

class Triang2D {
 
  class VertexOnCellIterator {
    // typedefs omitted ...
  private:
    cell_handle c;
    int         lv;
    grid_type   g;
  public:
    VertexOnCellIterator(Cell const& cc) : c(CC.c), lv(0), g(CC.g) {}
     
    self&  operator++() { ++lv; return *this; }
    Vertex operator*() const { return Vertex(*g, handle());}
    bool   IsDone() const { return (lv >= 3);}
    
    vertex_handle       handle()  const { return g->cells[3*c+lv];}
    Cell                TheCell() const { return Cell(*g,c);}
    grid_type const&    TheGrid() const { return *g;}
  };

}

Now we have defined enough types to satisfy the Cell-VertexInputGridRange model. It suffices to serve as argument for a generic copy operation, see the ConstructGrid operation.

It remains to fill in the grid_types<> template:

template<>
struct grid_types<Triang2D> {
  typedef Triang2D grid_type;

  typedef Triang2D::vertex_handle vertex_handle;
  
  //  etc.
};

We can, however, support a larger portion of the GrAL interface. Even if we do not store edges, they are implicitly present as sides of triangles, and we can define an Edge type for Triang2D, as well as EdgeIterator and EdgeOnCellIterator. The types Facet etc. will be just typedefs to the corresponding edge concepts.

First, we note that we have two possibilities of representing an edge: First, by its two vertices, and second, by an incident cell and the local edge number in that cell, much like in Triang2D::VertexOnCellIterator.

Now, the second representation gives readily access to the two vertices of the edge, given our data structure, but not the other way around. So, we choose to represent edges (and edge-on-cell-iterators) by (cell, local edge number) pairs. If we have to compare two edges for equality, we must compare the vertex sets, as our primary representation is not unique (each interior edge is visited from both sides).

As the internal representations for Triang2D::Edge and Triang2D::EdgeOnCellIterator are identical, we will implement the first in terms of the latter. Starting with Triang2D::EdgeOnCellIterator, we see that it is very similar to Triang2D::VertexOnCellIterator:

class Triang2D {
 
  class EdgeOnCellIterator {
    typedef  EdgeOnCellIterator self;
    // other typedefs omitted ...
  private:
    cell_handle c;
    int         le;
    grid_type   g;
  public:
    EdgeOnCellIterator(Cell const& cc) : c(CC.c), le(0), g(CC.g) {}
     
    self&  operator++() { ++le; return *this; }
    Edge   operator*() const { return Edge(*this);}
    bool   IsDone() const { return (le >= 3);}

    Vertex V1() const { return Vertex(*g,g->cells[3*c+le]);}
    Vertex V2() const { return Vertex(*g,g->cells[3*c+(le+1)%3]);}
    

    edge_handle         handle()  const { return edge_handle(c,le);}
    Cell                TheCell() const { return Cell(*g,c);}
    grid_type const&    TheGrid() const { return *g;}

    friend bool operator==(self a, self b)
    { return (a.c == b.c && a.le == b.le);}
  };

}

The type Triang2D::edge_handle is a simple struct containing cell handle and local number.

The class Triang2D::Edge just contains an object of type Triang2D::EdgeOnCellIterator. The only major difference is the comparison operator ==:

class Triang2D {
  
  class Edge {
    EdgeOnCellIterator e;
  public:
    Edge(EdgeOnCellIterator ee) : e(ee) {}
    // ... most things forwarded to e

    friend bool operator==(Edge a, Edge b)
     { return (  (a.V1() == b.V1() && a.V2() == b.V2())
               || (a.V1() == b.V2() && a.V2() == b.V1())); }
  };
};

So far, so good. Now, how can we iterate over all edges of Triang2D grids? We can do so be iterating over all cells and iteration over all edges incident to those cells. As every internal edge is visited twice, we must keep track of already visited ones. This can be done with an hash table, taking edges as keys (for the definition of a hash function on edges, see below). So, Triang2D::EdgeIterator contains Triang2D::EdgeOnCellIterator and this hash table. As this strategy is often employed, it has been implemented in a generic component:

class Triang2D {

   class EdgeIterator 
    : public facet_set_of_cells_iterator<Triang2D::Cell>
   {
     typedef facet_set_of_cells_iterator<Triang2D::Cell> base;
    public: 
     EdgeIterator() : base(CellIterator()) {}
     EdgeIterator(CellIterator  const& c) : base(c) {}
     EdgeIterator(grid_type const& g)  : base(CellIterator(g)) {}
  };
};

BTW, the class facet_iterator_of_cell_set can also be used for finding the facets of (the closure of) an arbitrary subset of the cells of any grid.

To be written. See ArchetypedGrid.

Now, our combinatorial layer is complete, the iterators that are missing cannot be supported efficiently with the given data. If we want to use algorithms which require more (e.g. cell-cell adjacencies), we have to calculate the additional information. Probably, the most comfortable thing to do would be to use a view on Triang2D which supports the required functionality, reusing what is already there. At present, there is no view in GrAL for adding cell-cell adjacencies, but it would certainly be a useful addition.

If we want to use geometric algorithms, we need to define a geometry for Triang2D. At first sight, this looks rather simple, but there is a subtle problem, relating to the fact that vertex coordinates are stored in a plain array of doubles:

class geom_Triang2D_base 
   double* coords;
   typedef ... coord_type;  // 2D point type

   ??? const& coord(Vertex v) const;
   ???      & coord(Vertex v);
;

What should we return from the non-const coord(Vertex v)? On the one hand, we want allow code like coord_type p = geom.coord(v), on the other hand, an assignment like geom.coord(v) = coord_type(1,0) should change the array coords. Thus, we cannot return coord_type &.

The solution is to define a proxy type:

class geom_Triang2D_base 

   class coord_proxy 
     double * coo;
     void operator=(coord_type const& p) 
       coo[0] = p[0]; coo[1] = p[1]; 
   ;

   coord_proxy        coord(Vertex v) 
    return coord_proxy(coords[2*v.handle()]);
;

For the `real' coord_type, we are free to choose any suitable implementation of a 2D geometric point. In order to get a bunch of useful geometric functionality (like volumes or centers), we can use a generic 2D geometry:

 typedef geometry2d<geom_Triang2D_base>  geom_Triang2D;

The efficiency of this generic implementation will depend on how well the compiler can optimize out the overhead. If it turns out to be too inefficient, parts have to be specialized for Triang2D using direct access to its data. This does not break the generic paradigm, as grid geometries are part of the kernel.

This part is actually quite easy, as predefined components for the common cases exist. In the case of vertices and cells, which are numbered consecutively, we take the class template grid_function_vector as our implementation:

template<class T>
class grid_function<Triang2D::Vertex, T>
 : public grid_function_vector<Triang2D::Vertex, T> 

 // repeat constructors here.
;

The same is done for Triang2D::Vertex. The actual code can be found in grid-functions.h.

For edges, there is no consecutive numbering we can exploit. Thus, we have to use a different approach for grid functions on edges. One viable way is to use hash tables. This is so common (it is also used for partial grid functions), that we can again use a generic component:

template<class T>
class grid_function<Triang2D::Edge, T> 
  : public grid_function_hash<Triang2D::Edge, T> 

  // repeat constructors here.
;

In order to make it work this way, we have to somehow pass a hash function to the grid_function_hash<> template. This is done by specializing the traits class element_traits<> for Triang2D::Edge and providing a hasher_type:

template<>
struct element_traits<Triang2D::Edge> 
  : public element_traits_edge_base<Triang2D> 

  struct hasher_type : public grid_types_Triang2D 
    enum  p= 17; // somewhat arbitrary
    size_t operator()(Triang2D::Edge const& e) const  
      vertex_handle v1 = e.V1().handle();
      vertex_handle v2 = e.V2().handle();
      return (v1 > v2 ? p*v1 + (v2%p) : p*v2 + (v1%p));
    
  ;
;

The actual code is found in element-traits.h. We do the same for Triang2D::Vertex and Triang2D::Cell. Here the hash function just returns the handle.

For defining partial grid function, there effectively does not remain anything to do (given the provision of a hasher_type by element_traits above), because they are already defined in a generic way. Thus, in partial-grid-functions.h, we just include the generic partial-grid-function-hash.h.

If we just need a wrapper for read-only access to an existing data structure, we are fine with what we have achieved so far. With the interface created above, we are able to use all generic algorithms on Triang2D which do not change the grid. Provided, of course, the incidence information supplied by Triang2D is sufficient for the algorithm. If not, we would have to create additional data structures - the interface we have defined gives access to all information inherent in the data structure.

Now, we might for example want to read from a file format for which there is a GrAL input adapter:

Triang2D t;
InputWeirdFileFmt In("mygrid.weird");
ConstructGrid0(t, In);

This will copy the connectivity information from the file mygrid.weird to the triangulation t. In order for this to work, we need to specialize the ConstructGrid0 primitive for Triang2D. (Suffixes like 0 at the end of ConstructGrid0 are just technical hacks for distinguishing different overloaded versions.)

template<class GRID>
void ConstructGrid0(Triang2D     & t,
                    GRID    const& g_src);

In the general case, we need some relationship between g_src and t (associative copy), so we use that version as a master implementation which all other versions of the ConstructGrid family use:

template<class GRID, class VTX_CORR, class CELL_CORR>
void ConstructGrid0(Triang2D     & t,
                    GRID    const& g_src,
                    VTX_CORR     & src2t_v,
                    CELL_COOR    & src2t_c);

Afterwards, src2t_v will map the vertices of g_src to those of t in a 1:1 fashion, and src2t_c will do the same for cells.

The actual implementation is rather simple in our case:

1. The vertices of g_src are mapped to those of t, i.e., numbered consecutively from 0 to nv-1, where nv = g_src.NumOfVertices(). This creates the mapping src2t_v.
2. Next, storage is allocated for the cells array, using g_src.NumOfCells().
3. Finally, the cells of t are created, by using the cell and the vertex-on-cell iterator types of GRID, and the vertex correspondence of src2t_v. At the same time, the mapping src2t_c is filled in.

The complete code can be found in construct.C and copy.C.

To be written.

Most parts of the GrAL interface (iterators, grid functions) can be tested using the testing module.

For examples, have a look at test-triang2d-construct.C, test-triang2d-geometry.C, and the Triang2D-specific test-triang2d.C / test-triang2d.h.

Creating an Interface Adaption to Your Own Grid Data Structure