Alberta subIndex returns the same index for different entities
- The code below shows the issue. In the first loop I store the vertex coordinates of an entity using the
index
method. Later, I loop over subEntity
sub-entities and ask for its index with the subIndex
method. The results between index
and subIndex
are inconsistent and the assertion at the end fails.
using Index = typename Grid::LeafIndexSet::IndexType;
using Coordinate = typename Grid::template Codim<0>::Geometry::GlobalCoordinate;
auto gv = g->leafGridView();
for(const auto& entity : elements(gv)) {
std::map<Index, Coordinate> corners;
// store all the corners of the grid using the index set
for (const auto& subEntity : subEntities(entity, Dune::Codim<Grid::dimension>{}))
corners[gv.indexSet().index(subEntity)] = subEntity.geometry().center();
Dune::Hybrid::forEach(std::make_index_sequence<Grid::dimension+1>{}, [&](auto cd) {
for (const auto& subEntity : subEntities(entity, Dune::Codim<cd>{})) {
// loop over all vertex sub-entities of subEntity
for (std::size_t s = 0; s != subEntity.subEntities(Grid::dimension); ++s) {
// compare the s-th vertex coordinate with the one stored above
auto index = gv.indexSet().subIndex(subEntity, s, Grid::dimension);
Coordinate e_corner = corners[index];
Coordinate se_corner = subEntity.geometry().corner(s);
if (e_corner != se_corner)
DUNE_THROW(Dune::RangeError, e_corner << " != " << se_corner);
}
}
});
}
- The image below shows a P1-2D finite element interpolation using the
subIndex
method to index its degrees of freedom on a 3D grid. The discontinuity is given by the problem mentioned above, where the wrong coefficient is chosen to interpolate the element.
- Indeed, the whole issue is marked as a TODO in the commit that introduced the method: dae993c7