Volume method for multilinear geometries

For non-affine multilinear geometry mappings the volume method doesn't calculate the correct result. Instead of using numerical integration with higher order it uses midpoint rule as approximation:

/** \brief obtain the volume of the mapping's image
 *
 *  \note The current implementation just returns
 *  \code
 *  integrationElement( refElement().position( 0, 0 ) ) * refElement().volume()
 *  \endcode
 *  which is wrong for n-linear surface maps and other nonlinear maps.
 */
ctype volume () const
{
  return integrationElement( refElement().position( 0, 0 ) ) * refElement().volume();
}

Compare also to:

https://dune-project.org/doxygen/pdelab/master/group__GeometryReferenceElements.html

The question is how to fix it:

• Estimate the correct quadrature order and use an appropriate quadrature rule.
• Introduce a parameter that can be set by the user.
• Do not fix it at all since the error in the discretization is probably not too large.

I would be interested to hear some opinions.