From f298ea2de0343cdf29468ad3cc766774f82025c8 Mon Sep 17 00:00:00 2001
From: Marc Fehling <mafehling.git@gmail.com>
Date: Sat, 21 Jan 2023 15:54:16 -0700
Subject: [PATCH] Minor changes to PoisseuilleFlow.

---
 include/deal.II/base/flow_function.h |  2 +-
 source/base/flow_function.cc         | 25 +++++++++++--------------
 2 files changed, 12 insertions(+), 15 deletions(-)

diff --git a/include/deal.II/base/flow_function.h b/include/deal.II/base/flow_function.h
index 33b213ca6b..6de3c60247 100644
--- a/include/deal.II/base/flow_function.h
+++ b/include/deal.II/base/flow_function.h
@@ -174,7 +174,7 @@ namespace Functions
                       std::vector<std::vector<double>> &values) const override;
 
   private:
-    const double radius;
+    const double inv_sqr_radius;
     const double Reynolds;
   };
 
diff --git a/source/base/flow_function.cc b/source/base/flow_function.cc
index c70b48b86a..08674fa824 100644
--- a/source/base/flow_function.cc
+++ b/source/base/flow_function.cc
@@ -189,7 +189,7 @@ namespace Functions
 
   template <int dim>
   PoisseuilleFlow<dim>::PoisseuilleFlow(const double r, const double Re)
-    : radius(r)
+    : inv_sqr_radius(1 / r / r)
     , Reynolds(Re)
   {
     Assert(Reynolds != 0., ExcMessage("Reynolds number cannot be zero"));
@@ -203,8 +203,7 @@ namespace Functions
     const std::vector<Point<dim>> &   points,
     std::vector<std::vector<double>> &values) const
   {
-    unsigned int n       = points.size();
-    double       stretch = 1. / radius;
+    const unsigned int n = points.size();
 
     Assert(values.size() == dim + 1,
            ExcDimensionMismatch(values.size(), dim + 1));
@@ -214,13 +213,12 @@ namespace Functions
     for (unsigned int k = 0; k < n; ++k)
       {
         const Point<dim> &p = points[k];
-        // First, compute the
-        // square of the distance to
-        // the x-axis divided by the
-        // radius.
+        // First, compute the square of the distance to the x-axis divided by
+        // the radius.
         double r2 = 0;
         for (unsigned int d = 1; d < dim; ++d)
-          r2 += p(d) * p(d) * stretch * stretch;
+          r2 += p(d) * p(d);
+        r2 *= inv_sqr_radius;
 
         // x-velocity
         values[0][k] = 1. - r2;
@@ -228,7 +226,7 @@ namespace Functions
         for (unsigned int d = 1; d < dim; ++d)
           values[d][k] = 0.;
         // pressure
-        values[dim][k] = -2 * (dim - 1) * stretch * stretch * p(0) / Reynolds +
+        values[dim][k] = -2 * (dim - 1) * inv_sqr_radius * p(0) / Reynolds +
                          this->mean_pressure;
       }
   }
@@ -241,8 +239,7 @@ namespace Functions
     const std::vector<Point<dim>> &           points,
     std::vector<std::vector<Tensor<1, dim>>> &values) const
   {
-    unsigned int n       = points.size();
-    double       stretch = 1. / radius;
+    const unsigned int n = points.size();
 
     Assert(values.size() == dim + 1,
            ExcDimensionMismatch(values.size(), dim + 1));
@@ -255,12 +252,12 @@ namespace Functions
         // x-velocity
         values[0][k][0] = 0.;
         for (unsigned int d = 1; d < dim; ++d)
-          values[0][k][d] = -2. * p(d) * stretch * stretch;
+          values[0][k][d] = -2. * p(d) * inv_sqr_radius;
         // other velocities
         for (unsigned int d = 1; d < dim; ++d)
           values[d][k] = 0.;
         // pressure
-        values[dim][k][0] = -2 * (dim - 1) * stretch * stretch / Reynolds;
+        values[dim][k][0] = -2 * (dim - 1) * inv_sqr_radius / Reynolds;
         for (unsigned int d = 1; d < dim; ++d)
           values[dim][k][d] = 0.;
       }
@@ -274,7 +271,7 @@ namespace Functions
     const std::vector<Point<dim>> &   points,
     std::vector<std::vector<double>> &values) const
   {
-    unsigned int n = points.size();
+    const unsigned int n = points.size();
     (void)n;
     Assert(values.size() == dim + 1,
            ExcDimensionMismatch(values.size(), dim + 1));
-- 
2.39.5