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
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
}
+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,
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>
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>&,
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>&);
<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
--- /dev/null
+//---------------------------- 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 ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 ();
+}
--- /dev/null
+
+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