]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add tests to check that the projection of an element yields the correct approximation.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 11 Aug 2006 17:11:55 +0000 (17:11 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 11 Aug 2006 17:11:55 +0000 (17:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@13672 0785d39b-7218-0410-832d-ea1e28bc413d

tests/deal.II/Makefile
tests/deal.II/project_common.cc [new file with mode: 0644]
tests/deal.II/project_q_01.cc [new file with mode: 0644]
tests/deal.II/project_q_01/cmp/generic [new file with mode: 0644]
tests/deal.II/project_rt_01.cc

index 3fdea2819ae97b96e612554a250b734015929700..8a469d27bf3f26bc4f3a0e09e69a7c0b6e1aac43 100644 (file)
@@ -54,7 +54,7 @@ tests_x = block_matrices \
        interpolate_* \
        have_same_coarse_mesh_* \
        get_finest_common_cells_* \
-       project_*
+       project_*[0-9]
 
 # from above list of regular expressions, generate the real set of
 # tests
diff --git a/tests/deal.II/project_common.cc b/tests/deal.II/project_common.cc
new file mode 100644 (file)
index 0000000..faebd59
--- /dev/null
@@ -0,0 +1,156 @@
+//----------------------------  project_common.cc  ---------------------------
+//    $Id: project_common.cc 12732 2006-03-28 23:15:45Z wolf $
+//    Version: $Name$ 
+//
+//    Copyright (C) 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  project_common.cc  ---------------------------
+
+
+// common framework to check whether an element of polynomial order p can
+// represent functions of order q
+
+#include "../tests.h"
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <numerics/vectors.h>
+#include <fe/fe_abf.h>
+#include <fe/fe_dgp.h>
+#include <fe/fe_dgp_monomial.h>
+#include <fe/fe_dgp_nonparametric.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_nedelec.h>
+#include <fe/fe_q.h>
+#include <fe/fe_q_hierarchical.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_system.h>
+
+#include <fstream>
+#include <vector>
+
+
+template <int dim>
+void test ();
+
+
+template <int dim>
+class F :  public Function<dim>
+{
+  public:
+    F (const unsigned int q,
+       const unsigned int n_components)
+                   :
+                   Function<dim>(n_components),
+                   q(q)
+      {}
+    
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component) const
+      {
+       Assert ((component == 0) && (this->n_components == 1),
+               ExcInternalError());
+       double val = 0;
+       for (unsigned int d=0; d<dim; ++d)
+         for (unsigned int i=0; i<=q; ++i)
+           val += (d+1)*(i+1)*std::pow (p[d], 1.*i);
+       return val;
+      }
+
+
+    virtual void vector_value (const Point<dim> &p,
+                              Vector<double>   &v) const
+      {
+       for (unsigned int c=0; c<v.size(); ++c)
+         {
+           v(c) = 0;
+           for (unsigned int d=0; d<dim; ++d)
+             for (unsigned int i=0; i<=q; ++i)
+               v(c) += (d+1)*(i+1)*std::pow (p[d], 1.*i)+c;
+         }
+      }
+    
+  private:
+    const unsigned int q;
+};
+
+
+
+template <int dim>
+void test_no_hanging_nodes (const FiniteElement<dim> &fe,
+                           const unsigned int        p)
+{
+  Triangulation<dim>     triangulation;
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (3);
+
+  DoFHandler<dim>        dof_handler(triangulation);
+  dof_handler.distribute_dofs (fe);
+
+  deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl;
+
+                                  // there are no constraints here
+  ConstraintMatrix constraints;
+  constraints.close ();
+
+  Vector<double> interpolant (dof_handler.n_dofs());
+  Vector<float>  error (triangulation.n_active_cells());
+  for (unsigned int q=0; q<=p+2; ++q)
+    {
+                                      // project the function
+      VectorTools::project (dof_handler,
+                           constraints,
+                           QGauss<dim>(p+2),
+                           F<dim> (q, fe.n_components()),
+                           interpolant);
+      
+                                      // then compute the interpolation error
+      VectorTools::integrate_difference (dof_handler,
+                                        interpolant,
+                                        F<dim> (q, fe.n_components()),
+                                        error,
+                                        QGauss<dim>(std::max(p,q)+1),
+                                        VectorTools::L2_norm);
+      deallog << fe.get_name() << ", P_" << q
+             << ", rel. error=" << error.l2_norm() / interpolant.l2_norm()
+             << std::endl;
+         
+      if (q<=p)
+       Assert (error.l2_norm() < 1e-12*interpolant.l2_norm(),
+               ExcInternalError());
+    }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile(logname);
+  logfile.precision (3);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<1>();
+  test<2>();
+  test<3>();
+}
+
diff --git a/tests/deal.II/project_q_01.cc b/tests/deal.II/project_q_01.cc
new file mode 100644 (file)
index 0000000..80042c6
--- /dev/null
@@ -0,0 +1,29 @@
+//----------------------------  project_q_01.cc  ---------------------------
+//    $Id: project_q_01.cc 12732 2006-03-28 23:15:45Z wolf $
+//    Version: $Name$ 
+//
+//    Copyright (C) 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  project_q_01.cc  ---------------------------
+
+
+// check that VectorTools::project works for Q elements correctly on
+// a uniformly refined mesh for functions of degree q
+
+char logname[] = "project_q_01/output";
+
+
+#include "project_common.cc"
+
+
+template <int dim>
+void test ()
+{
+  for (unsigned int p=1; p<6-dim; ++p)
+    test_no_hanging_nodes (FE_Q<dim>(p), p);
+}
diff --git a/tests/deal.II/project_q_01/cmp/generic b/tests/deal.II/project_q_01/cmp/generic
new file mode 100644 (file)
index 0000000..35cf8fe
--- /dev/null
@@ -0,0 +1,47 @@
+
+DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], P_0, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], P_1, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], P_2, rel. error=0.000278
+DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], P_3, rel. error=0.000654
+DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_0, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_1, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_2, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_3, rel. error=5.15e-06
+DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_4, rel. error=1.52e-05
+DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_0, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_1, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_2, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_3, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_4, rel. error=1.20e-07
+DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_5, rel. error=4.27e-07
+DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_0, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_1, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_2, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_3, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_4, rel. error=0
+DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_5, rel. error=3.21e-09
+DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_6, rel. error=1.33e-08
+DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], P_0, rel. error=0
+DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], P_1, rel. error=0
+DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], P_2, rel. error=8.74e-05
+DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], P_3, rel. error=0.000201
+DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_0, rel. error=0
+DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_1, rel. error=0
+DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_2, rel. error=0
+DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_3, rel. error=1.02e-06
+DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_4, rel. error=3.01e-06
+DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_0, rel. error=0
+DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_1, rel. error=0
+DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_2, rel. error=0
+DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_3, rel. error=0
+DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_4, rel. error=2.22e-08
+DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_5, rel. error=7.82e-08
+DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], P_0, rel. error=0
+DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], P_1, rel. error=0
+DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], P_2, rel. error=2.52e-05
+DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], P_3, rel. error=5.74e-05
+DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_0, rel. error=0
+DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_1, rel. error=0
+DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_2, rel. error=0
+DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_3, rel. error=2.04e-07
+DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_4, rel. error=6.07e-07
index b2b2040226892862fb24546d58efd89f54b7fec1..abd3beb019b11b123f47c9c29fd622f896a4d067 100644 (file)
 //----------------------------  project_rt_01.cc  ---------------------------
 
 
-// check that VectorTools::project works for RT elements correctly on
+// check that VectorTools::project works for Q elements correctly on
 // a uniformly refined mesh for functions of degree q
 
-#include "../tests.h"
-#include <base/function.h>
-#include <base/logstream.h>
-#include <base/quadrature_lib.h>
-#include <lac/vector.h>
+char logname[] = "project_rt_01/output";
 
-#include <grid/tria.h>
-#include <dofs/dof_handler.h>
-#include <dofs/dof_constraints.h>
-#include <grid/grid_generator.h>
-#include <grid/grid_refinement.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary_lib.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_tools.h>
-#include <numerics/vectors.h>
-#include <fe/fe_raviart_thomas.h>
-#include <fe/fe_system.h>
-
-#include <fstream>
-#include <vector>
-
-
-template <int dim>
-class F :  public Function<dim>
-{
-  public:
-    F (const unsigned int q,
-       const unsigned int n_components)
-                   :
-                   Function<dim>(n_components),
-                   q(q)
-      {}
-    
-    virtual void vector_value (const Point<dim> &p,
-                              Vector<double>   &v) const
-      {
-       for (unsigned int c=0; c<v.size(); ++c)
-         {
-           v(c) = 0;
-           for (unsigned int d=0; d<dim; ++d)
-             for (unsigned int i=0; i<=q; ++i)
-               v(c) += (d+1)*(i+1)*std::pow (p[d], 1.*i)+c;
-         }
-      }
-
-  private:
-    const unsigned int q;
-};
 
+#include "project_common.cc"
 
 
 template <int dim>
 void test ()
 {
-  Triangulation<dim>     triangulation;
-  GridGenerator::hyper_cube (triangulation);
-  triangulation.refine_global (3);
-
-  for (unsigned int p=0; p<6-dim; ++p)
-    {
-      FE_RaviartThomas<dim> fe(p);
-      DoFHandler<dim>        dof_handler(triangulation);
-      dof_handler.distribute_dofs (fe);
-
-      deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl;
-      
-      ConstraintMatrix constraints;
-      constraints.close ();
-
-      Vector<double> interpolant (dof_handler.n_dofs());
-      Vector<float>  error (triangulation.n_active_cells());
-      for (unsigned int q=0; q<=p+2; ++q)
-       {
-                                          // project the function
-         VectorTools::project (dof_handler,
-                               constraints,
-                               QGauss<dim>(p+2),
-                               F<dim> (q, fe.n_components()),
-                               interpolant);
-      
-                                          // then compute the interpolation error
-         VectorTools::integrate_difference (dof_handler,
-                                            interpolant,
-                                            F<dim> (q, fe.n_components()),
-                                            error,
-                                            QGauss<dim>(std::max(p,q)+1),
-                                            VectorTools::L2_norm);
-         deallog << fe.get_name() << ", P_" << q
-                 << ", rel. error=" << error.l2_norm() / interpolant.l2_norm()
-                 << std::endl;
-         
-         if (q<=p)
-           Assert (error.l2_norm() < 1e-12*interpolant.l2_norm(),
-                   ExcInternalError());
-       }
-    }
-}
-
-
-
-int main ()
-{
-  std::ofstream logfile("project_rt_01/output");
-  logfile.precision (3);
-  
-  deallog.attach(logfile);
-  deallog.depth_console(0);
-  deallog.threshold_double(1.e-10);
-
-  test<2>();
-  test<3>();
+  if (dim != 1)
+    for (unsigned int p=0; p<6-dim; ++p)
+      test_no_hanging_nodes (FE_RaviartThomas<dim>(p), p);
 }
-

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.