From aa8997d23048671136c099ba87d5d39c44b4a470 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 18 Nov 2013 01:31:26 +0000 Subject: [PATCH] Fix the failure of find_cell_10a by implementing a better stopping criterion for the Newton iteration, using a variable metric Newton method. git-svn-id: https://svn.dealii.org/trunk@31693 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 27 ++++-- deal.II/source/fe/mapping_q1.cc | 119 +++++++++++++++++--------- tests/bits/find_cell_10a.debug.output | 16 ---- tests/bits/find_cell_10a.output | 4 + 4 files changed, 102 insertions(+), 64 deletions(-) delete mode 100644 tests/bits/find_cell_10a.debug.output create mode 100644 tests/bits/find_cell_10a.output diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 04671e8690..0c189bee36 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -94,15 +94,6 @@ inconvenience this causes.
    -
  1. - 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. -
    - (Wolfgang Bangerth, 2013/11/17) -
    2. -
  2. 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.

    Specific improvements

      +
    1. + 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. +
      + (Wolfgang Bangerth, 2013/11/17) +
      2. + +
    2. + 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. +
      + (Wolfgang Bangerth, 2013/11/17) +
      3. +
    3. Fixed: dealii::FETools::interpolation_difference was not working for TrilinosWrappers::MPI::Vectors with ghost entries. The TrilinosWrappers::VectorBase class has now a diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index 530a199602..7c708daea1 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -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 p_unit = initial_p_unit; compute_shapes(std::vector > (1, p_unit), mdata); - Point p_real(transform_unit_to_real_cell_internal(mdata)); + Point p_real = transform_unit_to_real_cell_internal(mdata); Point 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 delta; - contract (delta, invert(df), static_cast&>(f)); + Tensor<2,spacedim> df_inverse = invert(df); + contract (delta, df_inverse, static_cast&>(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 p_unit_trial = p_unit; for (unsigned int i=0; i p_real_trial = transform_unit_to_real_cell_internal(mdata); const Point 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::ExcTransformationFailed())); } while (true); ++newton_iteration; if (newton_iteration > newton_iteration_limit) - goto failure; + AssertThrow (false, + (typename Mapping::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::ExcTransformationFailed())); - - // ...the compiler wants us to return something, though we can - // of course never get here... - return Point(); } diff --git a/tests/bits/find_cell_10a.debug.output b/tests/bits/find_cell_10a.debug.output deleted file mode 100644 index d6d313c0bd..0000000000 --- a/tests/bits/find_cell_10a.debug.output +++ /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 index 0000000000..9100abfcf3 --- /dev/null +++ b/tests/bits/find_cell_10a.output @@ -0,0 +1,4 @@ + +DEAL::A: 7.0665049837221305e-01 8.8513170467503358e-02 +DEAL::B: 7.0665049836814942e-01 8.8513170469451799e-02 +DEAL::done -- 2.39.5