From ca97fb3bc1ca817fe4cf129a390705b7c2daa610 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth
Date: Mon, 6 Mar 2006 21:05:51 +0000
Subject: [PATCH] New functions to evaluate a finite element solution at an
arbitrary point.
git-svn-id: https://svn.dealii.org/trunk@12542 0785d39b-7218-0410-832d-ea1e28bc413d
---
deal.II/deal.II/include/numerics/vectors.h | 35 +++-
.../include/numerics/vectors.templates.h | 79 +++++++-
.../source/numerics/vectors.instance.h | 15 +-
deal.II/doc/news/changes.html | 19 +-
tests/bits/Makefile | 3 +-
tests/bits/point_value_01.cc | 168 ++++++++++++++++++
tests/bits/point_value_01/cmp/generic | 22 +++
7 files changed, 329 insertions(+), 12 deletions(-)
create mode 100644 tests/bits/point_value_01.cc
create mode 100644 tests/bits/point_value_01/cmp/generic
diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h
index c32620f264..edf74c1b3c 100644
--- a/deal.II/deal.II/include/numerics/vectors.h
+++ b/deal.II/deal.II/include/numerics/vectors.h
@@ -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
static void point_difference (const DoFHandler& dof,
@@ -947,6 +947,39 @@ class VectorTools
Vector& difference,
const Point& 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
+ static
+ void
+ point_value (const DoFHandler &dof,
+ const InVector &fe_function,
+ const Point &point,
+ Vector &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
+ static
+ double
+ point_value (const DoFHandler &dof,
+ const InVector &fe_function,
+ const Point &point);
+
/**
* Subtract the (algebraic) mean value
* from a vector. This function is most
diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h
index 29cc10c8e6..30daa10fd6 100644
--- a/deal.II/deal.II/include/numerics/vectors.templates.h
+++ b/deal.II/deal.II/include/numerics/vectors.templates.h
@@ -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 &dof,
Assert (GeometryInfo::is_inside_unit_cell (unit_point),
ExcInternalError());
- const Quadrature quadrature (std::vector > (1, unit_point),
- std::vector (1, 1.));
+ const Quadrature quadrature (unit_point);
FEValues fe_values(mapping, fe, quadrature, update_values);
fe_values.reinit(cell);
@@ -1644,6 +1643,80 @@ VectorTools::point_difference (const DoFHandler &dof,
}
+template
+void
+VectorTools::point_value (const DoFHandler &dof,
+ const InVector &fe_function,
+ const Point &point,
+ Vector &value)
+{
+ static const MappingQ1 mapping;
+ const FiniteElement& 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::active_cell_iterator
+ cell = GridTools::find_active_cell_around_point (dof, point);
+
+ const Point unit_point
+ = mapping.transform_real_to_unit_cell(cell, point);
+ Assert (GeometryInfo::is_inside_unit_cell (unit_point),
+ ExcInternalError());
+
+ const Quadrature quadrature (unit_point);
+ FEValues 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 > u_value(1, Vector (fe.n_components()));
+ fe_values.get_function_values(fe_function, u_value);
+
+ value = u_value[0];
+}
+
+
+
+template
+double
+VectorTools::point_value (const DoFHandler &dof,
+ const InVector &fe_function,
+ const Point &point)
+{
+ static const MappingQ1 mapping;
+ const FiniteElement& 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::active_cell_iterator
+ cell = GridTools::find_active_cell_around_point (dof, point);
+
+ const Point unit_point
+ = mapping.transform_real_to_unit_cell(cell, point);
+ Assert (GeometryInfo::is_inside_unit_cell (unit_point),
+ ExcInternalError());
+
+ const Quadrature quadrature (unit_point);
+ FEValues 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 u_value(1);
+ fe_values.get_function_values(fe_function, u_value);
+
+ return u_value[0];
+}
+
+
template
double
diff --git a/deal.II/deal.II/source/numerics/vectors.instance.h b/deal.II/deal.II/source/numerics/vectors.instance.h
index a2d1464b22..4faf6c45a1 100644
--- a/deal.II/deal.II/source/numerics/vectors.instance.h
+++ b/deal.II/deal.II/source/numerics/vectors.instance.h
@@ -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 (
Vector&,
const Point&);
+template
+void VectorTools::point_value (
+ const DoFHandler&,
+ const VEC&,
+ const Point&,
+ Vector&);
+
+template
+double VectorTools::point_value (
+ const DoFHandler&,
+ const VEC&,
+ const Point&);
+
template
double VectorTools::compute_mean_value
(const Mapping&,
diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html
index 4a2de2b820..0b3c42caaf 100644
--- a/deal.II/doc/news/changes.html
+++ b/deal.II/doc/news/changes.html
@@ -495,14 +495,21 @@ inconvenience this causes.
deal.II
+ -
+ New: There are now functions VectorTools
::point_value
that evaluate the value of a
+ finite element function at a given point inside the domain.
+
+ (WB 2006/03/06)
+
-
Improved: GridOut
now has functions for
- declaring and parsing parameters in a ParameterHandler
, like used in DataOut_Interface
for a long time.
-
- (GK 2006/03/06)
-
+ declaring and parsing parameters in a ParameterHandler
, like used in DataOut_Interface
for a long time.
+
+ (GK 2006/03/06)
+
Improved: The Triangulation
,
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+
+template
+class MySquareFunction : public Function
+{
+ public:
+ virtual double value (const Point &p,
+ const unsigned int component) const
+ { return (component+1)*p.square()+1; };
+
+ virtual void vector_value (const Point &p,
+ Vector &values) const
+ { values(0) = value(p,0); };
+};
+
+
+template
+class MyExpFunction : public Function
+{
+ public:
+ virtual double value (const Point &p,
+ const unsigned int component) const
+ { return std::exp (p(0)); };
+
+ virtual void vector_value (const Point &p,
+ Vector &values) const
+ { values(0) = value(p,0); };
+};
+
+
+
+template
+void make_mesh (Triangulation &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::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
+void
+check ()
+{
+ Triangulation tria;
+ make_mesh (tria);
+
+ FE_Q element(3);
+ DoFHandler 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 function_1;
+ static const Functions::CosineFunction function_2;
+
+ const Function &
+ function = (i==0 ?
+ static_cast&>(function_1) :
+ static_cast&>(function_2));
+
+ Vector 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 p[3];
+ for (unsigned int d=0; d 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
index 0000000000..763c6337fe
--- /dev/null
+++ b/tests/bits/point_value_01/cmp/generic
@@ -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
--
2.39.5