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.