From 049ba4f438afde3751c4fe22bb0f93c64a7cb07c Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 12 Jul 2004 14:09:55 +0000 Subject: [PATCH] Use a logarithmic algorithm instead of a linear on in point_difference. git-svn-id: https://svn.dealii.org/trunk@9511 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/numerics/vectors.cc | 84 ++++++++++------------ 1 file changed, 38 insertions(+), 46 deletions(-) diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 73ae7cb54c..08baddb511 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -13,23 +13,24 @@ #include -#include -#include -#include -#include -#include -#include #include -#include -#include -#include #include #include #include #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include #include +#include +#include #include #include @@ -1506,11 +1507,11 @@ VectorTools::integrate_difference (const DoFHandler &dof, template void -VectorTools::point_difference (const DoFHandler& dof, - const InVector& fe_function, - const Function& exact_function, - Vector& difference, - const Point& point) +VectorTools::point_difference (const DoFHandler &dof, + const InVector &fe_function, + const Function &exact_function, + Vector &difference, + const Point &point) { static const MappingQ1 mapping; const FiniteElement& fe = dof.get_fe(); @@ -1518,43 +1519,34 @@ VectorTools::point_difference (const DoFHandler& dof, Assert(difference.size() == fe.n_components(), ExcDimensionMismatch(difference.size(), fe.n_components())); - typename DoFHandler::active_cell_iterator - cell = dof.begin_active(); + // 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 - endc = dof.end(); + 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 (std::vector > (1, unit_point), + std::vector (1, 1.)); + 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); -//TODO: replace this linear search over all cells by using a logarithmic: first find teh correct cell on the coarse grid, then loop recursively over its children - for (;cell != endc; ++cell) - { - if (cell->point_inside(point)) - { - Point unit_point = mapping.transform_real_to_unit_cell( - cell, point); - for (unsigned int i=0; i > quadrature_point(1, unit_point); - std::vector weight(1, 1.); - Quadrature quadrature(quadrature_point, weight); - FEValues fe_values( - mapping, fe, quadrature, UpdateFlags(update_values)); - fe_values.reinit(cell); - fe_values.get_function_values(fe_function, u_value); - break; - } - } if (fe.n_components() == 1) - { - const double value = exact_function.value(point); - difference(0) = value; - } else { - exact_function.vector_value(point, difference); - } - for (unsigned int i=0;i