]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement a line search strategy to deal with cases where the nonlinear iteration...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 27 Jun 2012 18:44:24 +0000 (18:44 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 27 Jun 2012 18:44:24 +0000 (18:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@25651 0785d39b-7218-0410-832d-ea1e28bc413d

14 files changed:
deal.II/doc/news/changes.h
deal.II/include/deal.II/fe/mapping.h
deal.II/include/deal.II/fe/mapping_q.h
deal.II/include/deal.II/fe/mapping_q1.h
deal.II/source/fe/mapping_q.cc
deal.II/source/fe/mapping_q1.cc
deal.II/source/grid/tria_accessor.cc
tests/fe/mapping_real_to_unit_q4_sphere.cc
tests/fe/mapping_real_to_unit_q4_sphere/cmp/generic
tests/fe/mapping_real_to_unit_q4_sphere_x.cc
tests/fe/mapping_real_to_unit_q4_sphere_x/cmp/generic
tests/fe/mapping_real_to_unit_q4_sphere_xx/cmp/generic [deleted file]
tests/fe/mapping_real_to_unit_q4_sphere_y.cc [moved from tests/fe/mapping_real_to_unit_q4_sphere_xx.cc with 65% similarity]
tests/fe/mapping_real_to_unit_q4_sphere_y/cmp/generic [new file with mode: 0644]

index b998e7f4a82ec439615657364e44edd7442d60fd..b30926173351db32b58bbd0ec9c40efd2e8c105b 100644 (file)
@@ -47,6 +47,22 @@ used for boundary indicators.
 
 
 <ol>
+<li> 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.
+<br>
+(Wolfgang Bangerth, Eric Heien, Sebastian Pauletti, 2012/06/27)
+
 <li> 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
index 53c4e2947c8aa40027729be99d08bcbd0f885b78..ed8ad6077004aec4a71688d8ec3feb12010ca0ca 100644 (file)
@@ -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<dim>
     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:
 
                                      /**
index 75d2122f54806a001780393142133324f0d5f1d7..d5278c09a0c915306bd77c5d7150e34cdbc9b05a 100644 (file)
@@ -116,6 +116,31 @@ class MappingQ : public MappingQ1<dim,spacedim>
                                       * 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<dim>
     transform_real_to_unit_cell (
index ce7d36d57422a5487fe792e0806e3d5fd7ae1adf..5c319e3839d1a39359bdf024ae2aff23d734bb13 100644 (file)
@@ -78,6 +78,31 @@ class MappingQ1 : public Mapping<dim,spacedim>
                                       * 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<dim>
     transform_real_to_unit_cell (
@@ -556,19 +581,6 @@ class MappingQ1 : public Mapping<dim,spacedim>
     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<dim, dim+1> */
     template<int dim_>
index 1d12cbbd1673cac6d97234051c25b14aa241d72f..4105d61f0289dec7783a68f5c44377326e75cc5f 100644 (file)
@@ -1445,19 +1445,47 @@ MappingQ<dim,spacedim>::
 transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                              const Point<spacedim>                            &p) const
 {
-                                   // first a Newton iteration based on a Q1
-                                   // mapping.
-  Point<dim> initial_p_unit =
-    MappingQ1<dim,spacedim>::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<dim> initial_p_unit;
+  bool       initial_p_unit_is_valid = false;
+  try
+    {
+      initial_p_unit
+       = MappingQ1<dim,spacedim>::transform_real_to_unit_cell(cell, p);
+
+      initial_p_unit_is_valid = true;
+    }
+  catch (const typename Mapping<dim,spacedim>::ExcTransformationFailed &)
+    {
+      for (unsigned int d=0; d<dim; ++d)
+       initial_p_unit[d] = 0.5;
+    }
 
                                    // then a Newton iteration based on the
                                    // full MappingQ if we need this. note that
                                    // for interior cells with dim==spacedim,
                                    // the mapping used is in fact a Q1
                                    // mapping, so there is nothing we need to
-                                   // do
-  if (cell->has_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
index cb8f8106cc9c565b28bd4181c9a8f090bbd8264f..4ca013da9a222c80042910573114ddb214433b81 100644 (file)
@@ -31,9 +31,6 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-
-
-
 template <int dim, int spacedim>
 const unsigned int MappingQ1<dim,spacedim>::n_shape_functions;
 
@@ -1678,8 +1675,10 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
                                        // Ignore non vertex support points.
       mdata->mapping_support_points.resize(GeometryInfo<dim>::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<spacedim> 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++<loop_limit)
+  unsigned int newton_iteration=0;
+  while (f.square()>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<const Tensor<1,spacedim>&>(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<dim; ++i)
-        p_unit[i] -= d[i];
-
-                                       // shape values and derivatives
-                                       // at new p_unit point
-      compute_shapes(std::vector<Point<dim> > (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<const Tensor<1,spacedim>&>(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<loop_limit, ExcTransformationFailed());
+      // do a line search
+      double step_length = 1;
+      do
+        {
+          // 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.
+          Point<dim> p_unit_trial = p_unit;
+          for(unsigned int i=0;i<dim; ++i)
+            p_unit_trial[i] -= step_length * delta[i];
+
+          // shape values and derivatives
+          // at new p_unit point
+          compute_shapes(std::vector<Point<dim> > (1, p_unit_trial), mdata);
+
+          // f(x)
+          Point<spacedim> p_real_trial = transform_unit_to_real_cell_internal(mdata);
+          const Point<spacedim> 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<dim,spacedim>::ExcTransformationFailed()));
+
+  // ...the compiler wants us to return something, though we can
+  // of course never get here...
+  return Point<dim>();
 }
 
 
@@ -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<loop_limit, ExcTransformationFailed());
+  AssertThrow (loop<loop_limit, (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
 
   return p_unit;
 }
index db2550087fad713fb90fafaa812e26af91664c64..8f8255323c6d61342a4bdf93196e910a49b7b479 100644 (file)
@@ -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<dim,spacedim> > cell_iterator (*this);
-  return (GeometryInfo<dim>::is_inside_unit_cell (
-              StaticMappingQ1<dim,spacedim>::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<dim,spacedim> > cell_iterator (*this);
+    return (GeometryInfo<dim>::is_inside_unit_cell
+             (StaticMappingQ1<dim,spacedim>::mapping.transform_real_to_unit_cell(cell_iterator, p)));
+  }
+  catch (const Mapping<dim,spacedim>::ExcTransformationFailed &e)
+  {
+    return false;
+  }
 }
 
 
index 16e57afa1de996269bd8ca4426ec594ddcbacdd3..4b53072d42157634f260448eedd0d51aa5ad218a 100644 (file)
@@ -81,7 +81,15 @@ void test_real_to_unit_cell()
   MappingQ<dim> map(4);
   typename Triangulation<dim >::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<dim>::ExcTransformationFailed &)
+    {
+      deallog << "Point is outside!" << std::endl;
+    }
 }
 
 
index 7bea8b1c89d6ac3123605160c7b020fda14fbfd5..d193d0581ee13dac85bc89508803949c51e807c1 100644 (file)
@@ -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!
index 8b018af8b58653ffb60abc31ee65bff5fa0687e0..5fcdbf2c9d51ac8ff3fb161652a06df12ad49b5d 100644 (file)
@@ -74,7 +74,14 @@ void test_real_to_unit_cell()
   MappingQ1<dim> map;
   typename Triangulation<dim >::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<dim>::ExcTransformationFailed &)
+    {
+      deallog << "Point is outside!" << std::endl;
+    }
 }
 
 
