]> https://gitweb.dealii.org/ - dealii.git/commitdiff
debugging RaviartThomasNodal
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 25 May 2005 09:14:10 +0000 (09:14 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 25 May 2005 09:14:10 +0000 (09:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@10740 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_poly_tensor.cc
tests/fe/rtn_1.cc

index d18bb445ed84c757ec30142f46b3a25b64d98755..fa50303cfcd34ecb35aed2203f12c6529dc4753a 100644 (file)
@@ -228,10 +228,10 @@ FE_PolyTensor<POLY,dim>::get_data (const UpdateFlags      update_flags,
          if (inverse_node_matrix.n_cols() == 0)
            for (unsigned int i=0; i<this->dofs_per_cell; ++i)
              data->shape_values[i][k] = values[i];
-         else
-           for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-             for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-               data->shape_values[i][k] = inverse_node_matrix(j,i) * values[j];
+         else
+           for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+             for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+               data->shape_values[i][k] += inverse_node_matrix(j,i) * values[j];
        
        if (flags & update_gradients)
          if (inverse_node_matrix.n_cols() == 0)
@@ -240,7 +240,7 @@ FE_PolyTensor<POLY,dim>::get_data (const UpdateFlags      update_flags,
          else
            for (unsigned int i=0; i<this->dofs_per_cell; ++i)
              for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-               data->shape_grads[i][k] = inverse_node_matrix(j,i) * grads[j];
+               data->shape_grads[i][k] += inverse_node_matrix(j,i) * grads[j];
       }
   return data;
 }
@@ -282,7 +282,7 @@ FE_PolyTensor<POLY,dim>::fill_fe_values (
   
   for (unsigned int i=0; i<n_dofs; ++i)
     {
-      const unsigned int first = i*dim;//data.shape_function_to_row_table[i];
+      const unsigned int first = data.shape_function_to_row_table[i];
       
       if (flags & update_values)
        for (unsigned int k=0; k<n_quad; ++k)
index 170ef7218ec26d7d73dc2671757cc91ccd882bdd..75f7ce463a8c213e093c6dec075f1efb275057f7 100644 (file)
@@ -48,7 +48,7 @@ check_support_points (const FiniteElement<dim>& fe)
   MappingCartesian<dim> mapping;
   
   Quadrature<dim> support_quadrature(points, weights);
-  UpdateFlags flags = update_values;
+  UpdateFlags flags = update_values | update_gradients;
   FEValues<dim> vals(mapping, fe, support_quadrature, flags);
   vals.reinit(dof.begin_active());
   
@@ -62,13 +62,15 @@ check_support_points (const FiniteElement<dim>& fe)
          for (unsigned int d=0;d<dim;++d)
            {
              const double sf = fe.shape_value_component(i,p,d);
-             const double sv = vals.shape_value_component(i,k,d);
-             deallog << ' ' << (int) round(60*sf)
-                     << '/' << (int) round(60*sv);
+             deallog << ' ' << (int) round(60*sf);
              
-//           const double diff = std::abs(sf - vals.shape_value_component(i,k,d));
-//           if (diff > 1.e-12)
-//             deallog << "Error" << vals.shape_value_component(i,k,d) << std::endl;
+             const double diff = std::abs(sf - vals.shape_value_component(i,k,d));
+             if (diff > 1.e-12)
+               deallog << "Error values" << std::endl;
+             Tensor<1,dim> grad = fe.shape_grad_component(i,p,d);
+             grad -= vals.shape_grad_component(i,k,d);
+             if (grad.norm() > 1.e-12)
+               deallog << "Error grads" << std::endl;
            }
        }
       deallog << std::endl;
@@ -77,6 +79,60 @@ check_support_points (const FiniteElement<dim>& fe)
 }
 
 
+template <int dim>
+void
+check_face_support_points (const FiniteElement<dim>& fe)
+{
+  deallog << "Face " << fe.get_name() << std::endl;
+  
+  const std::vector< Point<dim-1> > & sub_points = fe.get_unit_face_support_points();
+  std::vector<double> weights(sub_points.size());
+  Quadrature<dim-1> sub_quadrature(sub_points, weights);
+  std::vector< Point<dim> > points(sub_points.size());
+  
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+  
+  MappingCartesian<dim> mapping;
+  
+  UpdateFlags flags = update_values | update_gradients;
+  FEFaceValues<dim> vals(mapping, fe, sub_quadrature, flags);
+
+  for (unsigned int face=0;face<GeometryInfo<dim>::faces_per_cell;++face)
+    {
+      QProjector<dim>::project_to_face(sub_quadrature, face, points);
+      vals.reinit(dof.begin_active(), face);
+      
+      for (unsigned int k=0;k<points.size();++k)
+       {
+         const Point<dim>& p = points[k];
+         deallog << p;
+         for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+           {
+             deallog << '\t';
+             for (unsigned int d=0;d<dim;++d)
+               {
+                 const double sf = fe.shape_value_component(i,p,d);
+                 deallog << ' ' << (int) round(60*sf);
+                 
+                 const double diff = std::abs(sf - vals.shape_value_component(i,k,d));
+                 if (diff > 1.e-12)
+                   deallog << "Error values" << std::endl;
+                 Tensor<1,dim> grad = fe.shape_grad_component(i,p,d);
+                 grad -= vals.shape_grad_component(i,k,d);
+                 if (grad.norm() > 1.e-12)
+                   deallog << "Error grads" << std::endl;
+               }
+           }
+         deallog << std::endl;
+       }
+      deallog << std::endl;
+    }
+}
+
+
 int
 main()
 {
@@ -87,6 +143,8 @@ main()
 
   FE_RaviartThomasNodal<2> e20(0);
   check_support_points(e20);
+  check_face_support_points(e20);
   FE_RaviartThomasNodal<2> e21(1);
   check_support_points(e21);
+  check_face_support_points(e21);
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.