From f298ea2de0343cdf29468ad3cc766774f82025c8 Mon Sep 17 00:00:00 2001 From: Marc Fehling 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> &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 PoisseuilleFlow::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> & points, std::vector> &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 &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> & points, std::vector>> &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> & points, std::vector> &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