From: bangerth Date: Wed, 4 Jul 2012 13:40:33 +0000 (+0000) Subject: Fix the same bug as before in several more places. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d92fe230ad9d007f224b39fbce056524e1d0bd47;p=dealii-svn.git Fix the same bug as before in several more places. git-svn-id: https://svn.dealii.org/trunk@25682 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/numerics/fe_field_function.templates.h b/deal.II/include/deal.II/numerics/fe_field_function.templates.h index f438cafd12..500f61e770 100644 --- a/deal.II/include/deal.II/numerics/fe_field_function.templates.h +++ b/deal.II/include/deal.II/numerics/fe_field_function.templates.h @@ -61,10 +61,10 @@ namespace Functions typename DH::active_cell_iterator cell = cell_hint.get(); if (cell == dh->end()) cell = dh->begin_active(); - Point qp = mapping.transform_real_to_unit_cell(cell, p); - // Check if we already have all we need - if (!GeometryInfo::is_inside_unit_cell(qp)) + boost::optional > + qp = get_reference_coordinates (cell, p); + if (!qp) { const std::pair > my_pair = GridTools::find_active_cell_around_point (mapping, *dh, p); @@ -75,7 +75,7 @@ namespace Functions cell_hint.get() = cell; // Now we can find out about the point - Quadrature quad(qp); + Quadrature quad(qp.get()); FEValues fe_v(mapping, cell->get_fe(), quad, update_values); fe_v.reinit(cell); @@ -97,21 +97,22 @@ namespace Functions } + template void - FEFieldFunction::vector_gradient - (const Point &p, - std::vector > &gradients) const + FEFieldFunction:: + vector_gradient (const Point &p, + std::vector > &gradients) const { Assert (gradients.size() == n_components, ExcDimensionMismatch(gradients.size(), n_components)); typename DH::active_cell_iterator cell = cell_hint.get(); if(cell == dh->end()) cell = dh->begin_active(); - Point qp = mapping.transform_real_to_unit_cell(cell, p); - // Check if we already have all we need - if (!GeometryInfo::is_inside_unit_cell(qp)) + boost::optional > + qp = get_reference_coordinates (cell, p); + if (!qp) { std::pair > my_pair = GridTools::find_active_cell_around_point (mapping, *dh, p); @@ -122,7 +123,7 @@ namespace Functions cell_hint.get() = cell; // Now we can find out about the point - Quadrature quad(qp); + Quadrature quad(qp.get()); FEValues fe_v(mapping, cell->get_fe(), quad, update_gradients); fe_v.reinit(cell); @@ -135,8 +136,10 @@ namespace Functions template - Tensor<1,dim> FEFieldFunction::gradient - (const Point &p, unsigned int comp) const + Tensor<1,dim> + FEFieldFunction:: + gradient (const Point &p, + const unsigned int comp) const { std::vector > grads(n_components); vector_gradient(p, grads); @@ -146,18 +149,20 @@ namespace Functions template - void FEFieldFunction::vector_laplacian - (const Point &p,Vector &values) const + void + FEFieldFunction:: + vector_laplacian (const Point &p, + Vector &values) const { Assert (values.size() == n_components, ExcDimensionMismatch(values.size(), n_components)); typename DH::active_cell_iterator cell = cell_hint.get(); if (cell == dh->end()) cell = dh->begin_active(); - Point qp = mapping.transform_real_to_unit_cell(cell, p); - // Check if we already have all we need - if (!GeometryInfo::is_inside_unit_cell(qp)) + boost::optional > + qp = get_reference_coordinates (cell, p); + if (!qp) { const std::pair > my_pair = GridTools::find_active_cell_around_point (mapping, *dh, p); @@ -168,7 +173,7 @@ namespace Functions cell_hint.get() = cell; // Now we can find out about the point - Quadrature quad(qp); + Quadrature quad(qp.get()); FEValues fe_v(mapping, cell->get_fe(), quad, update_hessians); fe_v.reinit(cell); diff --git a/tests/bits/fe_field_function_04_vector.cc b/tests/bits/fe_field_function_04_vector.cc index 3586f30fc6..4e2c42ba5d 100644 --- a/tests/bits/fe_field_function_04_vector.cc +++ b/tests/bits/fe_field_function_04_vector.cc @@ -82,7 +82,7 @@ void test() points.push_back (Point(-0.7+i*0.07,-0.7+j*0.07)); points.push_back (Point(-0.27999999999999992, -0.62999999999999989)); - std::vector > m (points.size()); + std::vector > m (points.size(), Vector(2)); fe_function.vector_value_list (points, m); for (unsigned int i=0; i > m (points.size()); + std::vector > m (points.size(), Vector(2)); fe_function.vector_value_list (points, m); for (unsigned int i=0; i +#include +#include +#include +#include +#include +#include +#include +#include + + +template +class F : public Function +{ + public: + F() : Function(2) {} + virtual void vector_value (const Point &p, + Vector &v) const + { + v = 0; + v[0] = p.square(); + } +}; + + + +template +void test() +{ + const HyperBallBoundary boundary_description; + + Triangulation triangulation; + GridGenerator::hyper_ball (triangulation); + triangulation.set_boundary (0, boundary_description); + triangulation.refine_global (1); + + FESystem fe(FE_Q (2),2); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + // interpolate a quadratic function + // into the space; this function + // can be represented exactly, so + // that we can compare again later + Vector solution(dof_handler.n_dofs()); + VectorTools::interpolate (dof_handler, F(), solution); + + Functions::FEFieldFunction<2> fe_function (dof_handler, solution); + + // add only one points but also set + // the active cell to one that + // doesn't contain the current + // point. the problem happens + // because we walk over a bunch of + // cells in the process of finding + // all of these points and then + // realize when we get to the one + // at the end that the coordinates + // for this point can't be found in + // the cell we have touched last + // (it's too far away from that + // cell, and the inverse mapping + // does not converge + Point point(-0.27999999999999992, -0.62999999999999989); + fe_function.set_active_cell (typename DoFHandler::active_cell_iterator + (&triangulation, + 1, + 4, + &dof_handler)); + + Vector m(2); + fe_function.vector_value (point, m); + + { + Assert (std::fabs(m(0) - point.square()) + < + 1e-10 * std::fabs(m(0) + point.square()), + ExcInternalError()); + + Assert (std::fabs(m(1)) + < + 1e-10, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile("fe_field_function_06_vector/output"); + deallog.attach(logfile); + deallog.depth_console(0); + + test<2>(); + + return 0; +} + diff --git a/tests/bits/fe_field_function_06_vector/cmp/generic b/tests/bits/fe_field_function_06_vector/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/bits/fe_field_function_06_vector/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/bits/fe_field_function_07_vector.cc b/tests/bits/fe_field_function_07_vector.cc new file mode 100644 index 0000000000..55dd72deed --- /dev/null +++ b/tests/bits/fe_field_function_07_vector.cc @@ -0,0 +1,135 @@ +//------------------------------------------------------- +// $Id: fe_field_function_05_vector.cc $ +// Version: $Name$ +// +// Copyright (C) 2005, 2011, 2012 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. +// +//------------------------------------------------------- + + +// FEFieldFunction ran into an assertion after +// Mapping::transform_real_to_unit_cell started throwing exceptions +// when it couldn't find the point on the reference cell that belongs +// to a given point, rather than just silently giving up +// +// this is yet another variant of _04 and _05 and _06 found by +// inspecting the code in FEFieldFunction but the (same) exception is +// triggered in a different place + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +class F : public Function +{ + public: + F() : Function(2) {} + + virtual void vector_value (const Point &p, + Vector &v) const + { + v[0] = p[0]; + v[1] = 0; + } + + virtual void vector_gradient (const Point &p, + std::vector > &v) const + { + v[0] = 0; + v[1] = 0; + v[0][0] = 1; + } +}; + + + +template +void test() +{ + const HyperBallBoundary boundary_description; + + Triangulation triangulation; + GridGenerator::hyper_ball (triangulation); + triangulation.set_boundary (0, boundary_description); + triangulation.refine_global (1); + + FESystem fe(FE_Q (2),2); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + // interpolate a quadratic function + // into the space; this function + // can be represented exactly, so + // that we can compare again later + Vector solution(dof_handler.n_dofs()); + VectorTools::interpolate (dof_handler, F(), solution); + + Functions::FEFieldFunction<2> fe_function (dof_handler, solution); + + // add only one points but also set + // the active cell to one that + // doesn't contain the current + // point. the problem happens + // because we walk over a bunch of + // cells in the process of finding + // all of these points and then + // realize when we get to the one + // at the end that the coordinates + // for this point can't be found in + // the cell we have touched last + // (it's too far away from that + // cell, and the inverse mapping + // does not converge + Point point(-0.27999999999999992, -0.62999999999999989); + fe_function.set_active_cell (typename DoFHandler::active_cell_iterator + (&triangulation, + 1, + 4, + &dof_handler)); + + std::vector > m(2); + fe_function.vector_gradient (point, m); + + Assert (std::fabs(m[0][0] - 1) + < + 1e-10 * std::fabs(m[0][0] + 1), + ExcInternalError()); + Assert (std::fabs(m[0][1]) + < + 1e-10, + ExcInternalError()); + Assert (m[1].norm() + < + 1e-10, + ExcInternalError()); + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile("fe_field_function_07_vector/output"); + deallog.attach(logfile); + deallog.depth_console(0); + + test<2>(); + + return 0; +} + diff --git a/tests/bits/fe_field_function_07_vector/cmp/generic b/tests/bits/fe_field_function_07_vector/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/bits/fe_field_function_07_vector/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/bits/fe_field_function_08_vector.cc b/tests/bits/fe_field_function_08_vector.cc new file mode 100644 index 0000000000..52ca8a09e3 --- /dev/null +++ b/tests/bits/fe_field_function_08_vector.cc @@ -0,0 +1,134 @@ +//------------------------------------------------------- +// $Id: fe_field_function_08_vector.cc $ +// Version: $Name$ +// +// Copyright (C) 2005, 2011, 2012 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. +// +//------------------------------------------------------- + + +// FEFieldFunction ran into an assertion after +// Mapping::transform_real_to_unit_cell started throwing exceptions +// when it couldn't find the point on the reference cell that belongs +// to a given point, rather than just silently giving up +// +// this is yet another variant of _04 through _07 found by inspecting the +// code in FEFieldFunction but the (same) exception is triggered in a +// different place + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +class F : public Function +{ + public: + F() : Function(2) {} + + virtual void vector_value (const Point &p, + Vector &v) const + { + v = 0; + for (unsigned int i=0; i &p, + Vector &v) const + { + v = 0; + v[0] = p.square() * 4 * 3 * dim; + } +}; + + + +template +void test() +{ + const HyperBallBoundary boundary_description; + + Triangulation triangulation; + GridGenerator::hyper_ball (triangulation); + triangulation.set_boundary (0, boundary_description); + triangulation.refine_global (1); + + FESystem fe(FE_Q (4),2); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + // interpolate a quadratic function + // into the space; this function + // can be represented exactly, so + // that we can compare again later + Vector solution(dof_handler.n_dofs()); + VectorTools::interpolate (dof_handler, F(), solution); + + Functions::FEFieldFunction<2> fe_function (dof_handler, solution); + + // add only one points but also set + // the active cell to one that + // doesn't contain the current + // point. the problem happens + // because we walk over a bunch of + // cells in the process of finding + // all of these points and then + // realize when we get to the one + // at the end that the coordinates + // for this point can't be found in + // the cell we have touched last + // (it's too far away from that + // cell, and the inverse mapping + // does not converge + Point point(-0.27999999999999992, -0.62999999999999989); + fe_function.set_active_cell (typename DoFHandler::active_cell_iterator + (&triangulation, + 1, + 4, + &dof_handler)); + + Vector m(2); + fe_function.vector_laplacian (point, m); + + { + Assert (std::fabs(m(0) - point.square()*4*3) + < + 1e-8 * std::fabs(m(0) + point.square()*4*3), + ExcInternalError()); + + Assert (std::fabs(m(1)) + < + 1e-10, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile("fe_field_function_08_vector/output"); + deallog.attach(logfile); + deallog.depth_console(0); + + test<2>(); + + return 0; +} + diff --git a/tests/bits/fe_field_function_08_vector/cmp/generic b/tests/bits/fe_field_function_08_vector/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/bits/fe_field_function_08_vector/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK