[Dune] YaspGrid, periodicity and boundary segment indices

[Atgeirr Rasmussen]

Dear Dune,

I have some problems understanding what the proper behaviour of boundary segment indices should be.
I am currently updating our grid class CpGrid to properly support periodicity, and so I am looking at YaspGrid to get an impression of what's correct.
(I am currently on revision 7633 of dune-grid, in case that matters)

The program listed at the bottom constructs a 3x3x3 YaspGrid with user-controlled overlap,
that is periodic in all directions. Then I iterate over the interior partition, and print segment index and ids for all boundaries.
When I run the program (called periodic_grid_test in my system) with overlap 0, like this:

./periodic_grid_test 0 | sort -n

... it is easy to see that I get 54 unique indices, [0, ..., 53]. However, with overlap 1, like this:

./periodic_grid_test 1 | sort -n

... I get only 9 different indices, each repeated 6 times, once for each boundary id.
The indices are also not consecutive, but numbered 6, 7, 8, 11, 12, 13, 16, 17, 18.
The same applies with larger overlap (just different actual indices).

I assumed that these runs would produce identical results, but the do not.
Have I totally misunderstood the boundary segment indices? Or is there a bug in YaspGrid?

Any information will be greatly appreciated!

Atgeirr Flø Rasmussen
SINTEF ICT


----------------- Test program below ---------------------

#include <dune/grid/yaspgrid.hh>
#include <cstdlib>

template <class GridView>
void investigateBoundary(const GridView& g)
{
    const Dune::PartitionIteratorType ptype = Dune::Interior_Partition;
    typedef typename GridView::template Codim<0>::template Partition<ptype>::Iterator CI;
    typedef typename GridView::IntersectionIterator FI;
    for (CI c = g.template begin<0, ptype>(); c != g.template end<0, ptype>(); ++c) {
        for (FI f = g.ibegin(*c); f != g.iend(*c); ++f) {
            if (f->boundary()) {
                std::cout << f->boundarySegmentIndex() << ' ' << f->boundaryId() << '\n';
            }
        }
    }
}


int main(int argc, char** argv)
{
    int overlap = (argc > 1 ? std::atoi(argv[1]) : 0);
    Dune::FieldVector<double, 3> dim(1.0);
    Dune::FieldVector<int, 3> num(3);
    Dune::FieldVector<bool, 3> periodic(true);
    Dune::YaspGrid<3> grid(dim, num, periodic, overlap);
    investigateBoundary(grid.leafView());
}