]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New functions to evaluate a finite element solution at an arbitrary point.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 6 Mar 2006 21:05:51 +0000 (21:05 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 6 Mar 2006 21:05:51 +0000 (21:05 +0000)
git-svn-id: https://svn.dealii.org/trunk@12542 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/Makefile
tests/bits/point_value_01.cc [new file with mode: 0644]
tests/bits/point_value_01/cmp/generic [new file with mode: 0644]

index c32620f264b08fa4b383fca55ed6037eb0e83441..edf74c1b3c71036103e7517ce7829575b6afce2b 100644 (file)
@@ -938,7 +938,7 @@ class VectorTools
                                      * simple test for checking
                                      * whether a point is in a cell,
                                      * it is only implemented for
-                                     * Q1-mapping yet.
+                                     * Q1-mapping at this time.
                                      */
     template <int dim, class InVector>
     static void point_difference (const DoFHandler<dim>& dof,
@@ -947,6 +947,39 @@ class VectorTools
                                  Vector<double>&        difference,
                                  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.
+                                     */
+    template <int dim, class InVector>
+    static
+    void
+    point_value (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.
+                                     */
+    template <int dim, class InVector>
+    static
+    double
+    point_value (const DoFHandler<dim> &dof,
+                const InVector        &fe_function,
+                const Point<dim>      &point);
+    
                                     /**
                                      * Subtract the (algebraic) mean value
                                      * from a vector. This function is most
index 29cc10c8e6dd6f7c224bf7a30b7eaade97830be4..30daa10fd6e13a21d70475a6ff0ab53b65ef1371 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005 by the deal.II authors
+//    Copyright (C) 2005, 2006 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -1624,8 +1624,7 @@ VectorTools::point_difference (const DoFHandler<dim> &dof,
   Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
           ExcInternalError());
 
-  const Quadrature<dim> quadrature (std::vector<Point<dim> > (1, unit_point),
-                                    std::vector<double> (1, 1.));
+  const Quadrature<dim> quadrature (unit_point);
   FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
   fe_values.reinit(cell);
 
@@ -1644,6 +1643,80 @@ VectorTools::point_difference (const DoFHandler<dim> &dof,
 }
 
 
+template <int dim, class InVector>
+void
+VectorTools::point_value (const DoFHandler<dim> &dof,
+                         const InVector        &fe_function,
+                         const Point<dim>      &point,
+                         Vector<double>        &value)
+{
+  static const MappingQ1<dim> mapping;
+  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
+  const typename DoFHandler<dim>::active_cell_iterator
+    cell = GridTools::find_active_cell_around_point (dof, point);
+
+  const Point<dim> unit_point
+    = mapping.transform_real_to_unit_cell(cell, point);
+  Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
+          ExcInternalError());
+
+  const Quadrature<dim> quadrature (unit_point);
+  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+  fe_values.reinit(cell);
+
+                                   // 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 DoFHandler<dim> &dof,
+                         const InVector        &fe_function,
+                         const Point<dim>      &point)
+{
+  static const MappingQ1<dim> mapping;
+  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
+  const typename DoFHandler<dim>::active_cell_iterator
+    cell = GridTools::find_active_cell_around_point (dof, point);
+
+  const Point<dim> unit_point
+    = mapping.transform_real_to_unit_cell(cell, point);
+  Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
+          ExcInternalError());
+
+  const Quadrature<dim> quadrature (unit_point);
+  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+  fe_values.reinit(cell);
+
+                                   // 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>
 double
index a2d1464b22049b17207bdfdb159ee67d2d3c249c..4faf6c45a153ea321b31859d1c1cb64448f0676c 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005 by the deal.II authors
+//    Copyright (C) 2005, 2006 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -88,6 +88,19 @@ void VectorTools::point_difference<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
 double VectorTools::compute_mean_value<deal_II_dimension>
 (const Mapping<deal_II_dimension>&,
index 4a2de2b8203fef9110bdff8cbe2ab5eeb6a41127..0b3c42caafede669cb8dbb849503aed166858a03 100644 (file)
@@ -495,14 +495,21 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       New: There are now functions <code class="class">VectorTools</code>::<code
+       class="member">point_value</code> that evaluate the value of a
+       finite element function at a given point inside the domain.
+       <br> 
+       (WB 2006/03/06)
+       </p>
 
   <li> <p> Improved: <code class="class">GridOut</code> now has functions for
-  declaring and parsing parameters in a <code
-  class="class">ParameterHandler</code>, like used in <code
-  class="class">DataOut_Interface</code> for a long time.
-  <br>
-  (GK 2006/03/06)
-  </p>
+       declaring and parsing parameters in a <code
+       class="class">ParameterHandler</code>, like used in <code
+       class="class">DataOut_Interface</code> for a long time.
+       <br>
+       (GK 2006/03/06)
+       </p>
 
   <li> <p>
        Improved: The <code class="class">Triangulation</code>, <code
index f4a3f99f28773f315a49a835ba5a18671ac32c21..f8a973c8b0106600914e153b1603b7a5b47a79e4 100644 (file)
@@ -1,6 +1,6 @@
 ############################################################
 # Makefile,v 1.15 2002/06/13 12:51:13 hartmann Exp
-# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
+# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
 ############################################################
 
 ############################################################
@@ -65,6 +65,7 @@ tests_x = geometry_info_* \
          cylinder_shell_* \
          compressed_sparsity_pattern_* \
          point_difference_* \
+         point_value_* \
          fe_q_3d* \
           rt_* \
          nedelec_* \
diff --git a/tests/bits/point_value_01.cc b/tests/bits/point_value_01.cc
new file mode 100644 (file)
index 0000000..6ced721
--- /dev/null
@@ -0,0 +1,168 @@
+//----------------------------  point_value_01.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_01.cc  ---------------------------
+
+// check VectorTools::point_value
+
+
+
+#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 <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);
+  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 (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 (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_01/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_01/cmp/generic b/tests/bits/point_value_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

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.