]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix the same bug as before in several more places.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Jul 2012 13:40:33 +0000 (13:40 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Jul 2012 13:40:33 +0000 (13:40 +0000)
git-svn-id: https://svn.dealii.org/trunk@25682 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/numerics/fe_field_function.templates.h
tests/bits/fe_field_function_04_vector.cc
tests/bits/fe_field_function_05_vector.cc
tests/bits/fe_field_function_06_vector.cc [new file with mode: 0644]
tests/bits/fe_field_function_06_vector/cmp/generic [new file with mode: 0644]
tests/bits/fe_field_function_07_vector.cc [new file with mode: 0644]
tests/bits/fe_field_function_07_vector/cmp/generic [new file with mode: 0644]
tests/bits/fe_field_function_08_vector.cc [new file with mode: 0644]
tests/bits/fe_field_function_08_vector/cmp/generic [new file with mode: 0644]

index f438cafd12b987b82f845ff9fae815cc355d5fea..500f61e77024e3dae24ec91f2ad53c04e3622c96 100644 (file)
@@ -61,10 +61,10 @@ namespace Functions
     typename DH::active_cell_iterator cell = cell_hint.get();
     if (cell == dh->end())
       cell = dh->begin_active();
-    Point<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
 
-                                     // Check if we already have all we need
-    if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+    boost::optional<Point<dim> >
+      qp = get_reference_coordinates (cell, p);
+    if (!qp)
       {
         const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
           = GridTools::find_active_cell_around_point (mapping, *dh, p);
@@ -75,7 +75,7 @@ namespace Functions
     cell_hint.get() = cell;
 
                                      // Now we can find out about the point
-    Quadrature<dim> quad(qp);
+    Quadrature<dim> quad(qp.get());
     FEValues<dim> fe_v(mapping, cell->get_fe(), quad,
                        update_values);
     fe_v.reinit(cell);
@@ -97,21 +97,22 @@ namespace Functions
   }
 
 
+
   template <int dim, typename DH, typename VECTOR>
   void
-  FEFieldFunction<dim, DH, VECTOR>::vector_gradient
-  (const Point<dim> &p,
-   std::vector<Tensor<1,dim> > &gradients) const
+  FEFieldFunction<dim, DH, VECTOR>::
+  vector_gradient (const Point<dim>            &p,
+                  std::vector<Tensor<1,dim> > &gradients) const
   {
     Assert (gradients.size() == n_components,
             ExcDimensionMismatch(gradients.size(), n_components));
     typename DH::active_cell_iterator cell = cell_hint.get();
     if(cell == dh->end())
       cell = dh->begin_active();
-    Point<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
 
-                                     // Check if we already have all we need
-    if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+    boost::optional<Point<dim> >
+      qp = get_reference_coordinates (cell, p);
+    if (!qp)
       {
         std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
           = GridTools::find_active_cell_around_point (mapping, *dh, p);
@@ -122,7 +123,7 @@ namespace Functions
     cell_hint.get() = cell;
 
                                      // Now we can find out about the point
-    Quadrature<dim> quad(qp);
+    Quadrature<dim> quad(qp.get());
     FEValues<dim> fe_v(mapping, cell->get_fe(), quad,
                        update_gradients);
     fe_v.reinit(cell);
@@ -135,8 +136,10 @@ namespace Functions
 
 
   template <int dim, typename DH, typename VECTOR>
-  Tensor<1,dim> FEFieldFunction<dim, DH, VECTOR>::gradient
-  (const Point<dim>   &p, unsigned int comp) const
+  Tensor<1,dim>
+  FEFieldFunction<dim, DH, VECTOR>::
+  gradient (const Point<dim>   &p,
+           const unsigned int comp) const
   {
     std::vector<Tensor<1,dim> > grads(n_components);
     vector_gradient(p, grads);
@@ -146,18 +149,20 @@ namespace Functions
 
 
   template <int dim, typename DH, typename VECTOR>
-  void FEFieldFunction<dim, DH, VECTOR>::vector_laplacian
-  (const Point<dim> &p,Vector<double>   &values) const
+  void
+  FEFieldFunction<dim, DH, VECTOR>::
+  vector_laplacian (const Point<dim> &p,
+                   Vector<double>   &values) const
   {
     Assert (values.size() == n_components,
             ExcDimensionMismatch(values.size(), n_components));
     typename DH::active_cell_iterator cell = cell_hint.get();
     if (cell == dh->end())
       cell = dh->begin_active();
-    Point<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
 
-                                     // Check if we already have all we need
-    if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+    boost::optional<Point<dim> >
+      qp = get_reference_coordinates (cell, p);
+    if (!qp)
       {
         const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
           = GridTools::find_active_cell_around_point (mapping, *dh, p);
@@ -168,7 +173,7 @@ namespace Functions
     cell_hint.get() = cell;
 
                                      // Now we can find out about the point
-    Quadrature<dim> quad(qp);
+    Quadrature<dim> quad(qp.get());
     FEValues<dim> fe_v(mapping, cell->get_fe(), quad,
                        update_hessians);
     fe_v.reinit(cell);
index 3586f30fc609c739364a4c72790258d55ecdcb53..4e2c42ba5d48470f8fd6b9e580e26bb69af883a4 100644 (file)
@@ -82,7 +82,7 @@ void test()
       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());
+  std::vector<Vector<double> > m (points.size(), Vector<double>(2));
   fe_function.vector_value_list (points, m);
 
   for (unsigned int i=0; i<m.size(); ++i)
index 5a1dde289cd22ac01faba03bd7337a592df37ba0..1cca11c83de8849d9933fbb974f9d9c263cbc909 100644 (file)
@@ -94,7 +94,7 @@ void test()
                                4,
                                &dof_handler));
 
-  std::vector<Vector<double> > m (points.size());
+  std::vector<Vector<double> > m (points.size(), Vector<double>(2));
   fe_function.vector_value_list (points, m);
 
   for (unsigned int i=0; i<m.size(); ++i)
diff --git a/tests/bits/fe_field_function_06_vector.cc b/tests/bits/fe_field_function_06_vector.cc
new file mode 100644 (file)
index 0000000..79d334a
--- /dev/null
@@ -0,0 +1,125 @@
+//-------------------------------------------------------
+//    $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 yet another variant of _04 and _05 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);
+
+                                  // 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
+  Point<dim> point(-0.27999999999999992, -0.62999999999999989);
+  fe_function.set_active_cell (typename DoFHandler<dim>::active_cell_iterator
+                              (&triangulation,
+                               1,
+                               4,
+                               &dof_handler));
+
+  Vector<double> m(2);
+  fe_function.vector_value (point, m);
+
+  {
+    Assert (std::fabs(m(0) - point.square())
+           <
+           1e-10 * std::fabs(m(0) + point.square()),
+           ExcInternalError());
+
+    Assert (std::fabs(m(1))
+           <
+           1e-10,
+           ExcInternalError());
+  }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile("fe_field_function_06_vector/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<2>();
+
+  return 0;
+}
+
diff --git a/tests/bits/fe_field_function_06_vector/cmp/generic b/tests/bits/fe_field_function_06_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_07_vector.cc b/tests/bits/fe_field_function_07_vector.cc
new file mode 100644 (file)
index 0000000..55dd72d
--- /dev/null
@@ -0,0 +1,135 @@
+//-------------------------------------------------------
+//    $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 yet another variant of _04 and _05 and _06 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] = p[0];
+       v[1] = 0;
+      }
+
+    virtual void vector_gradient (const Point<dim> &p,
+                                 std::vector<Tensor<1,dim> > &v) const
+      {
+       v[0] = 0;
+       v[1] = 0;
+       v[0][0] = 1;
+      }
+};
+
+
+
+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);
+
+                                  // 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
+  Point<dim> point(-0.27999999999999992, -0.62999999999999989);
+  fe_function.set_active_cell (typename DoFHandler<dim>::active_cell_iterator
+                              (&triangulation,
+                               1,
+                               4,
+                               &dof_handler));
+
+  std::vector<Tensor<1,dim> > m(2);
+  fe_function.vector_gradient (point, m);
+
+  Assert (std::fabs(m[0][0] - 1)
+         <
+         1e-10 * std::fabs(m[0][0] + 1),
+         ExcInternalError());
+  Assert (std::fabs(m[0][1])
+         <
+         1e-10,
+         ExcInternalError());
+  Assert (m[1].norm()
+         <
+         1e-10,
+         ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile("fe_field_function_07_vector/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<2>();
+
+  return 0;
+}
+
diff --git a/tests/bits/fe_field_function_07_vector/cmp/generic b/tests/bits/fe_field_function_07_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_08_vector.cc b/tests/bits/fe_field_function_08_vector.cc
new file mode 100644 (file)
index 0000000..52ca8a0
--- /dev/null
@@ -0,0 +1,134 @@
+//-------------------------------------------------------
+//    $Id: fe_field_function_08_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 yet another variant of _04 through _07 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;
+       for (unsigned int i=0; i<dim; ++i)
+         v[0] += p[i]*p[i]*p[i]*p[i];
+      }
+
+    virtual void vector_laplacian (const Point<dim> &p,
+                                  Vector<double> &v) const
+      {
+       v = 0;
+       v[0] = p.square() * 4 * 3 * dim;
+      }
+};
+
+
+
+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> (4),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);
+
+                                  // 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
+  Point<dim> point(-0.27999999999999992, -0.62999999999999989);
+  fe_function.set_active_cell (typename DoFHandler<dim>::active_cell_iterator
+                              (&triangulation,
+                               1,
+                               4,
+                               &dof_handler));
+
+  Vector<double> m(2);
+  fe_function.vector_laplacian (point, m);
+
+  {
+    Assert (std::fabs(m(0) - point.square()*4*3)
+           <
+           1e-8 * std::fabs(m(0) + point.square()*4*3),
+           ExcInternalError());
+
+    Assert (std::fabs(m(1))
+           <
+           1e-10,
+           ExcInternalError());
+  }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile("fe_field_function_08_vector/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<2>();
+
+  return 0;
+}
+
diff --git a/tests/bits/fe_field_function_08_vector/cmp/generic b/tests/bits/fe_field_function_08_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.