From 4eaaee0b66867a2ab338e7da440fb3f9a184823e Mon Sep 17 00:00:00 2001 From: rschulz Date: Thu, 11 May 2006 11:46:17 +0000 Subject: [PATCH] Changed VectorTools::point_difference and ::point_value to support the new cell finding routines also git-svn-id: https://svn.dealii.org/trunk@13095 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/vectors.h | 68 +++++++ .../include/numerics/vectors.templates.h | 112 +++++++++++ .../source/numerics/vectors.instance.h | 24 +++ deal.II/doc/news/changes.html | 9 + tests/bits/point_difference_02.cc | 174 ++++++++++++++++++ tests/bits/point_difference_02/cmp/generic | 40 ++++ tests/bits/point_value_02.cc | 170 +++++++++++++++++ tests/bits/point_value_02/cmp/generic | 22 +++ 8 files changed, 619 insertions(+) create mode 100644 tests/bits/point_difference_02.cc create mode 100644 tests/bits/point_difference_02/cmp/generic create mode 100644 tests/bits/point_value_02.cc create mode 100644 tests/bits/point_value_02/cmp/generic diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index f1eed46568..267aa2f7ea 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -978,6 +978,31 @@ class VectorTools const Point& 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 + static void point_difference (const Mapping &mapping, + const DoFHandler& dof, + const InVector& fe_function, + const Function& exact_solution, + Vector& difference, + const Point& point); + + /** * Evaluate a possibly * vector-valued finite element * function defined by the given @@ -1011,6 +1036,49 @@ class VectorTools 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. 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 + static + void + point_value (const Mapping &mapping, + 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. 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 + static + double + point_value (const Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Point &point); + + /** * Subtract the (algebraic) mean value * from a vector. This function is most * frequently used as a mean-value filter diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index 32e9174973..f23a7c8f4d 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -1808,6 +1808,48 @@ VectorTools::point_difference (const DoFHandler &dof, } +template +void +VectorTools::point_difference (const Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Function &exact_function, + Vector &difference, + const Point &point) +{ + const FiniteElement& 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::active_cell_iterator, Point > + cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + + Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); + GeometryInfo::project_to_unit_cell(cell_point.second); + + const Quadrature quadrature (cell_point.second); + FEValues 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 > u_value(1, Vector (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 void VectorTools::point_value (const DoFHandler &dof, @@ -1881,6 +1923,76 @@ VectorTools::point_value (const DoFHandler &dof, return u_value[0]; } +template +void +VectorTools::point_value (const Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Point &point, + Vector &value) +{ + 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 + std::pair::active_cell_iterator, Point > + cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + + Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); + + GeometryInfo::project_to_unit_cell(cell_point.second); + + const Quadrature quadrature (cell_point.second); + FEValues 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 > 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 Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Point &point) +{ + 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 + std::pair::active_cell_iterator, Point > + cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + + Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); + + GeometryInfo::project_to_unit_cell(cell_point.second); + + const Quadrature quadrature (cell_point.second); + FEValues 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 u_value(1); + fe_values.get_function_values(fe_function, u_value); + + return u_value[0]; +} template diff --git a/deal.II/deal.II/source/numerics/vectors.instance.h b/deal.II/deal.II/source/numerics/vectors.instance.h index 643dd8057b..e83be170f4 100644 --- a/deal.II/deal.II/source/numerics/vectors.instance.h +++ b/deal.II/deal.II/source/numerics/vectors.instance.h @@ -131,8 +131,31 @@ void VectorTools::point_difference ( Vector&, const Point&); +template +void VectorTools::point_difference ( + const Mapping&, + const DoFHandler&, + const VEC&, + const Function&, + 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 void VectorTools::point_value ( + const Mapping&, const DoFHandler&, const VEC&, const Point&, @@ -140,6 +163,7 @@ void VectorTools::point_value ( template double VectorTools::point_value ( + const Mapping&, const DoFHandler&, const VEC&, const Point&); diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index 32346a41a5..50f155819d 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -646,6 +646,15 @@ inconvenience this causes.

deal.II

    +
  1. + Improved: The functions VectorTools::point_value and VectorTools::point_difference now can also use arbitrary + mappings, using the new GridTools::find_active_cell_around_point algorithm. +
    + (Ralf B. Schulz, 2006/05/11) +

  2. New: In GridTools, several functions have been added. GridTools::find_closest_vertex searches for the vertex diff --git a/tests/bits/point_difference_02.cc b/tests/bits/point_difference_02.cc new file mode 100644 index 0000000000..6163cc0794 --- /dev/null +++ b/tests/bits/point_difference_02.cc @@ -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 +#include +#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 () +{ + MappingQ1 mapping; + 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 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(), + 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 index 0000000000..b83ea63bab --- /dev/null +++ b/tests/bits/point_difference_02/cmp/generic @@ -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 index 0000000000..073d324e18 --- /dev/null +++ b/tests/bits/point_value_02.cc @@ -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 +#include +#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); + MappingQ1 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 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 (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 index 0000000000..763c6337fe --- /dev/null +++ b/tests/bits/point_value_02/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