From ef531126c6f72b7a7e64ff00f58d7662ccaa9a69 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Tue, 17 Dec 2024 21:18:32 -0700
Subject: [PATCH] Avoid the use of infinities.

---
 include/deal.II/fe/mapping.h            | 3 ++-
 include/deal.II/fe/mapping_q_internal.h | 2 +-
 source/fe/mapping.cc                    | 4 +++-
 source/fe/mapping_q.cc                  | 4 ++--
 4 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h
index 3978d43e61..e578c85354 100644
--- a/include/deal.II/fe/mapping.h
+++ b/include/deal.II/fe/mapping.h
@@ -498,7 +498,8 @@ public:
    * MappingQ. The only difference in behavior is that this function
    * will never throw an ExcTransformationFailed() exception. If the
    * transformation fails for `real_points[i]`, the returned `unit_points[i]`
-   * contains std::numeric_limits<double>::infinity() as the first entry.
+   * contains std::numeric_limits<double>::lowest() as the first component
+   * of the point, marking this one point as invalid.
    */
   virtual void
   transform_points_real_to_unit_cell(
diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h
index e4299356bc..210f90565f 100644
--- a/include/deal.II/fe/mapping_q_internal.h
+++ b/include/deal.II/fe/mapping_q_internal.h
@@ -552,7 +552,7 @@ namespace internal
       const unsigned int newton_iteration_limit = 20;
 
       Point<dim, Number> invalid_point;
-      invalid_point[0]                = std::numeric_limits<double>::infinity();
+      invalid_point[0]                = std::numeric_limits<double>::lowest();
       bool tried_project_to_unit_cell = false;
 
       unsigned int newton_iteration            = 0;
diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc
index 314a2d18a8..c7fb0891e2 100644
--- a/source/fe/mapping.cc
+++ b/source/fe/mapping.cc
@@ -142,8 +142,10 @@ Mapping<dim, spacedim>::transform_points_real_to_unit_cell(
         }
       catch (typename Mapping<dim>::ExcTransformationFailed &)
         {
+          // If the transformation for this one point failed, mark it
+          // as invalid as described in the documentation.
           unit_points[i]    = Point<dim>();
-          unit_points[i][0] = std::numeric_limits<double>::infinity();
+          unit_points[i][0] = std::numeric_limits<double>::lowest();
         }
     }
 }
diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc
index 60bdd76cce..734fb98646 100644
--- a/source/fe/mapping_q.cc
+++ b/source/fe/mapping_q.cc
@@ -615,7 +615,7 @@ MappingQ<dim, spacedim>::transform_real_to_unit_cell(
   // statement may throw an exception, which we simply pass up to the caller
   const Point<dim> p_unit =
     this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
-  AssertThrow(numbers::is_finite(p_unit[0]),
+  AssertThrow(p_unit[0] != std::numeric_limits<double>::lowest(),
               (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
   return p_unit;
 }
@@ -690,7 +690,7 @@ MappingQ<dim, spacedim>::transform_points_real_to_unit_cell(
         // determinants) from other SIMD lanes. Repeat the computation in this
         // unlikely case with scalar arguments.
         for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-          if (numbers::is_finite(unit_point[0][j]))
+          if (unit_point[0][j] != std::numeric_limits<double>::lowest())
             for (unsigned int d = 0; d < dim; ++d)
               unit_points[i + j][d] = unit_point[d][j];
           else
-- 
2.39.5