Testing sphere surface quadrature

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=14;
   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<<"deg=" << i<<" ";
      for(unsigned int j=i+1; j<=i+1; ++j)
      {
         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_DGQArbitraryNodes<2,3> fe(QGaussLobatto<1>(i+1));
         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<<"ndofs=" << 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);
         
         
         QGaussLobatto<2> quad(j);
         QGaussLobatto<1> face_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);
         FEFaceValues<2,3> fe_face_values (mapping, fe, face_quad,
                                           update_values |
                                           update_quadrature_points |
                                           update_normal_vectors |
                                           update_JxW_values);
         FEFaceValues<2,3> fe_face_values_nbr (mapping, fe, face_quad,
                                               update_values |
                                               update_quadrature_points |
                                               update_normal_vectors |
                                               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;
         //     }
         
         double max_dw= 0, max_dn = 0;
         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);
            }
            for(unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
            {
               fe_face_values.reinit(cell, f);
               fe_face_values_nbr.reinit (cell->neighbor(f), cell->neighbor_of_neighbor(f));
               for (unsigned int q_point=0; q_point<face_quad.size(); ++q_point)
               {
                  double dw = std::fabs(fe_face_values.JxW(q_point) - fe_face_values_nbr.JxW(q_point));
                  dw /= fe_face_values.JxW(q_point);
                  max_dw = std::max(max_dw, dw);
                  Tensor<1,3> dn = fe_face_values.normal_vector(q_point) +
                  fe_face_values_nbr.normal_vector(q_point);
                  max_dn = std::max(max_dn, dn.norm());
               }
            }
            
         }
         // 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<<", area error =" << error<<std::endl;
         std::cout << "max dweight =" << max_dw << "  max dnormal=" << max_dn << std::endl;
         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;
   }
}

@nicola-giuliani @luca-heltai I modified this gist

https://gist.github.com/nicola-giuliani/86f535ddfd83485150767d123c95c242

by replacing FE_Q with FE_DGQArbitraryNodes(QGaussLobatto<1>(degree+1)) and QGauss with QGaussLobatto. For degree=i, I use (i+1)-GaussLobatto quadrature. This is less accurate than QGauss so I run upto degree=14. The result is this

Testing the surface computation of a sphere
exact area : 12.5664
1 1536 0.00823725 
2 3456 1.21078e-05 
3 6144 3.48582e-08 
4 9600 8.77328e-09 
5 13824 9.8126e-09 
6 18816 9.66812e-09 
7 24576 1.27884e-08 
8 31104 1.12567e-08 
9 38400 1.28967e-08 
10 46464 1.24116e-08 
11 55296 1.36827e-08 
12 64896 1.05048e-08 
13 75264 1.4113e-08 
14 86400 1.37372e-08

Convergence is stalled around 1e-8.

@nicola-giuliani @luca-heltai In the above code, if I use FE_DGQArbitraryNodes(QGauss<1>(degree+1)) and QGauss(degree+1) for quadrature, I get

Testing the surface computation of a sphere
exact area : 12.5664
1 1536 ERROR!!! : 0.00813246
0.00833776 
2 3456 ERROR!!! : 9.19648e-07
4.86606e-06 
3 6144 ERROR!!! : 2.56938e-06
1.22582e-09 
4 9600 ERROR!!! : 6.87922e-07
6.78803e-10 
5 13824 ERROR!!! : 2.48012e-08
2.69808e-09 
6 18816 ERROR!!! : 4.63677e-07
3.1846e-09 
7 24576 ERROR!!! : 3.90752e-07
4.70775e-09 
8 31104 ERROR!!! : 7.55285e-08
4.34411e-09 
9 38400 ERROR!!! : 1.28744e-08
5.09134e-09 
10 46464 ERROR!!! : 1.30113e-07
5.1399e-09 
11 55296 ERROR!!! : 1.30596e-07
5.78457e-09 
12 64896 ERROR!!! : 2.34806e-08
4.55343e-09 
13 75264 ERROR!!! : 6.54167e-09
6.14423e-09 
14 86400 ERROR!!! : 5.53619e-08
6.08515e-09

Since this issue seems to disappear when changing w into 1-w, as you indicate in one of the issues, my guess is that there is an internal problem in the computation of the "intermediate support points" of the finite elements.

I'm looking into this...