]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Several more tests.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Jul 2012 12:48:34 +0000 (12:48 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Jul 2012 12:48:34 +0000 (12:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@25679 0785d39b-7218-0410-832d-ea1e28bc413d

tests/bits/fe_field_function_01_vector.cc [new file with mode: 0644]
tests/bits/fe_field_function_01_vector/cmp/generic [new file with mode: 0644]
tests/bits/fe_field_function_02_vector.cc [new file with mode: 0644]
tests/bits/fe_field_function_02_vector/cmp/generic [new file with mode: 0644]
tests/bits/fe_field_function_03_vector.cc [new file with mode: 0644]
tests/bits/fe_field_function_03_vector/cmp/generic [new file with mode: 0644]
tests/bits/fe_field_function_04_vector.cc [new file with mode: 0644]
tests/bits/fe_field_function_04_vector/cmp/generic [new file with mode: 0644]
tests/bits/fe_field_function_05_vector.cc [new file with mode: 0644]
tests/bits/fe_field_function_05_vector/cmp/generic [new file with mode: 0644]

diff --git a/tests/bits/fe_field_function_01_vector.cc b/tests/bits/fe_field_function_01_vector.cc
new file mode 100644 (file)
index 0000000..391a699
--- /dev/null
@@ -0,0 +1,110 @@
+//----------------------------  template.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005, 2011, 2012 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.
+//
+//----------------------------  template.cc  ---------------------------
+
+
+// Test the FEFieldFunction class
+
+#include "../tests.h"
+#include <fstream>
+
+// all include files you need here
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/vectors.h>
+
+
+template <int dim>
+class F : public Function<dim>
+{
+  public:
+    F() : Function<dim>(2) {}
+    virtual void vector_value (const Point<dim> &p,
+                              Vector<double> &v) const
+      {
+       v = 0;
+       v[0] = p.square();
+      }
+};
+
+double abs_zero(double a)
+{
+  if( std::abs(a) < 1e-10)
+    return 0;
+  else
+    return a;
+}
+
+template <int dim>
+void test() {
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(8/dim);
+
+  FESystem<dim> fe(FE_Q<dim> (2),2);
+  DoFHandler<dim> dh(tria);
+
+  dh.distribute_dofs(fe);
+  deallog << "Ndofs :" << dh.n_dofs() << std::endl;
+
+
+
+  Vector<double> v1(dh.n_dofs()), v2(dh.n_dofs());
+
+  VectorTools::interpolate(dh, F<dim>(), v1);
+  deallog << "V norm: " << v1.l2_norm() << std::endl;
+
+  Functions::FEFieldFunction<dim, DoFHandler<dim>, Vector<double> >
+    fef(dh, v1);
+
+  VectorTools::interpolate(dh, fef, v2);
+
+  v2.add(-1, v1);
+  deallog << "Interpolation error: " << abs_zero(v2.l2_norm())
+         << std::endl;
+
+  Vector<double> error(tria.n_active_cells());
+  QGauss<dim> quad(2);
+  VectorTools::integrate_difference(dh, v1, F<dim>(),
+                                   error, quad,
+                                   VectorTools::L2_norm);
+  deallog << "L2 Interpolation error: "
+         << abs_zero(error.l2_norm()) << std::endl;
+  error = 0;
+  VectorTools::integrate_difference(dh, v1, fef,
+                                   error, quad,
+                                   VectorTools::L2_norm);
+  deallog << "L2 Interpolation error with fef: "
+         << abs_zero(error.l2_norm()) << std::endl;
+
+}
+
+int main ()
+{
+  std::ofstream logfile("fe_field_function_01_vector/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<1>();
+  test<2>();
+  test<3>();
+
+  return 0;
+}
+
diff --git a/tests/bits/fe_field_function_01_vector/cmp/generic b/tests/bits/fe_field_function_01_vector/cmp/generic
new file mode 100644 (file)
index 0000000..3115fec
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL::Ndofs :1026
+DEAL::V norm: 10.1440
+DEAL::Interpolation error: 0.00000
+DEAL::L2 Interpolation error: 0.00000
+DEAL::L2 Interpolation error with fef: 0.00000
+DEAL::Ndofs :2178
+DEAL::V norm: 26.5652
+DEAL::Interpolation error: 0.00000
+DEAL::L2 Interpolation error: 0.00000
+DEAL::L2 Interpolation error with fef: 0.00000
+DEAL::Ndofs :1458
+DEAL::V norm: 32.6964
+DEAL::Interpolation error: 0.00000
+DEAL::L2 Interpolation error: 0.00000
+DEAL::L2 Interpolation error with fef: 0.00000
diff --git a/tests/bits/fe_field_function_02_vector.cc b/tests/bits/fe_field_function_02_vector.cc
new file mode 100644 (file)
index 0000000..c5db68e
--- /dev/null
@@ -0,0 +1,109 @@
+//----------------------------  template.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005, 2011, 2012 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.
+//
+//----------------------------  template.cc  ---------------------------
+
+
+// Test the FEFieldFunction class. This class was not thread-safe at one point
+// because it keeps a cache on the side that was invalidated when different
+// threads kept pouncing on it. Patrick Sodre wrote a fix for that that's
+// based on a thread-local variable instead of the single cache.
+
+#include "../tests.h"
+#include <fstream>
+
+// all include files you need here
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/vectors.h>
+#include <deal.II/fe/fe_system.h>
+
+
+template <int dim>
+class F : public Function<dim>
+{
+  public:
+    F() : Function<dim>(2) {}
+    virtual void vector_value (const Point<dim> &p,
+                              Vector<double> &v) const
+      {
+       v = 0;
+       v[0] = p.square();
+      }
+};
+
+
+double abs_zero(double a)
+{
+  if( std::abs(a) < 1e-10)
+    return 0;
+  else
+    return a;
+}
+
+template <int dim>
+void test() {
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(8/dim);
+
+  FESystem<dim> fe(FE_Q<dim> (2),2);
+  DoFHandler<dim> dh(tria);
+
+  dh.distribute_dofs(fe);
+  deallog << "Ndofs :" << dh.n_dofs() << std::endl;
+
+  F<dim> ff;
+
+  Vector<double> v1(dh.n_dofs()), v2(dh.n_dofs());
+
+  VectorTools::interpolate(dh, ff, v1);
+  deallog << "V norm: " << v1.l2_norm() << std::endl;
+
+  Functions::FEFieldFunction<dim, DoFHandler<dim>, Vector<double> >
+    fef(dh, v1);
+
+                                  // project the discrete function fef back
+                                  // onto the finite element space. this
+                                  // should be the identity operation and
+                                  // consequently subtracting the vector from
+                                  // itself should yield zero
+  {
+    ConstraintMatrix cm;
+    cm.close();
+    VectorTools::project(dh, cm, QGauss<dim>(3), fef, v2);
+  }
+  v2.add(-1, v1);
+
+  deallog << "Projection error: " << abs_zero(v2.l2_norm())
+         << std::endl;
+  Assert (v2.l2_norm() < 1e-10, ExcInternalError());
+}
+
+int main ()
+{
+  std::ofstream logfile("fe_field_function_02_vector/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<1>();
+  test<2>();
+  test<3>();
+
+  return 0;
+}
+
diff --git a/tests/bits/fe_field_function_02_vector/cmp/generic b/tests/bits/fe_field_function_02_vector/cmp/generic
new file mode 100644 (file)
index 0000000..021c0ab
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Ndofs :1026
+DEAL::V norm: 10.1440
+DEAL::Projection error: 0.00000
+DEAL::Ndofs :2178
+DEAL::V norm: 26.5652
+DEAL::Projection error: 0.00000
+DEAL::Ndofs :1458
+DEAL::V norm: 32.6964
+DEAL::Projection error: 0.00000
diff --git a/tests/bits/fe_field_function_03_vector.cc b/tests/bits/fe_field_function_03_vector.cc
new file mode 100644 (file)
index 0000000..834205d
--- /dev/null
@@ -0,0 +1,88 @@
+//-------------------------------------------------------
+//    $Id: fe_field_function_03.cc $
+//    Version: $Name$
+//
+//    Copyright (C) 2005, 2011 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.
+//
+//-------------------------------------------------------
+
+
+// Test the functionality of the laplacian in the FEFieldFunction class. 
+
+#include "../tests.h"
+#include <fstream>
+
+// all include files you need here
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/vectors.h>
+
+template <int dim>
+void test() {
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(9/dim);
+
+  FE_Q<dim> fe(2);
+  DoFHandler<dim> dh(tria);
+
+  dh.distribute_dofs(fe);
+  deallog << "Ndofs :" << dh.n_dofs() << std::endl;
+
+  Functions::CosineFunction<dim> ff;
+
+  Vector<double> v1(dh.n_dofs()), v2(dh.n_dofs());
+
+  VectorTools::interpolate(dh, ff, v1);
+  deallog << "V norm: " << v1.l2_norm() << std::endl;
+
+  Functions::FEFieldFunction<dim, DoFHandler<dim>, Vector<double> >
+    fef(dh, v1);
+
+       //create the origin
+       Point<dim> p;
+       //compute the error of the laplacian in this point
+  deallog << "Value of the laplacian in 0:" << std::endl;
+  deallog << "correct value: " << ff.laplacian(p) <<", approximation: "<< fef.laplacian(p) << std::endl;
+  
+  //now we want to test the list version
+  Point<dim> p1 = Point<dim>::unit_vector(0);
+  p1 = p1 * 0.5;
+  Point<dim> p2 = p1 * 0.5;
+  std::vector<Point<dim> > vec;
+  vec.push_back(p1);
+  vec.push_back(p2);
+  std::vector<double> values_c(2);
+  std::vector<double> values_a(2);
+  
+  //get the laplacians at these two points
+  ff.laplacian_list(vec, values_c);
+  fef.laplacian_list(vec, values_a);
+    deallog << "Value of the laplacian in 0.5*e1 and 0.25 * e1:" << std::endl;
+    deallog << " correct values: " <<values_c[0] << " " << values_c[1]  <<", approximations: "<<values_a[0] << " " << values_a[1] << std::endl;
+}
+
+int main ()
+{
+  std::ofstream logfile("fe_field_function_03/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<1>();
+  test<2>();
+  test<3>();
+
+  return 0;
+}
+
diff --git a/tests/bits/fe_field_function_03_vector/cmp/generic b/tests/bits/fe_field_function_03_vector/cmp/generic
new file mode 100644 (file)
index 0000000..cd95fc4
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL::Ndofs :1025
+DEAL::V norm: 22.6385
+DEAL::Value of the laplacian in 0:
+DEAL::correct value: -2.46740, approximation: -2.46743
+DEAL::Value of the laplacian in 0.5*e1 and 0.25 * e1:
+DEAL:: correct values: -1.74472 -2.27958, approximations: -1.74739 -2.28103
+DEAL::Ndofs :1089
+DEAL::V norm: 16.5000
+DEAL::Value of the laplacian in 0:
+DEAL::correct value: -4.93480, approximation: -4.92787
+DEAL::Value of the laplacian in 0.5*e1 and 0.25 * e1:
+DEAL:: correct values: -3.48943 -4.55916, approximations: -3.57012 -4.59908
+DEAL::Ndofs :4913
+DEAL::V norm: 24.7815
+DEAL::Value of the laplacian in 0:
+DEAL::correct value: -7.40220, approximation: -7.36064
+DEAL::Value of the laplacian in 0.5*e1 and 0.25 * e1:
+DEAL:: correct values: -5.23415 -6.83874, approximations: -5.37564 -6.89283
diff --git a/tests/bits/fe_field_function_04_vector.cc b/tests/bits/fe_field_function_04_vector.cc
new file mode 100644 (file)
index 0000000..3586f30
--- /dev/null
@@ -0,0 +1,115 @@
+//-------------------------------------------------------
+//    $Id: fe_field_function_04_vector.cc $
+//    Version: $Name$
+//
+//    Copyright (C) 2005, 2011, 2012 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.
+//
+//-------------------------------------------------------
+
+
+// FEFieldFunction ran into an assertion after
+// Mapping::transform_real_to_unit_cell started throwing exceptions
+// when it couldn't find the point on the reference cell that belongs
+// to a given point, rather than just silently giving up
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/numerics/vectors.h>
+#include <deal.II/fe/fe_system.h>
+
+
+template <int dim>
+class F : public Function<dim>
+{
+  public:
+    F() : Function<dim>(2) {}
+    virtual void vector_value (const Point<dim> &p,
+                              Vector<double> &v) const
+      {
+       v = 0;
+       v[0] = p.square();
+      }
+};
+
+
+template<int dim>
+void test()
+{
+  const HyperBallBoundary<dim> boundary_description;
+
+  Triangulation<dim>   triangulation;
+  GridGenerator::hyper_ball (triangulation);
+  triangulation.set_boundary (0, boundary_description);
+  triangulation.refine_global (1);
+
+  FESystem<dim> fe(FE_Q<dim> (2),2);
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  Vector<double> solution(dof_handler.n_dofs());
+  VectorTools::interpolate (dof_handler, F<dim>(), solution);
+
+  Functions::FEFieldFunction<2> fe_function (dof_handler, solution);
+  std::vector<Point<dim> > points;
+
+                                  // add a bunch of points. all
+                                  // points are inside the circle.
+                                  // the problem happens because we
+                                  // walk over a bunch of cells in
+                                  // the process of finding all of
+                                  // these points and then realize
+                                  // when we get to the one at the
+                                  // end that the coordinates for
+                                  // this point can't be found in the
+                                  // cell we have touched last (it's
+                                  // too far away from that cell, and
+                                  // the inverse mapping does not
+                                  // converge
+  for (unsigned int i=0; i<20; ++i)
+    for (unsigned int j=0; j<20; ++j)
+      points.push_back (Point<dim>(-0.7+i*0.07,-0.7+j*0.07));
+  points.push_back (Point<dim>(-0.27999999999999992, -0.62999999999999989));
+
+  std::vector<Vector<double> > m (points.size());
+  fe_function.vector_value_list (points, m);
+
+  for (unsigned int i=0; i<m.size(); ++i)
+    {
+      Assert (std::fabs(m[i](0) - points[i].square())
+             <
+             1e-10 * std::fabs(m[i](0) + points[i].square()),
+             ExcInternalError());
+
+      Assert (std::fabs(m[i](1))
+             <
+             1e-10,
+             ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile("fe_field_function_04_vector/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<2>();
+
+  return 0;
+}
+
diff --git a/tests/bits/fe_field_function_04_vector/cmp/generic b/tests/bits/fe_field_function_04_vector/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/bits/fe_field_function_05_vector.cc b/tests/bits/fe_field_function_05_vector.cc
new file mode 100644 (file)
index 0000000..5a1dde2
--- /dev/null
@@ -0,0 +1,127 @@
+//-------------------------------------------------------
+//    $Id: fe_field_function_05_vector.cc $
+//    Version: $Name$
+//
+//    Copyright (C) 2005, 2011, 2012 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.
+//
+//-------------------------------------------------------
+
+
+// FEFieldFunction ran into an assertion after
+// Mapping::transform_real_to_unit_cell started throwing exceptions
+// when it couldn't find the point on the reference cell that belongs
+// to a given point, rather than just silently giving up
+//
+// this is a variant of _04 found by inspecting the code in
+// FEFieldFunction but the (same) exception is triggered in a
+// different place
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/numerics/vectors.h>
+#include <deal.II/fe/fe_system.h>
+
+
+template <int dim>
+class F : public Function<dim>
+{
+  public:
+    F() : Function<dim>(2) {}
+    virtual void vector_value (const Point<dim> &p,
+                              Vector<double> &v) const
+      {
+       v = 0;
+       v[0] = p.square();
+      }
+};
+
+
+
+template<int dim>
+void test()
+{
+  const HyperBallBoundary<dim> boundary_description;
+
+  Triangulation<dim>   triangulation;
+  GridGenerator::hyper_ball (triangulation);
+  triangulation.set_boundary (0, boundary_description);
+  triangulation.refine_global (1);
+
+  FESystem<dim> fe(FE_Q<dim> (2),2);
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+                                  // interpolate a quadratic function
+                                  // into the space; this function
+                                  // can be represented exactly, so
+                                  // that we can compare again later
+  Vector<double> solution(dof_handler.n_dofs());
+  VectorTools::interpolate (dof_handler, F<dim>(), solution);
+
+  Functions::FEFieldFunction<2> fe_function (dof_handler, solution);
+  std::vector<Point<dim> > points;
+
+                                  // add only one points but also set
+                                  // the active cell to one that
+                                  // doesn't contain the current
+                                  // point.  the problem happens
+                                  // because we walk over a bunch of
+                                  // cells in the process of finding
+                                  // all of these points and then
+                                  // realize when we get to the one
+                                  // at the end that the coordinates
+                                  // for this point can't be found in
+                                  // the cell we have touched last
+                                  // (it's too far away from that
+                                  // cell, and the inverse mapping
+                                  // does not converge
+  points.push_back (Point<dim>(-0.27999999999999992, -0.62999999999999989));
+  fe_function.set_active_cell (typename DoFHandler<dim>::active_cell_iterator
+                              (&triangulation,
+                               1,
+                               4,
+                               &dof_handler));
+
+  std::vector<Vector<double> > m (points.size());
+  fe_function.vector_value_list (points, m);
+
+  for (unsigned int i=0; i<m.size(); ++i)
+    {
+      Assert (std::fabs(m[i](0) - points[i].square())
+             <
+             1e-10 * std::fabs(m[i](0) + points[i].square()),
+             ExcInternalError());
+
+      Assert (std::fabs(m[i](1))
+             <
+             1e-10,
+             ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile("fe_field_function_05_vector/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<2>();
+
+  return 0;
+}
+
diff --git a/tests/bits/fe_field_function_05_vector/cmp/generic b/tests/bits/fe_field_function_05_vector/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::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.