surface_quad.cc

#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/manifold_lib.h>

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>


#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/mapping_q_eulerian.h>
#include <deal.II/fe/mapping_fe_field.h>
#include <deal.II/fe/fe_system.h>

#include <deal.II/numerics/vector_tools.h>

#include <cmath>
#include <iostream>
#include <fstream>
#include <string>

using namespace dealii;
int main(int argc, char ** argv)
{
  const unsigned int quad_order = 35;
  const unsigned int quad_max = 50;
  const double exact_area=4.*numbers::PI;
  const unsigned int ref_max=2;
  const unsigned int degree_max=9;
  const double tol=1e-10;
  std::cout<<"Testing the surface computation of a sphere"<<std::endl;
  std::cout<<"exact area : "<<exact_area<<std::endl;
  for(unsigned int i=1; i<=degree_max; ++i)
  {
    std::cout<<i<<" ";
    for(unsigned int j=5; j<=quad_max; j=j+5)
    {
      double surface_area=0.;

      Triangulation<2,3> tria;
      Point<3> center;

      GridGenerator::hyper_sphere	(	tria,center,1. );
      SphericalManifold<2, 3> manifold(center);
      // Manifold<2,3>
      // manifold = std::shared_ptr<SphericalManifold<2, 3> > (new SphericalManifold<2, 3> (pippo));
      tria.set_all_manifold_ids(0);
      tria.set_manifold(0, manifold);
      tria.refine_global(3);

      FE_Q<2,3> fe(i);
      FESystem<2, 3> fee(FE_Q<2, 3> (i), 3);

      DoFHandler<2,3> dh(tria);
      DoFHandler<2,3> dhh(tria);
      dh.distribute_dofs(fe);
      dhh.distribute_dofs(fee);
      std::cout<<dh.n_dofs()<<" ";
      Vector<double> euler_vec(dhh.n_dofs());
      VectorTools::get_position_vector(dhh,euler_vec);
      // MappingFEField<2,3> mapping(dhh,euler_vec);
      MappingQ<2,3> mapping(i);


      QGauss<2> quad(j);
      const unsigned int n_q_points    = quad.size();

      FEValues<2,3> fe_values (mapping, fe, quad,
                                   update_values |
                                   update_quadrature_points |
                                   update_JxW_values);

     std::vector<Point<3> > support_points(dh.n_dofs());
     DoFTools::map_dofs_to_support_points<2, 3>(mapping,
             dh, support_points);

      for(unsigned int i=0;i<support_points.size();++i)
        if(std::fabs(support_points[i].square()-1.)>tol)
        {
          std::cout<<"ERROR!!! : "<<std::fabs(support_points[i].square()-1.)<<std::endl;
          break;
        }
    //  for(unsigned int i=0;i<support_points.size();++i)
    //  if(std::fabs(support_points[i].square()-1.)<tol)
    //     {
    //       std::cout<<"OK!!! : "<<std::fabs(support_points[i].square()-1.)<<std::endl;
    //     }

      for (typename DoFHandler<2,3>::active_cell_iterator
        cell = dh.begin_active(),
        endc = dh.end();
        cell!=endc; ++cell)
      {
        fe_values.reinit(cell);
          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
        {
           surface_area += fe_values.JxW(q_point);
        }

      }
      // std::cout<<surface_area<<" "<<exact_area<<std::endl;
      double error = std::fabs(surface_area-exact_area)/exact_area;
      // std::cout<<"degree "<<i<<" error "<<error<<std::endl;
      std::cout<<error<<" ";
      std::string output_filename;
      output_filename="sphere_order_"+Utilities::int_to_string(i)+"_refinement_"+Utilities::int_to_string(j)+".vtu";
      std::ofstream output_file(output_filename);
      GridOut().write_vtu (tria, output_file);
      tria.set_manifold(0);
    }
    std::cout<<std::endl;
  }
}

@luca-heltai When I download and run this code on my computer, I get this

Testing the surface computation of a sphere
exact area : 12.5664
1 386 0.00833771 386 0.00833771 386 0.00833771 386 0.00833771 386 0.00833771 386 0.00833771 386 0.00833771 386 0.00833771 386 0.00833771 386 0.00833771 
2 1538 4.86182e-06 1538 4.86182e-06 1538 4.86182e-06 1538 4.86182e-06 1538 4.86182e-06 1538 4.86182e-06 1538 4.86182e-06 1538 4.86182e-06 1538 4.86182e-06 1538 4.86182e-06 
3 3458 3.34359e-09 3458 3.346e-09 3458 3.34599e-09 3458 3.346e-09 3458 3.34599e-09 3458 3.34598e-09 3458 3.34598e-09 3458 3.34603e-09 3458 3.34604e-09 3458 3.34597e-09 
4 6146 6.78804e-10 6146 9.84731e-10 6146 9.84725e-10 6146 9.84734e-10 6146 9.84694e-10 6146 9.84735e-10 6146 9.84723e-10 6146 9.847e-10 6146 9.84723e-10 6146 9.84764e-10 
5 9602 2.40824e-09 9602 2.94218e-09 9602 2.94219e-09 9602 2.94218e-09 9602 2.94222e-09 9602 2.94215e-09 9602 2.94217e-09 9602 2.94218e-09 9602 2.94218e-09 9602 2.94214e-09 
6 13826 1.06568e-07 13826 3.41674e-09 13826 3.41674e-09 13826 3.41675e-09 13826 3.41679e-09 13826 3.41676e-09 13826 3.41678e-09 13826 3.41676e-09 13826 3.41679e-09 13826 3.41675e-09 
7 18818 1.76091e-07 18818 4.99957e-09 18818 4.9868e-09 18818 4.9868e-09 18818 4.98679e-09 18818 4.98682e-09 18818 4.9868e-09 18818 4.98681e-09 18818 4.9868e-09 18818 4.98684e-09 
8 24578 5.91816e-07 24578 4.56815e-09 24578 4.56721e-09 24578 4.56721e-09 24578 4.56727e-09 24578 4.56722e-09 24578 4.56725e-09 24578 4.56724e-09 24578 4.5672e-09 24578 4.56716e-09 
9 31106 8.31487e-06 31106 5.09134e-09 31106 5.3417e-09 31106 5.34175e-09 31106 5.34174e-09 31106 5.34173e-09 31106 5.3417e-09 31106 5.34168e-09 31106 5.3417e-09 31106 5.34174e-09 

These errors are much larger than what Luca posted here dealii/dealii#2992 (comment)

I am using this deal.II

commit f9196e343617dd5aafea4abff43704f35f7ac7b0