From d3682d333c174d658f0316281e922fed28e71c21 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 27 Jun 2012 18:44:24 +0000 Subject: [PATCH] Implement a line search strategy to deal with cases where the nonlinear iteration in MappingQ1::transform_real_to_unit_cell does not converge. If the line search doesn't converge either, then throw an exception and assume that in fact this means that the point lies outside. Deal with this result in all of the places where we call the function. git-svn-id: https://svn.dealii.org/trunk@25651 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 16 +++ deal.II/include/deal.II/fe/mapping.h | 43 ++++++- deal.II/include/deal.II/fe/mapping_q.h | 25 ++++ deal.II/include/deal.II/fe/mapping_q1.h | 38 +++--- deal.II/source/fe/mapping_q.cc | 42 +++++-- deal.II/source/fe/mapping_q1.cc | 109 +++++++++++------- deal.II/source/grid/tria_accessor.cc | 24 ++-- tests/fe/mapping_real_to_unit_q4_sphere.cc | 10 +- .../cmp/generic | 9 +- tests/fe/mapping_real_to_unit_q4_sphere_x.cc | 9 +- .../cmp/generic | 9 +- .../cmp/generic | 9 -- ...cc => mapping_real_to_unit_q4_sphere_y.cc} | 59 ++++++---- .../cmp/generic | 2 + 14 files changed, 282 insertions(+), 122 deletions(-) delete mode 100644 tests/fe/mapping_real_to_unit_q4_sphere_xx/cmp/generic rename tests/fe/{mapping_real_to_unit_q4_sphere_xx.cc => mapping_real_to_unit_q4_sphere_y.cc} (65%) create mode 100644 tests/fe/mapping_real_to_unit_q4_sphere_y/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index b998e7f4a8..b309261733 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -47,6 +47,22 @@ used for boundary indicators.
    +
  1. Fixed: The MappingQ1::transform_real_to_unit_cell function as +well as the equivalent ones in derived classes sometimes get into +trouble if they are asked to compute the preimage of this point +in reference cell coordinates. This is because for points outside +the reference cell, the mapping from unit to real cell is not +necessarily invertible, and consequently the Newton iteration to +find the preimage did not always converge, leading to an exception. +While this is not entirely wrong (we could, after all, not compute +the desired quantity), not all callers of this function were prepared +to accept this result -- in particular, the function +CellAccessor<3>::point_inside should have really just returned false +in such cases but instead let the exception so generated propagate +through. This should now be fixed. +
    +(Wolfgang Bangerth, Eric Heien, Sebastian Pauletti, 2012/06/27) +
  2. Fixed: The function VectorTools::compute_no_normal_flux_constraints had a bug that led to an exception whenever we were computing constraints for vector fields located on edges shared between two faces of a 3d cell if those diff --git a/deal.II/include/deal.II/fe/mapping.h b/deal.II/include/deal.II/fe/mapping.h index 53c4e2947c..ed8ad60770 100644 --- a/deal.II/include/deal.II/fe/mapping.h +++ b/deal.II/include/deal.II/fe/mapping.h @@ -238,9 +238,9 @@ class Mapping : public Subscriptor /** * Transforms the point @p p on - * the real cell to the point - * @p p_unit on the unit cell - * @p cell and returns @p p_unit. + * the real @p cell to the corresponding + * point on the unit cell, and + * return its coordinates. * * In the codimension one case, * this function returns the @@ -248,6 +248,31 @@ class Mapping : public Subscriptor * point @p p on the curve or * surface identified by the @p * cell. + * + * @note Polynomial mappings from + * the reference (unit) cell coordinates + * to the coordinate system of a real + * cell are not always invertible if + * the point for which the inverse + * mapping is to be computed lies + * outside the cell's boundaries. + * In such cases, the current function + * may fail to compute a point on + * the reference cell whose image + * under the mapping equals the given + * point @p p. If this is the case + * then this function throws an + * exception of type + * Mapping::ExcTransformationFailed . + * Whether the given point @p p lies + * outside the cell can therefore be + * determined by checking whether the + * return reference coordinates lie + * inside of outside the reference + * cell (e.g., using + * GeometryInfo::is_inside_unit_cell) + * or whether the exception mentioned + * above has been thrown. */ virtual Point transform_real_to_unit_cell ( @@ -681,6 +706,18 @@ class Mapping : public Subscriptor */ DeclException0 (ExcInvalidData); + + /** + * Computing the mapping between a + * real space point and a point + * in reference space failed, typically because the given point + * lies outside the cell where the inverse mapping is not + * unique. + * + * @ingroup Exceptions + */ + DeclException0(ExcTransformationFailed); + private: /** diff --git a/deal.II/include/deal.II/fe/mapping_q.h b/deal.II/include/deal.II/fe/mapping_q.h index 75d2122f54..d5278c09a0 100644 --- a/deal.II/include/deal.II/fe/mapping_q.h +++ b/deal.II/include/deal.II/fe/mapping_q.h @@ -116,6 +116,31 @@ class MappingQ : public MappingQ1 * point @p p on the curve or * surface identified by the @p * cell. + * + * @note Polynomial mappings from + * the reference (unit) cell coordinates + * to the coordinate system of a real + * cell are not always invertible if + * the point for which the inverse + * mapping is to be computed lies + * outside the cell's boundaries. + * In such cases, the current function + * may fail to compute a point on + * the reference cell whose image + * under the mapping equals the given + * point @p p. If this is the case + * then this function throws an + * exception of type + * Mapping::ExcTransformationFailed . + * Whether the given point @p p lies + * outside the cell can therefore be + * determined by checking whether the + * return reference coordinates lie + * inside of outside the reference + * cell (e.g., using + * GeometryInfo::is_inside_unit_cell) + * or whether the exception mentioned + * above has been thrown. */ virtual Point transform_real_to_unit_cell ( diff --git a/deal.II/include/deal.II/fe/mapping_q1.h b/deal.II/include/deal.II/fe/mapping_q1.h index ce7d36d574..5c319e3839 100644 --- a/deal.II/include/deal.II/fe/mapping_q1.h +++ b/deal.II/include/deal.II/fe/mapping_q1.h @@ -78,6 +78,31 @@ class MappingQ1 : public Mapping * point @p p on the curve or * surface identified by the @p * cell. + * + * @note Polynomial mappings from + * the reference (unit) cell coordinates + * to the coordinate system of a real + * cell are not always invertible if + * the point for which the inverse + * mapping is to be computed lies + * outside the cell's boundaries. + * In such cases, the current function + * may fail to compute a point on + * the reference cell whose image + * under the mapping equals the given + * point @p p. If this is the case + * then this function throws an + * exception of type + * Mapping::ExcTransformationFailed . + * Whether the given point @p p lies + * outside the cell can therefore be + * determined by checking whether the + * return reference coordinates lie + * inside of outside the reference + * cell (e.g., using + * GeometryInfo::is_inside_unit_cell) + * or whether the exception mentioned + * above has been thrown. */ virtual Point transform_real_to_unit_cell ( @@ -556,19 +581,6 @@ class MappingQ1 : public Mapping virtual bool preserves_vertex_locations () const; - /** - * The Newton iteration computing - * the mapping between a - * real space point and a point - * in reference space did not - * converge, thus resulting in a - * wrong mapping. - * - * @ingroup Exceptions - */ - DeclException0(ExcTransformationFailed); - - protected: /* Trick to templatize transform_real_to_unit_cell */ template diff --git a/deal.II/source/fe/mapping_q.cc b/deal.II/source/fe/mapping_q.cc index 1d12cbbd16..4105d61f02 100644 --- a/deal.II/source/fe/mapping_q.cc +++ b/deal.II/source/fe/mapping_q.cc @@ -1445,19 +1445,47 @@ MappingQ:: transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, const Point &p) const { - // first a Newton iteration based on a Q1 - // mapping. - Point initial_p_unit = - MappingQ1::transform_real_to_unit_cell(cell, p); + // first a Newton iteration based + // on a Q1 mapping to get a good + // starting point, the idea being + // that this is cheaper than trying + // to start with the real mapping + // and likely also more robust. + // + // that said, this doesn't always + // work: there are cases where the + // point is outside the cell and + // the inverse mapping doesn't + // converge. in that case, use the + // center point of the cell as a + // starting point + Point initial_p_unit; + bool initial_p_unit_is_valid = false; + try + { + initial_p_unit + = MappingQ1::transform_real_to_unit_cell(cell, p); + + initial_p_unit_is_valid = true; + } + catch (const typename Mapping::ExcTransformationFailed &) + { + for (unsigned int d=0; dhas_boundary_lines() || - use_mapping_q_on_all_cells || + // do unless the iteration above failed + if ((initial_p_unit_is_valid == false) + || + cell->has_boundary_lines() + || + use_mapping_q_on_all_cells + || (dim!=spacedim) ) { // use the full mapping. in case the diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index cb8f8106cc..4ca013da9a 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -31,9 +31,6 @@ DEAL_II_NAMESPACE_OPEN - - - template const unsigned int MappingQ1::n_shape_functions; @@ -1678,8 +1675,10 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it // Ignore non vertex support points. mdata->mapping_support_points.resize(GeometryInfo::vertices_per_cell); - // perform the Newton iteration and - // return the result + // perform the Newton iteration and + // return the result. note that this + // statement may throw an exception, which + // we simply pass up to the caller return transform_real_to_unit_cell_internal(cell, p, initial_p_unit, *mdata); } @@ -1727,10 +1726,10 @@ transform_real_to_unit_cell_internal Point f = p_real-p; const double eps = 1.e-12*cell->diameter(); - const unsigned int loop_limit = 50; + const unsigned int newton_iteration_limit = 20; - unsigned int loop=0; - while (f.square()>eps*eps && loop++eps*eps) { // f'(x) Tensor<2,spacedim> df; @@ -1745,41 +1744,67 @@ transform_real_to_unit_cell_internal } // Solve [f'(x)]d=f(x) - Tensor<1,spacedim> d; - contract (d, invert(df), static_cast&>(f)); - - // update of p_unit. The - // spacedimth component of - // transformed point is simply - // ignored in codimension one - // case. When this component is - // not zero, then we are - // projecting the point to the - // surface or curve identified - // by the cell. - for(unsigned int i=0;i > (1, p_unit), mdata); - - // f(x) - p_real = transform_unit_to_real_cell_internal(mdata); - f = p_real-p; - } + Tensor<1,spacedim> delta; + contract (delta, invert(df), static_cast&>(f)); - // Here we check that in the last - // execution of while the first - // condition was already wrong, - // meaning the residual was below - // eps. Only if the first condition - // failed, loop will have been - // increased and tested, and thus - // havereached the limit. - AssertThrow (loop p_unit_trial = p_unit; + for(unsigned int i=0;i > (1, p_unit_trial), mdata); + + // f(x) + Point p_real_trial = transform_unit_to_real_cell_internal(mdata); + const Point f_trial = p_real_trial-p; + + // see if we are making progress with the current step length + // and if not, reduce it by a factor of two and try again + if (f_trial.norm() < f.norm()) + { + p_real = p_real_trial; + p_unit = p_unit_trial; + f = f_trial; + break; + } + else if (step_length > 0.05) + step_length /= 2; + else + goto failure; + } + while (true); + + ++newton_iteration; + if (newton_iteration > newton_iteration_limit) + goto failure; + } return p_unit; + + // if we get to the following label, then we have either run out + // of Newton iterations, or the line search has not converged. + // in either case, we need to give up, so throw an exception that + // can then be caught + failure: + AssertThrow (false, (typename Mapping::ExcTransformationFailed())); + + // ...the compiler wants us to return something, though we can + // of course never get here... + return Point(); } @@ -1944,6 +1969,8 @@ transform_real_to_unit_cell_internal_codim1 } } +//TODO: implement a line search here in much the same way as for +// the corresponding function above that does so for codim==0. p_minus_F = p; p_minus_F -= transform_unit_to_real_cell_internal(mdata); @@ -1965,7 +1992,7 @@ transform_real_to_unit_cell_internal_codim1 // failed, loop will have been // increased and tested, and thus // have reached the limit. - AssertThrow (loop::ExcTransformationFailed())); return p_unit; } diff --git a/deal.II/source/grid/tria_accessor.cc b/deal.II/source/grid/tria_accessor.cc index db2550087f..8f8255323c 100644 --- a/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/source/grid/tria_accessor.cc @@ -1184,13 +1184,23 @@ bool CellAccessor<3>::point_inside (const Point<3> &p) const if ((p[d] < minp[d]) || (p[d] > maxp[d])) return false; - // now we need to check more - // carefully: transform to the - // unit cube - // and check there. - const TriaRawIterator< CellAccessor > cell_iterator (*this); - return (GeometryInfo::is_inside_unit_cell ( - StaticMappingQ1::mapping.transform_real_to_unit_cell(cell_iterator, p))); + // now we need to check more carefully: transform to the + // unit cube and check there. unfortunately, this isn't + // completely trivial since the transform_real_to_unit_cell + // function may throw an exception that indicates that the + // point given could not be inverted. we take this as a sign + // that the point actually lies outside, as also documented + // for that function + try + { + const TriaRawIterator< CellAccessor > cell_iterator (*this); + return (GeometryInfo::is_inside_unit_cell + (StaticMappingQ1::mapping.transform_real_to_unit_cell(cell_iterator, p))); + } + catch (const Mapping::ExcTransformationFailed &e) + { + return false; + } } diff --git a/tests/fe/mapping_real_to_unit_q4_sphere.cc b/tests/fe/mapping_real_to_unit_q4_sphere.cc index 16e57afa1d..4b53072d42 100644 --- a/tests/fe/mapping_real_to_unit_q4_sphere.cc +++ b/tests/fe/mapping_real_to_unit_q4_sphere.cc @@ -81,7 +81,15 @@ void test_real_to_unit_cell() MappingQ map(4); typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); - deallog << map.transform_real_to_unit_cell(cell,p) << std::endl; + + try + { + map.transform_real_to_unit_cell(cell,p); + } + catch (const typename Mapping::ExcTransformationFailed &) + { + deallog << "Point is outside!" << std::endl; + } } diff --git a/tests/fe/mapping_real_to_unit_q4_sphere/cmp/generic b/tests/fe/mapping_real_to_unit_q4_sphere/cmp/generic index 7bea8b1c89..d193d0581e 100644 --- a/tests/fe/mapping_real_to_unit_q4_sphere/cmp/generic +++ b/tests/fe/mapping_real_to_unit_q4_sphere/cmp/generic @@ -1,9 +1,2 @@ -DEAL::dim=1, spacedim=1 -DEAL::OK -DEAL::dim=2, spacedim=2 -DEAL::OK -DEAL::dim=3, spacedim=3 -DEAL::OK -DEAL::dim=1, spacedim=2 -DEAL::OK +DEAL::Point is outside! diff --git a/tests/fe/mapping_real_to_unit_q4_sphere_x.cc b/tests/fe/mapping_real_to_unit_q4_sphere_x.cc index 8b018af8b5..5fcdbf2c9d 100644 --- a/tests/fe/mapping_real_to_unit_q4_sphere_x.cc +++ b/tests/fe/mapping_real_to_unit_q4_sphere_x.cc @@ -74,7 +74,14 @@ void test_real_to_unit_cell() MappingQ1 map; typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); - deallog << map.transform_real_to_unit_cell(cell,p) << std::endl; + try + { + map.transform_real_to_unit_cell(cell,p); + } + catch (const typename Mapping::ExcTransformationFailed &) + { + deallog << "Point is outside!" << std::endl; + } } diff --git a/tests/fe/mapping_real_to_unit_q4_sphere_x/cmp/generic b/tests/fe/mapping_real_to_unit_q4_sphere_x/cmp/generic index 7bea8b1c89..d193d0581e 100644 --- a/tests/fe/mapping_real_to_unit_q4_sphere_x/cmp/generic +++ b/tests/fe/mapping_real_to_unit_q4_sphere_x/cmp/generic @@ -1,9 +1,2 @@ -DEAL::dim=1, spacedim=1 -DEAL::OK -DEAL::dim=2, spacedim=2 -DEAL::OK -DEAL::dim=3, spacedim=3 -DEAL::OK -DEAL::dim=1, spacedim=2 -DEAL::OK +DEAL::Point is outside! diff --git a/tests/fe/mapping_real_to_unit_q4_sphere_xx/cmp/generic b/tests/fe/mapping_real_to_unit_q4_sphere_xx/cmp/generic deleted file mode 100644 index 7bea8b1c89..0000000000 --- a/tests/fe/mapping_real_to_unit_q4_sphere_xx/cmp/generic +++ /dev/null @@ -1,9 +0,0 @@ - -DEAL::dim=1, spacedim=1 -DEAL::OK -DEAL::dim=2, spacedim=2 -DEAL::OK -DEAL::dim=3, spacedim=3 -DEAL::OK -DEAL::dim=1, spacedim=2 -DEAL::OK diff --git a/tests/fe/mapping_real_to_unit_q4_sphere_xx.cc b/tests/fe/mapping_real_to_unit_q4_sphere_y.cc similarity index 65% rename from tests/fe/mapping_real_to_unit_q4_sphere_xx.cc rename to tests/fe/mapping_real_to_unit_q4_sphere_y.cc index 483422e805..69f09bc72d 100644 --- a/tests/fe/mapping_real_to_unit_q4_sphere_xx.cc +++ b/tests/fe/mapping_real_to_unit_q4_sphere_y.cc @@ -12,9 +12,9 @@ //----------------------------------------------------------------------------- -// a redux of mapping_real_to_unit_q4_sphere_x. -// It seems that the point is outside the cell?? - +// a redux of mapping_real_to_unit_q4_sphere. one doesn't in fact need +// the Q4 mapping, the problem already happens in the initial guess +// generation using a Q1 mapping #include "../tests.h" @@ -24,7 +24,6 @@ #include #include -#include void test_real_to_unit_cell() { @@ -61,32 +60,44 @@ void test_real_to_unit_cell() triangulation.create_triangulation (vertices, cells, SubCellData()); - - FE_Q fe (1); - DoFHandler dh(triangulation); - dh.distribute_dofs (fe); - Vector dummy; - dummy.reinit(dh.n_dofs()); - DataOut > data_out; - data_out.attach_dof_handler (dh); - data_out.add_data_vector (dummy,"dummy"); - data_out.build_patches (); - std::ofstream output ("mapping_real_to_unit_q4_sphere_xx/plot.vtk"); - data_out.write_vtk (output); - - - // const Point p (-3.56413e+06, 1.74215e+06, 2.14762e+06); - // MappingQ1 map; - // typename Triangulation::active_cell_iterator - // cell = triangulation.begin_active(); - // deallog << map.transform_real_to_unit_cell(cell,p) << std::endl; + // set the boundary indicator for + // one face and adjacent edges of + // the single cell + triangulation.set_boundary (1, boundary); + triangulation.begin_active()->face(5)->set_all_boundary_indicators (1); + + // now try to find the coordinates + // of the following point in the + // reference coordinate system of + // the cell + const Point p (-3.56413e+06, 1.74215e+06, 2.14762e+06); + MappingQ1 map; + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + + // the following call will fail + // because the point is outside the + // bounds of the cell. catch the + // exception and make sure that + // there is output that documents + // that the function threw an + // exception + try + { + map.transform_real_to_unit_cell(cell,p); + } + catch (const std::exception &e) + { + deallog << "Caught exception of type " << typeid(e).name() + << std::endl; + } } int main() { - std::ofstream logfile ("mapping_real_to_unit_q4_sphere_xx/output"); + std::ofstream logfile ("mapping_real_to_unit_q4_sphere_y/output"); deallog.attach(logfile); deallog.depth_console(0); deallog.threshold_double(1.e-10); diff --git a/tests/fe/mapping_real_to_unit_q4_sphere_y/cmp/generic b/tests/fe/mapping_real_to_unit_q4_sphere_y/cmp/generic new file mode 100644 index 0000000000..37dc17a6cc --- /dev/null +++ b/tests/fe/mapping_real_to_unit_q4_sphere_y/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::Caught exception of type N6dealii7MappingILi3ELi3EE23ExcTransformationFailedE -- 2.39.5