]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Tests to test point_value in the hp setting.
authorgoll <goll@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 26 Jan 2012 09:46:44 +0000 (09:46 +0000)
committergoll <goll@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 26 Jan 2012 09:46:44 +0000 (09:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@24930 0785d39b-7218-0410-832d-ea1e28bc413d

tests/bits/point_value_hp_01.cc [new file with mode: 0644]
tests/bits/point_value_hp_01/cmp/generic [new file with mode: 0644]
tests/bits/point_value_hp_02.cc [new file with mode: 0644]
tests/bits/point_value_hp_02/cmp/generic [new file with mode: 0644]

diff --git a/tests/bits/point_value_hp_01.cc b/tests/bits/point_value_hp_01.cc
new file mode 100644 (file)
index 0000000..aebfe10
--- /dev/null
@@ -0,0 +1,183 @@
+//----------------------------  point_value_hp_01.cc  ---------------------------
+//    Copyright (C) 2004, 2006, 2007 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.
+//
+//----------------------------  point_value_hp_01.cc  ---------------------------
+
+// check VectorTools::point_value for a hp dofhandler.
+
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.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/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vectors.h>
+
+#include <fstream>
+#include <cmath>
+#include <iomanip>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+  public:
+    MySquareFunction () : Function<dim> () {}
+
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const
+      {        return (component+1)*p.square()+1; }
+    
+    virtual void   vector_value (const Point<dim>   &p,
+                                Vector<double>     &values) const
+      { values(0) = value(p,0); }
+};
+
+
+template<int dim>
+class MyExpFunction : public Function<dim>
+{
+  public:
+    MyExpFunction () : Function<dim> () {}
+
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const
+      {        return std::exp (p(0)); }
+    
+    virtual void   vector_value (const Point<dim>   &p,
+                                Vector<double>     &values) const
+      { values(0) = value(p,0); }
+};
+
+
+
+template <int dim>
+void make_mesh (Triangulation<dim> &tria)
+{
+  
+  GridGenerator::hyper_cube(tria, -1, 1);
+
+                                   // refine the mesh in a random way so as to
+                                   // generate as many cells with
+                                   // hanging nodes as possible
+  tria.refine_global (4-dim);
+  const double steps[4] = { /*d=0*/ 0, 7, 3, 3 };
+  for (unsigned int i=0; i<steps[dim]; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator
+        cell = tria.begin_active();
+      for (unsigned int index=0; cell != tria.end(); ++cell, ++index)
+        if (index % (3*dim) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement ();
+    }
+}
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tria;
+  make_mesh (tria);
+  
+  hp::FECollection<dim> fe;
+  fe.push_back (FE_Q<dim>(3));
+  fe.push_back (FE_Q<dim>(4));
+  fe.push_back (FE_Q<dim>(5));
+
+  hp::DoFHandler<dim> dof_handler (tria);
+
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell->set_active_fe_index (std::rand() % fe.size());
+    }
+
+  dof_handler.distribute_dofs (fe);
+
+                                   // test with two different functions: one
+                                   // that is exactly representable on the
+                                   // chosen finite element space, and one
+                                   // that isn't
+  for (unsigned int i=0; i<2; ++i)
+    {
+      static const MySquareFunction<dim>          function_1;
+      static const Functions::CosineFunction<dim> function_2;
+
+      const Function<dim> &
+        function = (i==0 ?
+                    static_cast<const Function<dim>&>(function_1) :
+                    static_cast<const Function<dim>&>(function_2));
+      
+      Vector<double> v (dof_handler.n_dofs());
+      VectorTools::interpolate (dof_handler, function, v);
+
+                                       // for the following points, check the
+                                       // function value, output it, and
+                                       // verify that the value retrieved from
+                                       // the interpolated function is close
+                                       // enough to that of the real function
+                                       //
+                                       // also verify that the actual value is
+                                       // roughly correct
+      Point<dim> p[3];
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          p[0][d] = 0;
+          p[1][d] = 0.5;
+          p[2][d] = 1./3.;
+        }
+      Vector<double> value(1);
+      for (unsigned int i=0; i<3; ++i)
+        {
+          VectorTools::point_value (dof_handler, v, p[i], value);
+          deallog << -value(0) << std::endl;
+
+          Assert (std::abs(value(0) - function.value(p[i])) < 1e-4,
+                  ExcInternalError());
+
+         const double scalar_value = VectorTools::point_value (dof_handler, v, p[i]);
+          Assert (std::abs(value(0) - scalar_value) < 1e-4,
+                  ExcInternalError());
+        }  
+    }
+  
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("point_value_hp_01/output");
+  deallog << std::setprecision (4);
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/bits/point_value_hp_01/cmp/generic b/tests/bits/point_value_hp_01/cmp/generic
new file mode 100644 (file)
index 0000000..763c633
--- /dev/null
@@ -0,0 +1,22 @@
+
+DEAL:1d::-1.000
+DEAL:1d::-1.250
+DEAL:1d::-1.111
+DEAL:1d::-1.000
+DEAL:1d::-0.7071
+DEAL:1d::-0.8660
+DEAL:1d::OK
+DEAL:2d::-1.000
+DEAL:2d::-1.500
+DEAL:2d::-1.222
+DEAL:2d::-1.000
+DEAL:2d::-0.5000
+DEAL:2d::-0.7500
+DEAL:2d::OK
+DEAL:3d::-1.000
+DEAL:3d::-1.750
+DEAL:3d::-1.333
+DEAL:3d::-1.000
+DEAL:3d::-0.3536
+DEAL:3d::-0.6495
+DEAL:3d::OK
diff --git a/tests/bits/point_value_hp_02.cc b/tests/bits/point_value_hp_02.cc
new file mode 100644 (file)
index 0000000..c14af9a
--- /dev/null
@@ -0,0 +1,201 @@
+//----------------------------  point_value_hp_02.cc  ---------------------------
+//    Copyright (C) 2004, 2006, 2007 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.
+//
+//----------------------------  point_value_hp_02.cc  ---------------------------
+
+// check VectorTools::point_value for hp, alternative algorithm with mapping
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.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/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/numerics/vectors.h>
+
+#include <fstream>
+#include <cmath>
+#include <iomanip>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+  public:
+    MySquareFunction () : Function<dim> () {}
+
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const
+      {        return (component+1)*p.square()+1; }
+    
+    virtual void   vector_value (const Point<dim>   &p,
+                                Vector<double>     &values) const
+      { values(0) = value(p,0); }
+};
+
+
+template<int dim>
+class MyExpFunction : public Function<dim>
+{
+  public:
+    MyExpFunction () : Function<dim> () {}
+
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const
+      {        return std::exp (p(0)); }
+    
+    virtual void   vector_value (const Point<dim>   &p,
+                                Vector<double>     &values) const
+      { values(0) = value(p,0); }
+};
+
+
+
+template <int dim>
+void make_mesh (Triangulation<dim> &tria)
+{
+  
+  GridGenerator::hyper_cube(tria, -1, 1);
+
+                                   // refine the mesh in a random way so as to
+                                   // generate as many cells with
+                                   // hanging nodes as possible
+  tria.refine_global (4-dim);
+  const double steps[4] = { /*d=0*/ 0, 7, 3, 3 };
+  for (unsigned int i=0; i<steps[dim]; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator
+        cell = tria.begin_active();
+      for (unsigned int index=0; cell != tria.end(); ++cell, ++index)
+        if (index % (3*dim) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement ();
+    }
+}
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tria;
+  make_mesh (tria);
+  
+  hp::FECollection<dim> fe;
+  fe.push_back (FE_Q<dim>(3));
+  fe.push_back (FE_Q<dim>(4));
+  fe.push_back (FE_Q<dim>(5));
+
+  hp::MappingCollection<dim> mapping_1;
+  mapping_1.push_back(MappingQ1<dim>());
+  mapping_1.push_back(MappingQ1<dim>());
+  mapping_1.push_back(MappingQ1<dim>());
+
+  hp::MappingCollection<dim> mapping_2;
+  mapping_2.push_back(MappingQ1<dim>());
+
+  hp::DoFHandler<dim> dof_handler (tria);
+
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell->set_active_fe_index (std::rand() % fe.size());
+    }
+
+  dof_handler.distribute_dofs (fe);
+
+                                   // test with two different functions: one
+                                   // that is exactly representable on the
+                                   // chosen finite element space, and one
+                                   // that isn't
+  for (unsigned int i=0; i<2; ++i)
+    {
+      static const MySquareFunction<dim>          function_1;
+      static const Functions::CosineFunction<dim> function_2;
+
+      const Function<dim> &
+        function = (i==0 ?
+                    static_cast<const Function<dim>&>(function_1) :
+                    static_cast<const Function<dim>&>(function_2));
+      
+      Vector<double> v (dof_handler.n_dofs());
+      VectorTools::interpolate (dof_handler, function, v);
+
+                                       // for the following points, check the
+                                       // function value, output it, and
+                                       // verify that the value retrieved from
+                                       // the interpolated function is close
+                                       // enough to that of the real function
+                                       //
+                                       // also verify that the actual value is
+                                       // roughly correct
+      Point<dim> p[3];
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          p[0][d] = 0;
+          p[1][d] = 0.5;
+          p[2][d] = 1./3.;
+        }
+      Vector<double> value(1);
+      for (unsigned int i=0; i<3; ++i)
+        {
+          VectorTools::point_value (mapping_1, dof_handler, v, p[i], value);
+          deallog << -value(0) << std::endl;
+
+          Assert (std::abs(value(0) - function.value(p[i])) < 1e-4,
+                   ExcInternalError());
+          VectorTools::point_value (mapping_2, dof_handler, v, p[i], value);
+          deallog << -value(0) << std::endl;
+
+          Assert (std::abs(value(0) - function.value(p[i])) < 1e-4,
+                  ExcInternalError());
+
+          const double scalar_value_1 = VectorTools::point_value (mapping_1, dof_handler, v, p[i]);
+          Assert (std::abs(value(0) - scalar_value_1) < 1e-4,
+              ExcInternalError());
+
+          const double scalar_value_2 = VectorTools::point_value (mapping_2, dof_handler, v, p[i]);
+          Assert (std::abs(value(0) - scalar_value_2) < 1e-4,
+              ExcInternalError());
+        }  
+    }
+  
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("point_value_hp_02/output");
+  deallog << std::setprecision (4);
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/bits/point_value_hp_02/cmp/generic b/tests/bits/point_value_hp_02/cmp/generic
new file mode 100644 (file)
index 0000000..2d85d56
--- /dev/null
@@ -0,0 +1,40 @@
+
+DEAL:1d::-1.000
+DEAL:1d::-1.000
+DEAL:1d::-1.250
+DEAL:1d::-1.250
+DEAL:1d::-1.111
+DEAL:1d::-1.111
+DEAL:1d::-1.000
+DEAL:1d::-1.000
+DEAL:1d::-0.7071
+DEAL:1d::-0.7071
+DEAL:1d::-0.8660
+DEAL:1d::-0.8660
+DEAL:1d::OK
+DEAL:2d::-1.000
+DEAL:2d::-1.000
+DEAL:2d::-1.500
+DEAL:2d::-1.500
+DEAL:2d::-1.222
+DEAL:2d::-1.222
+DEAL:2d::-1.000
+DEAL:2d::-1.000
+DEAL:2d::-0.5000
+DEAL:2d::-0.5000
+DEAL:2d::-0.7500
+DEAL:2d::-0.7500
+DEAL:2d::OK
+DEAL:3d::-1.000
+DEAL:3d::-1.000
+DEAL:3d::-1.750
+DEAL:3d::-1.750
+DEAL:3d::-1.333
+DEAL:3d::-1.333
+DEAL:3d::-1.000
+DEAL:3d::-1.000
+DEAL:3d::-0.3536
+DEAL:3d::-0.3536
+DEAL:3d::-0.6495
+DEAL:3d::-0.6495
+DEAL:3d::OK

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.