]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Changed VectorTools::point_difference and ::point_value to support the new cell findi...
authorrschulz <rschulz@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 11 May 2006 11:46:17 +0000 (11:46 +0000)
committerrschulz <rschulz@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 11 May 2006 11:46:17 +0000 (11:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@13095 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/include/numerics/vectors.templates.h
deal.II/deal.II/source/numerics/vectors.instance.h
deal.II/doc/news/changes.html
tests/bits/point_difference_02.cc [new file with mode: 0644]
tests/bits/point_difference_02/cmp/generic [new file with mode: 0644]
tests/bits/point_value_02.cc [new file with mode: 0644]
tests/bits/point_value_02/cmp/generic [new file with mode: 0644]

index f1eed465683ef59f1442f095cdf12395d6ced389..267aa2f7eab6260a952fd0abf0c5952f20a8867b 100644 (file)
@@ -978,6 +978,31 @@ class VectorTools
                                  const Point<dim>&      point);
 
                                     /**
+                                     * Point error evaluation. Find
+                                     * the first cell containing the
+                                     * given point and compute the
+                                     * difference of a (possibly
+                                     * vector-valued) finite element
+                                     * function and a continuous
+                                     * function (with as many vector
+                                     * components as the finite
+                                     * element) at this point.
+                                     *
+                                      * This function additionally
+                                      * expects a mapping and uses
+                                      * a different algorithm to
+                                      * determine the cell surrounding
+                                      * the given point.
+                                     */
+    template <int dim, class InVector>
+    static void point_difference (const Mapping<dim>    &mapping,
+                                  const DoFHandler<dim>& dof,
+                                 const InVector&        fe_function,
+                                 const Function<dim>&   exact_solution,
+                                 Vector<double>&        difference,
+                                 const Point<dim>&      point);
+
+                                     /**
                                      * Evaluate a possibly
                                      * vector-valued finite element
                                      * function defined by the given
@@ -1011,6 +1036,49 @@ class VectorTools
                 const Point<dim>      &point);
     
                                     /**
+                                     * Evaluate a possibly
+                                     * vector-valued finite element
+                                     * function defined by the given
+                                     * DoFHandler and nodal vector at
+                                     * the given point, and return
+                                     * the (vector) value of this
+                                     * function through the last
+                                     * argument. This function expects
+                                      * additionally a mapping, and
+                                      * uses a different algorithm to
+                                      * determine the location of the
+                                      * point (see GridTools::find_active_cell_around_point).
+                                     */
+    template <int dim, class InVector>
+    static
+    void
+    point_value (const Mapping<dim>    &mapping,
+                 const DoFHandler<dim> &dof,
+                const InVector        &fe_function,
+                const Point<dim>      &point,
+                Vector<double>        &value);
+
+                                    /**
+                                     * Evaluate a scalar finite
+                                     * element function defined by
+                                     * the given DoFHandler and nodal
+                                     * vector at the given point, and
+                                     * return the value of this
+                                     * function. This function expects
+                                      * additionally a mapping, and
+                                      * uses a different algorithm to
+                                      * determine the location of the
+                                      * point (see GridTools::find_active_cell_around_point).
+                                     */
+    template <int dim, class InVector>
+    static
+    double
+    point_value (const Mapping<dim>    &mapping,
+                 const DoFHandler<dim> &dof,
+                const InVector        &fe_function,
+                const Point<dim>      &point);
+
+                                     /**
                                      * Subtract the (algebraic) mean value
                                      * from a vector. This function is most
                                      * frequently used as a mean-value filter
index 32e91749739020ffa8e852a7b7455475023671a2..f23a7c8f4d65bf31cf86274b581a6200cb216997 100644 (file)
@@ -1808,6 +1808,48 @@ VectorTools::point_difference (const DoFHandler<dim> &dof,
 }
 
 
+template <int dim, class InVector>
+void
+VectorTools::point_difference (const Mapping<dim>    &mapping,
+                               const DoFHandler<dim> &dof,
+                              const InVector        &fe_function,
+                              const Function<dim>   &exact_function,
+                              Vector<double>        &difference,
+                              const Point<dim>      &point)
+{
+  const FiniteElement<dim>& fe = dof.get_fe();
+
+  Assert(difference.size() == fe.n_components(),
+        ExcDimensionMismatch(difference.size(), fe.n_components()));
+
+                                   // first find the cell in which this point
+                                   // is, initialize a quadrature rule with
+                                   // it, and then a FEValues object
+  std::pair<typename DoFHandler<dim>::active_cell_iterator, Point<dim> >
+    cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
+
+  Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError());
+  GeometryInfo<dim>::project_to_unit_cell(cell_point.second);
+  
+  const Quadrature<dim> quadrature (cell_point.second);
+  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+  fe_values.reinit(cell_point.first);
+
+                                   // then use this to get at the values of
+                                   // the given fe_function at this point
+  std::vector<Vector<double> > u_value(1, Vector<double> (fe.n_components()));
+  fe_values.get_function_values(fe_function, u_value);
+
+  if (fe.n_components() == 1)
+    difference(0) = exact_function.value(point);
+  else
+    exact_function.vector_value(point, difference);
+    
+  for (unsigned int i=0; i<difference.size(); ++i)
+    difference(i) -= u_value[0](i);
+}
+
+
 template <int dim, class InVector>
 void
 VectorTools::point_value (const DoFHandler<dim> &dof,
@@ -1881,6 +1923,76 @@ VectorTools::point_value (const DoFHandler<dim> &dof,
   return u_value[0];
 }
 
+template <int dim, class InVector>
+void
+VectorTools::point_value (const Mapping<dim>    &mapping,
+                          const DoFHandler<dim> &dof,
+                         const InVector        &fe_function,
+                         const Point<dim>      &point,
+                         Vector<double>        &value)
+{
+  const FiniteElement<dim>& fe = dof.get_fe();
+
+  Assert(value.size() == fe.n_components(),
+        ExcDimensionMismatch(value.size(), fe.n_components()));
+
+                                   // first find the cell in which this point
+                                   // is, initialize a quadrature rule with
+                                   // it, and then a FEValues object
+  std::pair<typename DoFHandler<dim>::active_cell_iterator, Point<dim> >
+    cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
+
+  Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError());
+
+  GeometryInfo<dim>::project_to_unit_cell(cell_point.second);
+  
+  const Quadrature<dim> quadrature (cell_point.second);
+  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+  fe_values.reinit(cell_point.first);
+
+                                   // then use this to get at the values of
+                                   // the given fe_function at this point
+  std::vector<Vector<double> > u_value(1, Vector<double> (fe.n_components()));
+  fe_values.get_function_values(fe_function, u_value);
+
+  value = u_value[0];
+}
+
+
+
+template <int dim, class InVector>
+double
+VectorTools::point_value (const Mapping<dim>    &mapping,
+                          const DoFHandler<dim> &dof,
+                         const InVector        &fe_function,
+                         const Point<dim>      &point)
+{
+  const FiniteElement<dim>& fe = dof.get_fe();
+
+  Assert(fe.n_components() == 1,
+        ExcMessage ("Finite element is not scalar as is necessary for this function"));
+
+                                   // first find the cell in which this point
+                                   // is, initialize a quadrature rule with
+                                   // it, and then a FEValues object
+  std::pair<typename DoFHandler<dim>::active_cell_iterator, Point<dim> >
+    cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
+
+  Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError());
+
+  GeometryInfo<dim>::project_to_unit_cell(cell_point.second);
+  
+  const Quadrature<dim> quadrature (cell_point.second);
+  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+  fe_values.reinit(cell_point.first);
+
+                                   // then use this to get at the values of
+                                   // the given fe_function at this point
+  std::vector<double> u_value(1);
+  fe_values.get_function_values(fe_function, u_value);
+
+  return u_value[0];
+}
 
 
 template <int dim, class InVector>
index 643dd8057bd50133cdd377d8a5cf0017258f03f0..e83be170f4a91c496febdf3da0e669d02377c8b9 100644 (file)
@@ -131,8 +131,31 @@ void VectorTools::point_difference<deal_II_dimension> (
   Vector<double>&,
   const Point<deal_II_dimension>&);
 
+template
+void VectorTools::point_difference<deal_II_dimension> (
+  const Mapping<deal_II_dimension>&,
+  const DoFHandler<deal_II_dimension>&,
+  const VEC&,
+  const Function<deal_II_dimension>&,
+  Vector<double>&,
+  const Point<deal_II_dimension>&);
+
+template
+void VectorTools::point_value<deal_II_dimension> (
+  const DoFHandler<deal_II_dimension>&,
+  const VEC&,
+  const Point<deal_II_dimension>&,
+  Vector<double>&);
+
+template
+double VectorTools::point_value<deal_II_dimension> (
+  const DoFHandler<deal_II_dimension>&,
+  const VEC&,
+  const Point<deal_II_dimension>&);
+
 template
 void VectorTools::point_value<deal_II_dimension> (
+  const Mapping<deal_II_dimension>&,
   const DoFHandler<deal_II_dimension>&,
   const VEC&,
   const Point<deal_II_dimension>&,
@@ -140,6 +163,7 @@ void VectorTools::point_value<deal_II_dimension> (
 
 template
 double VectorTools::point_value<deal_II_dimension> (
+  const Mapping<deal_II_dimension>&,
   const DoFHandler<deal_II_dimension>&,
   const VEC&,
   const Point<deal_II_dimension>&);
index 32346a41a53a6c5bd1412c8e21ba3f18a5ab7123..50f155819d6a051eadb174d29d8acc391d1b8b18 100644 (file)
@@ -646,6 +646,15 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       Improved: The functions <code class="class">VectorTools</code>::<code
+       class="member">point_value</code> and <code class="class">VectorTools</code>::<code
+       class="member">point_difference</code> now can also use arbitrary
+       mappings, using the new <code class="class">GridTools</code>::<code
+       class="member">find_active_cell_around_point</code> algorithm.
+       <br> 
+       (Ralf B. Schulz, 2006/05/11)
+       </p>
   <li> <p>
        New: In <code>GridTools</code>, several functions have been added.
        <code>GridTools::find_closest_vertex</code> searches for the vertex
diff --git a/tests/bits/point_difference_02.cc b/tests/bits/point_difference_02.cc
new file mode 100644 (file)
index 0000000..6163cc0
--- /dev/null
@@ -0,0 +1,174 @@
+//----------------------------  point_difference_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004 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_difference_02.cc  ---------------------------
+
+// the same as point_difference_01, but for the alternative point_difference
+// algorithm
+
+// check that VectorTools::point_difference returns the same before and after
+// the change to a logarithmic algorithm
+
+
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/function_lib.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/mapping_q1.h>
+#include <numerics/vectors.h>
+
+#include <fstream>
+#include <cmath>
+#include <iostream>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+  public:
+    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:
+    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 ()
+{
+  MappingQ1<dim> mapping;
+  Triangulation<dim> tria;
+  make_mesh (tria);
+  
+  FE_Q<dim> element(3);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(element);
+
+                                   // 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.n_dofs());
+      VectorTools::interpolate (dof, 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> difference(1);
+      for (unsigned int i=0; i<3; ++i)
+        {
+          VectorTools::point_difference (mapping, dof, v, function, difference, p[i]);
+          deallog << difference(0) << std::endl;
+          Assert (difference(0) < 1e-4, ExcInternalError());
+
+          VectorTools::point_difference (mapping, dof, v, ZeroFunction<dim>(),
+                                         difference, p[i]);
+          deallog << difference(0) << std::endl;
+          Assert (std::abs(-difference(0) - function.value(p[i])) < 1e-4,
+                  ExcInternalError());
+        }  
+    }
+  
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("point_difference_02/output");
+  logfile.precision (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_difference_02/cmp/generic b/tests/bits/point_difference_02/cmp/generic
new file mode 100644 (file)
index 0000000..b83ea63
--- /dev/null
@@ -0,0 +1,40 @@
+
+DEAL:1d::0
+DEAL:1d::-1.000
+DEAL:1d::0
+DEAL:1d::-1.250
+DEAL:1d::0
+DEAL:1d::-1.111
+DEAL:1d::0
+DEAL:1d::-1.000
+DEAL:1d::0
+DEAL:1d::-0.7071
+DEAL:1d::0
+DEAL:1d::-0.8660
+DEAL:1d::OK
+DEAL:2d::0
+DEAL:2d::-1.000
+DEAL:2d::0
+DEAL:2d::-1.500
+DEAL:2d::-6.661e-16
+DEAL:2d::-1.222
+DEAL:2d::0
+DEAL:2d::-1.000
+DEAL:2d::0
+DEAL:2d::-0.5000
+DEAL:2d::-3.331e-16
+DEAL:2d::-0.7500
+DEAL:2d::OK
+DEAL:3d::0
+DEAL:3d::-1.000
+DEAL:3d::-2.220e-16
+DEAL:3d::-1.750
+DEAL:3d::0
+DEAL:3d::-1.333
+DEAL:3d::0
+DEAL:3d::-1.000
+DEAL:3d::1.110e-16
+DEAL:3d::-0.3536
+DEAL:3d::0
+DEAL:3d::-0.6495
+DEAL:3d::OK
diff --git a/tests/bits/point_value_02.cc b/tests/bits/point_value_02.cc
new file mode 100644 (file)
index 0000000..073d324
--- /dev/null
@@ -0,0 +1,170 @@
+//----------------------------  point_value_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 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.
+//
+//----------------------------  point_value_02.cc  ---------------------------
+
+// check VectorTools::point_value, alternative algorithm with mapping
+
+
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/function_lib.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/mapping_q1.h>
+#include <numerics/vectors.h>
+
+#include <fstream>
+#include <cmath>
+#include <iostream>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+  public:
+    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:
+    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);
+  
+  FE_Q<dim> element(3);
+  DoFHandler<dim> dof(tria);
+  MappingQ1<dim> mapping;
+  dof.distribute_dofs(element);
+
+                                   // 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.n_dofs());
+      VectorTools::interpolate (dof, 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, dof, 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 (mapping, dof, v, p[i]);
+          Assert (std::abs(value(0) - scalar_value) < 1e-4,
+                  ExcInternalError());
+        }  
+    }
+  
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("point_value_02/output");
+  logfile.precision (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_02/cmp/generic b/tests/bits/point_value_02/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

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.