Add ReferenceElement::subEntities() method returning a range

For a reference element re the call

  re.subEntity(i,codim,j,c);

returns the index (wrt re) of the j-th sub-subentity of codim (wrt re) c of the subentity (i,codim). This allows e.g. to access all vertices of a face. However, there's a big problem with this method:

It is often misinterpreted by users that the passed index j is the index of one grid-subentity within another grid-subentity. To be more precise it is assumed that the subEntity() methods of the reference element and a grid element are consistent with each other. Unfortunately this is in general not true because it is impossible to construct such an consistent subentity embedding for any grid. As a consequence both may provide a different embedding and the index j is essentially useless despite allowing to iterate over all such subentities using

// Don't use j directly in such a loop, it's very likely that
// your assumptions on the meaning of j don't hold true!
for(std::size_t j=0; j<re.size(i,codim,c); ++j)
  foo(re.subEntity(i,codim,j,c));

While the behaviour is documented, users repeatedly fall into this trap and the issue is discussed on the list once in a while.

As a remedy this introduces a new method ReferenceElement::subEntities() returning all the subentities as an iterable range. The actual type of this range is implementation defined. Using this the above loop turns into:

// No j here!
for(auto k : re.subEntities(i,codim,c))
  foo(k);

This is not only more readable, but it especially avoids misinterpreting the 'index' j. You can still run into the trap on purpose by assuming a special ordering of the range. But this is very unlikely and you have to take special action to construct the index j.

This MR is based on !93 (merged) (contains) because it consistently tests for all codim/c combinations. But one could also restrict the allowed cases to c>=codim here. In a later stage we may want to deprecate the error prone subEntity() method.

Side note: Thanks to the implementer (@martin.nolte?) of the reference elements for caching those numbers in a consecutive way. Implementing this feature just amounts into wrapping pointers into the cache in a Dune::IteratorRange. However, I'd even more like to put them into gsl:span.

changed the description

added 3 commits

• edcac37f - Allow ReferenceElement::size(i,c,cc) for cc<c
• 4a45da2c - Add `ReferenceElement::subEntities()` method returning a range
• a387f1a6 - [test] Add test for new method `ReferenceElement::subEntities()`

Compare with previous version

I force pushed an updated version based on !93 (merged).

Any comments? I don't want to merge an interface change in a central part without prior discussion, even if it's just an interface addition.

Ping.

Another often needed feature is the reverse lookup: Does subentity (i,iCodim) of entity of the reference element contain subentity (j, jCodim) or, in terms of this proposal, is j contained in the range obtained by re.subEntities(i,iCodim,jCodim).

This could de implemented as re.contains(i,iCodim,j,jCodim). While this is currently not stored in the reference elements, it could simply be obtained and cached on construction.

A maybe more intuitive alternatively would be re.subEntities(i,iCodim,jCodim).contains(j). I could prepare an implementation once we at least decided on this merge request.

assigned to @smuething

added 2 commits

• fca01008 - Add contains() and size() method to subentity range
• 86542ec1 - Implement exact upper bound for number of subentities

Compare with previous version

I just added re.subEntities(i,iCodim,jCodim).contains(j) and re.subEntities(i,iCodim,jCodim).size(). Notice that the values of contains() are cached in a std::array<std::bitset>.

I compared this to std::array<std::vector<bool>> and std::vector<bool> with per-codim offset and the same with bool replaced by char. It turns out that the chosen solution performs slightly better and is also more readable than solutions with an offset instead of a nested container.

Notice that the new methods do not only increase readability of code and reduce sources of errors. The contains() method also provides a small performance improvement because it avoids a search loop that you had to do before.

@joe likes this and no one else seems to care. Feel free to merge, @smuething, after having a look at my changes to the interface code.

I also like this a lot. Will merge.

merged

mentioned in commit d721f3b8

I know I am late but can we please have a description of this in CHANGELOG.md?

Done: !116 (merged).