]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix the failure of find_cell_10a by implementing a better stopping criterion for...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 18 Nov 2013 01:31:26 +0000 (01:31 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 18 Nov 2013 01:31:26 +0000 (01:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@31693 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/source/fe/mapping_q1.cc
tests/bits/find_cell_10a.debug.output [deleted file]
tests/bits/find_cell_10a.output [new file with mode: 0644]

index 04671e86902c9ca96f0e34d168c4ff9ab5a3b329..0c189bee36e2f7f4e546d7f89141595ae35ed961 100644 (file)
@@ -94,15 +94,6 @@ inconvenience this causes.
 
 
 <ol>
-  <li>
-  Fixed: VectorTools::compute_no_normal_flux_constraints had a bug that
-  only appeared in rare cases at vertices of the domain if one adjacent
-  cell had two boundary indicators selected for no normal flux and another
-  had only one. This is now fixed.
-  <br>
-  (Wolfgang Bangerth, 2013/11/17)
-  </li>
-
   <li> New: introduced "make test" that runs a minimal set of tests. We
   encourage every user to run this, especially if they run in to problems.
   The tests are automatically picked depending on the configuration and
@@ -220,6 +211,24 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>
+  Fixed: MappingQ1::transform_real_to_unit_cell() could fail in some
+  cases with very elongated and twisted cells. This should now be fixed
+  with an algorithm that uses a better method of computing the Newton
+  convergence.
+  <br>
+  (Wolfgang Bangerth, 2013/11/17)
+  </li>
+
+  <li>
+  Fixed: VectorTools::compute_no_normal_flux_constraints had a bug that
+  only appeared in rare cases at vertices of the domain if one adjacent
+  cell had two boundary indicators selected for no normal flux and another
+  had only one. This is now fixed.
+  <br>
+  (Wolfgang Bangerth, 2013/11/17)
+  </li>
+
   <li> Fixed: dealii::FETools::interpolation_difference was
   not working for TrilinosWrappers::MPI::Vectors with ghost
   entries. The TrilinosWrappers::VectorBase class has now a
index 530a199602515c3c29e8cea08d7f33bbedacfc66..7c708daea1560159ed56a1d9de9432393072a969 100644 (file)
@@ -1684,33 +1684,69 @@ transform_real_to_unit_cell_internal
 
 
   // Newton iteration to solve
-  // f(x)=p(x)-p=0
-  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
+  //    f(x)=p(x)-p=0
+  // where we are looking for 'x' and p(x) is the forward transformation
+  // from unit to real cell. We solve this using a Newton iteration
+  //    x_{n+1}=x_n-[f'(x)]^{-1}f(x)
+  // The start value is set to be the linear approximation to the cell
 
-  // The start value was set to be the
-  // linear approximation to the cell
-
-  // The shape values and derivatives
-  // of the mapping at this point are
+  // The shape values and derivatives of the mapping at this point are
   // previously computed.
 
-  // For the <2,3>  case there is a
-  // template specialization.
-
-  // f(x)
-
   Point<dim> p_unit = initial_p_unit;
 
   compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
-  Point<spacedim> p_real(transform_unit_to_real_cell_internal(mdata));
+  Point<spacedim> p_real = transform_unit_to_real_cell_internal(mdata);
   Point<spacedim> f = p_real-p;
 
-  const double eps = 1.e-12*cell->diameter();
+  // early out if we already have our point
+  if (f.square() < 1e-24 * cell->diameter() * cell->diameter())
+    return p_unit;
+  
+  // we need to compare the position of the computed p(x) against the given
+  // point 'p'. We will terminate the iteration and return 'x' if they are
+  // less than eps apart. The question is how to choose eps -- or, put maybe
+  // more generally: in which norm we want these 'p' and 'p(x)' to be eps
+  // apart.
+  //
+  // the question is difficult since we may have to deal with very elongated
+  // cells where we may achieve 1e-12*h for the distance of these two points
+  // in the 'long' direction, but achieving this tolerance in the 'short'
+  // direction of the cell may not be possible
+  //
+  // what we do instead is then to terminate iterations if
+  //    \| p(x) - p \|_A < eps
+  // where the A-norm is somehow induced by the transformation of the cell.
+  // in particular, we want to measure distances relative to the sizes of
+  // the cell in its principal directions.
+  //
+  // to define what exactly A should be, note that to first order we have
+  // the following (assuming that x* is the solution of the problem, i.e.,
+  // p(x*)=p):
+  //    p(x) - p = p(x) - p(x*)
+  //             = -grad p(x) * (x*-x) + higher order terms
+  // This suggest to measure with a norm that corresponds to
+  //    A = {[grad p(x]^T [grad p(x)]}^{-1}
+  // because then
+  //    \| p(x) - p \|_A  \approx  \| x - x* \|
+  // Consequently, we will try to enforce that
+  //    \| p(x) - p \|_A  =  \| f \|  <=  eps
+  //
+  // Note that using this norm is a bit dangerous since the norm changes
+  // in every iteration (A isn't fixed by depends on xk). However, if the
+  // cell is not too deformed (it may be stretched, but not twisted) then
+  // the mapping is almost linear and A is indeed constant or nearly so.
+  const double eps = 1.e-11;
   const unsigned int newton_iteration_limit = 20;
 
   unsigned int newton_iteration=0;
-  while (f.square()>eps*eps)
+  double last_f_weighted_norm = -1./0.;
+  do
     {
+#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
+      std::cout << "Newton iteration " << newton_iteration << std::endl;
+#endif
+      
       // f'(x)
       Tensor<2,spacedim> df;
       for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
@@ -1725,21 +1761,21 @@ transform_real_to_unit_cell_internal
 
       // Solve  [f'(x)]d=f(x)
       Tensor<1,spacedim> delta;
-      contract (delta, invert(df), static_cast<const Tensor<1,spacedim>&>(f));
+      Tensor<2,spacedim> df_inverse = invert(df);
+      contract (delta, df_inverse, static_cast<const Tensor<1,spacedim>&>(f));
 
+#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
+      std::cout << "   delta=" << delta  << std::endl;
+#endif
+      
       // 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.
+          // update of p_unit. The spacedim-th 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];
@@ -1752,8 +1788,20 @@ transform_real_to_unit_cell_internal
           Point<spacedim> p_real_trial = transform_unit_to_real_cell_internal(mdata);
           const Point<spacedim> f_trial = p_real_trial-p;
 
+#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
+         std::cout << "     step_length=" << step_length << std::endl
+                   << "       ||f ||   =" << f.norm() << std::endl
+                   << "       ||f*||   =" << f_trial.norm() << std::endl
+                   << "       ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl;
+#endif
+         
           // see if we are making progress with the current step length
           // and if not, reduce it by a factor of two and try again
+         //
+         // strictly speaking, we should probably use the same norm as we use
+         // for the outer algorithm. in practice, line search is just a
+         // crutch to find a "reasonable" step length, and so using the l2
+         // norm is probably just fine
           if (f_trial.norm() < f.norm())
             {
               p_real = p_real_trial;
@@ -1764,27 +1812,20 @@ transform_real_to_unit_cell_internal
           else if (step_length > 0.05)
             step_length /= 2;
           else
-            goto failure;
+           AssertThrow (false,
+                        (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
         }
       while (true);
 
       ++newton_iteration;
       if (newton_iteration > newton_iteration_limit)
-        goto failure;
+       AssertThrow (false,
+                    (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
+      last_f_weighted_norm = (df_inverse * f).norm();
     }
-
+  while (last_f_weighted_norm > eps);
+    
   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>();
 }
 
 
diff --git a/tests/bits/find_cell_10a.debug.output b/tests/bits/find_cell_10a.debug.output
deleted file mode 100644 (file)
index d6d313c..0000000
+++ /dev/null
@@ -1,16 +0,0 @@
-
-DEAL::cells 31242
-DEAL::1:
-DEAL::cell 27707
-DEAL::vertex 20297
-DEAL::0.0969822 1125.56
-DEAL::vertex 30988
-DEAL::0.0612193 1125.61
-DEAL::vertex 1298
-DEAL::0.00000 1125.52
-DEAL::vertex 1297
-DEAL::0.00000 1125.60
-DEAL::point 0.706650 0.0885132
-DEAL::1b:
-DEAL::2:
-DEAL::done
diff --git a/tests/bits/find_cell_10a.output b/tests/bits/find_cell_10a.output
new file mode 100644 (file)
index 0000000..9100abf
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::A: 7.0665049837221305e-01 8.8513170467503358e-02
+DEAL::B: 7.0665049836814942e-01 8.8513170469451799e-02
+DEAL::done

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.