MappingManifold does not support update_hessians and similar higher order mappings

[tamiko]

Related to #2478:

Currently, the following code fails:

constexpr int fe_order = 1;
dealii::FE_Q<dim, spacedim> finite_element(fe_order);
dealii::MappingManifold<dim, spacedim> mapping;
dealii::QGauss<dim> quadrature(fe_order);

dealii::FEValues<dim, spacedim> fe_values(
    mapping, finite_element, quadrature, dealii::update_hessians);

because the required transformations are not implemented: source/fe/mapping_manifold.cc

   253   // Now throw an exception if we stumble upon something that was not            
   254   // implemented yet                                                             
   255   Assert(!(out & (update_jacobian_grads | update_jacobian_pushed_forward_grads | 
   256                   update_jacobian_pushed_forward_2nd_derivatives |               
   257                   update_jacobian_pushed_forward_3rd_derivatives)),              
   258          ExcNotImplemented());   

@luca-heltai ping

[bangerth]

I believe that we talk about that in the ACM TOMS paper with @luca-heltai and me as authors. Implementing this requires computing second derivatives of the function that describes the manifold. (Equivalently, we need derivatives of the tangent vectors.) The Manifold class does not provide an interface for this.

[luca-heltai]

I believe that we talk about that in the ACM TOMS paper with @luca-heltai and me as authors. Implementing this requires computing second derivatives of the function that describes the manifold. (Equivalently, we need derivatives of the tangent vectors.) The Manifold class does not provide an interface for this.

This is correct. In this case a reasonable approach to provide a good approximation would be to build a MappingQ(p) with support points inserted according to the manifold, and ask higher derivatives to the MappingQ instead of to MappingManifold. This is what happens currently if you build a MappingQ. If you want to mimic what a MappingManifold would do, then you should use a MappingQ with support points coinciding with your quadrature points. This guarantees you have quadrature points exactly on the geometry, and gives you a good polynomial approximation of the Hessian.