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

    +
  1. + 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) +

  2. 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) +

  3. 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