Transfinite interpolation: Fix refresh of Jacobian

diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index a2110f226b79d49eef8d2d8f8155768a86d2b40b..d1f0813f342a73420a4fabf6cae526e9ba60586e 100644 (file)
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -2154,13 +2154,14 @@ TransfiniteInterpolationManifold<dim, spacedim>::pull_back(
       const Tensor<1, spacedim> old_residual = residual;
       while (alpha > 1e-4)
         {
-          Point<dim> guess = chart_point + alpha * update;
-          residual =
+          Point<dim>                guess = chart_point + alpha * update;
+          const Tensor<1, spacedim> residual_guess =
             point - compute_transfinite_interpolation(
                       *cell, guess, coarse_cell_is_flat[cell->index()]);
-          const double residual_norm_new = residual.norm_square();
+          const double residual_norm_new = residual_guess.norm_square();
           if (residual_norm_new < residual_norm_square)
             {
+              residual             = residual_guess;
               residual_norm_square = residual_norm_new;
               chart_point += alpha * update;
               break;
@@ -2198,7 +2199,7 @@ TransfiniteInterpolationManifold<dim, spacedim>::pull_back(
       // prevent division by zero. This number should be scale-invariant
       // because Jinv_deltaf carries no units and x is in reference
       // coordinates.
-      if (std::abs(delta_x * Jinv_deltaf) > 1e-12)
+      if (std::abs(delta_x * Jinv_deltaf) > 1e-12 && !must_recompute_jacobian)
         {
           const Tensor<1, dim> factor =
             (delta_x - Jinv_deltaf) / (delta_x * Jinv_deltaf);