index 7bea8b1c89d6ac3123605160c7b020fda14fbfd5..d193d0581ee13dac85bc89508803949c51e807c1 100644 (file)
@@ -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 (file)
index 7bea8b1..0000000
+++ /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
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 483422e8051b75ff1ec787cced775c4b767765c3..69f09bc72d5f728b0fbfe65f332a438d4dd9b3df 100644 (file)
@@ -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 <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/fe/mapping_q.h>
 
-#include <deal.II/numerics/data_out.h>
 
 void test_real_to_unit_cell()
 {
@@ -61,32 +60,44 @@ void test_real_to_unit_cell()
   triangulation.create_triangulation (vertices, cells,
                                      SubCellData());
 
-
-  FE_Q<dim> fe (1);
-  DoFHandler<dim> dh(triangulation);
-  dh.distribute_dofs (fe);
-  Vector<double> dummy;
-  dummy.reinit(dh.n_dofs());
-  DataOut<dim, DoFHandler<dim> > 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<dim> p (-3.56413e+06, 1.74215e+06, 2.14762e+06);
-  // MappingQ1<dim> map;
-  // typename Triangulation<dim >::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<dim> p (-3.56413e+06, 1.74215e+06, 2.14762e+06);
+  MappingQ1<dim> map;
+  typename Triangulation<dim >::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 (file)
index 0000000..37dc17a
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::Caught exception of type N6dealii7MappingILi3ELi3EE23ExcTransformationFailedE

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.