From 38aa8c2d2f9005a6eba838d6dbfcfff317263c26 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Tue, 16 Jul 2024 21:44:32 +0200
Subject: [PATCH] Use std::numbers::pi replacing boost constants

Remove unused boost/... headers.
---
 .../Utils/GeoTools/addDataToRaster.cpp        |  8 ++---
 GeoLib/Polygon.cpp                            |  1 -
 MaterialLib/MPL/Properties/Exponential.cpp    |  1 -
 MeshToolsLib/MeshQuality/AngleSkewMetric.cpp  | 34 +++++++++----------
 MeshToolsLib/MeshSurfaceExtraction.cpp        |  5 ++-
 .../Fem/FiniteElement/TemplateIsoparametric.h |  4 +--
 ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp    |  4 +--
 ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp    |  4 +--
 ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp    |  4 +--
 .../HeatTransportBHE/BHE/BoreholeGeometry.h   |  8 ++---
 ProcessLib/HeatTransportBHE/BHE/Pipe.h        | 10 +++---
 .../BHE/ThermalResistancesCoaxial.h           | 13 +++----
 ProcessLib/PhaseField/PhaseFieldFEM-impl.h    |  6 ++--
 .../GeoLib/TestEquidistantPointGeneration.cpp |  4 +--
 Tests/GeoLib/TestSurfaceIsPointInSurface.cpp  |  5 ++-
 .../TestGasPressureDependentPermeability.cpp  |  1 -
 .../TestStrainDependentPermeability.cpp       |  1 -
 Tests/MeshLib/TestTetQualityCriteria.cpp      |  7 ++--
 Tests/MeshLib/TestTriQualityCriteria.cpp      | 10 +++---
 Tests/NumLib/ODEs.h                           |  8 ++---
 20 files changed, 62 insertions(+), 76 deletions(-)

diff --git a/Applications/Utils/GeoTools/addDataToRaster.cpp b/Applications/Utils/GeoTools/addDataToRaster.cpp
index 607403bc972..425bfb5d689 100644
--- a/Applications/Utils/GeoTools/addDataToRaster.cpp
+++ b/Applications/Utils/GeoTools/addDataToRaster.cpp
@@ -16,9 +16,9 @@
 #endif
 
 #include <algorithm>
-#include <boost/math/constants/constants.hpp>
 #include <cmath>
 #include <memory>
+#include <numbers>
 #include <numeric>
 
 #include "GeoLib/AABB.h"
@@ -45,10 +45,8 @@ double computeSinXSinY(GeoLib::Point const& point, GeoLib::AABB const& aabb)
     auto const aabb_size = aabb.getMaxPoint() - aabb.getMinPoint();
     auto const offset = aabb.getMinPoint();
 
-    return std::sin((point[0] - offset[0]) / aabb_size[0] *
-                    boost::math::double_constants::pi) *
-           std::sin((point[1] - offset[1]) / aabb_size[1] *
-                    boost::math::double_constants::pi);
+    return std::sin((point[0] - offset[0]) / aabb_size[0] * std::numbers::pi) *
+           std::sin((point[1] - offset[1]) / aabb_size[1] * std::numbers::pi);
 }
 
 double computeStepFunction(GeoLib::Point const& point, GeoLib::AABB const& aabb)
diff --git a/GeoLib/Polygon.cpp b/GeoLib/Polygon.cpp
index ab2616f1ab7..b2a686412c2 100644
--- a/GeoLib/Polygon.cpp
+++ b/GeoLib/Polygon.cpp
@@ -15,7 +15,6 @@
 #include "Polygon.h"
 
 #include <algorithm>
-#include <boost/math/constants/constants.hpp>
 
 #include "AnalyticalGeometry.h"
 #include "BaseLib/quicksort.h"
diff --git a/MaterialLib/MPL/Properties/Exponential.cpp b/MaterialLib/MPL/Properties/Exponential.cpp
index b5490c270c2..5e2f6a10a02 100644
--- a/MaterialLib/MPL/Properties/Exponential.cpp
+++ b/MaterialLib/MPL/Properties/Exponential.cpp
@@ -10,7 +10,6 @@
 
 #include "MaterialLib/MPL/Properties/Exponential.h"
 
-#include <boost/math/special_functions/pow.hpp>
 #include <cmath>
 
 namespace MaterialPropertyLib
diff --git a/MeshToolsLib/MeshQuality/AngleSkewMetric.cpp b/MeshToolsLib/MeshQuality/AngleSkewMetric.cpp
index 29ae02268ee..a639aa9e80e 100644
--- a/MeshToolsLib/MeshQuality/AngleSkewMetric.cpp
+++ b/MeshToolsLib/MeshQuality/AngleSkewMetric.cpp
@@ -14,14 +14,12 @@
 
 #include "AngleSkewMetric.h"
 
-#include <boost/math/constants/constants.hpp>
 #include <cmath>
+#include <numbers>
 
 #include "MathLib/MathTools.h"
 #include "MeshLib/Node.h"
 
-using namespace boost::math::double_constants;
-
 namespace MeshToolsLib
 {
 namespace
@@ -30,7 +28,7 @@ template <unsigned long N>
 std::tuple<double, double> getMinMaxAngle(
     std::array<MeshLib::Node, N> const& nodes)
 {
-    double min_angle(two_pi);
+    double min_angle(2 * std::numbers::pi);
     double max_angle(0.0);
 
     for (decltype(N) i = 0; i < N; ++i)
@@ -45,11 +43,11 @@ std::tuple<double, double> getMinMaxAngle(
 
 double checkTriangle(MeshLib::Element const& elem)
 {
+    using namespace std::numbers;
     std::array const nodes = {*elem.getNode(0), *elem.getNode(1),
                               *elem.getNode(2)};
     auto const& [min_angle, max_angle] = getMinMaxAngle(nodes);
-    return std::max((max_angle - third_pi) / two_thirds_pi,
-                    (third_pi - min_angle) / third_pi);
+    return std::max((max_angle - pi / 3) / 2, (pi / 3 - min_angle)) * 3 / pi;
 }
 
 double checkQuad(MeshLib::Element const& elem)
@@ -58,8 +56,8 @@ double checkQuad(MeshLib::Element const& elem)
                               *elem.getNode(2), *elem.getNode(3)};
     auto const& [min_angle, max_angle] = getMinMaxAngle(nodes);
 
-    return std::max((max_angle - half_pi) / half_pi,
-                    (half_pi - min_angle) / half_pi);
+    using namespace std::numbers;
+    return std::max((max_angle - pi / 2), (pi / 2 - min_angle)) * 2 / pi;
 }
 
 double checkTetrahedron(MeshLib::Element const& elem)
@@ -77,8 +75,8 @@ double checkTetrahedron(MeshLib::Element const& elem)
     double const min_angle = *std::min_element(min.begin(), min.end());
     double const max_angle = *std::max_element(max.begin(), max.end());
 
-    return std::max((max_angle - third_pi) / two_thirds_pi,
-                    (third_pi - min_angle) / third_pi);
+    using namespace std::numbers;
+    return std::max((max_angle - pi / 3) / 2, (pi / 3 - min_angle)) * 3 / pi;
 }
 
 double checkHexahedron(MeshLib::Element const& elem)
@@ -96,8 +94,8 @@ double checkHexahedron(MeshLib::Element const& elem)
     double const min_angle = *std::min_element(min.begin(), min.end());
     double const max_angle = *std::max_element(max.begin(), max.end());
 
-    return std::max((max_angle - half_pi) / half_pi,
-                    (half_pi - min_angle) / half_pi);
+    using namespace std::numbers;
+    return std::max((max_angle - pi / 2), (pi / 2 - min_angle)) * 2 / pi;
 }
 
 double checkPrism(MeshLib::Element const& elem)
@@ -117,9 +115,10 @@ double checkPrism(MeshLib::Element const& elem)
     auto const min_angle_tri = std::min(min_angle_tri0, min_angle_tri1);
     auto const max_angle_tri = std::max(max_angle_tri0, max_angle_tri1);
 
-    double const tri_criterion(
-        std::max((max_angle_tri - third_pi) / two_thirds_pi,
-                 (third_pi - min_angle_tri) / third_pi));
+    using namespace std::numbers;
+    double const tri_criterion =
+        std::max((max_angle_tri - pi / 3) / 2, (pi / 3 - min_angle_tri)) * 3 /
+        pi;
 
     std::array<double, 3> min;
     std::array<double, 3> max;
@@ -134,8 +133,9 @@ double checkPrism(MeshLib::Element const& elem)
     double const min_angle_quad = *std::min_element(min.begin(), min.end());
     double const max_angle_quad = *std::max_element(max.begin(), max.end());
 
-    double const quad_criterion(std::max((max_angle_quad - half_pi) / half_pi,
-                                         (half_pi - min_angle_quad) / half_pi));
+    using namespace std::numbers;
+    double const quad_criterion =
+        std::max((max_angle_quad - pi / 2), (pi / 2 - min_angle_quad)) * 2 / pi;
 
     return std::min(tri_criterion, quad_criterion);
 }
diff --git a/MeshToolsLib/MeshSurfaceExtraction.cpp b/MeshToolsLib/MeshSurfaceExtraction.cpp
index c606cc2d54f..763835f7d98 100644
--- a/MeshToolsLib/MeshSurfaceExtraction.cpp
+++ b/MeshToolsLib/MeshSurfaceExtraction.cpp
@@ -14,8 +14,8 @@
 
 #include "MeshSurfaceExtraction.h"
 
-#include <boost/math/constants/constants.hpp>
 #include <memory>
+#include <numbers>
 
 #include "BaseLib/Logging.h"
 #include "MeshLib/Elements/Line.h"
@@ -311,8 +311,7 @@ void MeshSurfaceExtraction::get2DSurfaceElements(
 
     bool const complete_surface = (dir.dot(dir) == 0);
 
-    double const pi(boost::math::constants::pi<double>());
-    double const cos_theta(std::cos(angle * pi / 180.0));
+    double const cos_theta(std::cos(angle * std::numbers::pi / 180.0));
     Eigen::Vector3d const norm_dir(dir.normalized());
 
     for (auto const* elem : all_elements)
diff --git a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
index da04325acb2..07dc174fab9 100644
--- a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
+++ b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
@@ -12,8 +12,8 @@
 
 #pragma once
 
-#include <boost/math/constants/constants.hpp>
 #include <cassert>
+#include <numbers>
 
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/Node.h"
@@ -148,7 +148,7 @@ private:
         // located at edge of the unit triangle, it is possible that
         // r becomes zero.
         auto const r = interpolateZerothCoordinate(shape.N);
-        shape.integralMeasure = boost::math::constants::two_pi<double>() * r;
+        shape.integralMeasure = 2. * std::numbers::pi * r;
     }
 
     const MeshLib::Element* _ele;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
index 089623cf2c4..bf075fd94fe 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
@@ -10,7 +10,7 @@
 
 #include "BHE_1P.h"
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 #include "FlowAndTemperatureControl.h"
 #include "Physics.h"
@@ -112,7 +112,7 @@ void BHE_1P::updateHeatTransferCoefficients(double const flow_rate)
 std::array<double, BHE_1P::number_of_unknowns> BHE_1P::calcThermalResistances(
     double const Nu)
 {
-    constexpr double pi = boost::math::constants::pi<double>();
+    constexpr double pi = std::numbers::pi;
 
     double const lambda_r = refrigerant.thermal_conductivity;
     double const lambda_g = grout.lambda_g;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
index e2d10c6847a..a6df8711e45 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -10,7 +10,7 @@
 
 #include "BHE_1U.h"
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 #include "FlowAndTemperatureControl.h"
 #include "Physics.h"
@@ -173,7 +173,7 @@ void BHE_1U::updateHeatTransferCoefficients(double const flow_rate)
 std::array<double, BHE_1U::number_of_unknowns> BHE_1U::calcThermalResistances(
     double const Nu)
 {
-    constexpr double pi = boost::math::constants::pi<double>();
+    constexpr double pi = std::numbers::pi;
 
     double const lambda_r = refrigerant.thermal_conductivity;
     double const lambda_g = grout.lambda_g;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
index 9034e31d680..4c269046e0b 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
@@ -10,7 +10,7 @@
 
 #include "BHE_2U.h"
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 #include "FlowAndTemperatureControl.h"
 #include "Physics.h"
@@ -193,7 +193,7 @@ void BHE_2U::updateHeatTransferCoefficients(double const flow_rate)
 std::array<double, BHE_2U::number_of_unknowns> BHE_2U::calcThermalResistances(
     double const Nu)
 {
-    constexpr double pi = boost::math::constants::pi<double>();
+    constexpr double pi = std::numbers::pi;
 
     double const lambda_r = refrigerant.thermal_conductivity;
     double const lambda_g = grout.lambda_g;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h b/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
index 0dd6b5d25a2..0b42c799d10 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
@@ -10,7 +10,7 @@
 
 #pragma once
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 namespace BaseLib
 {
@@ -36,11 +36,7 @@ struct BoreholeGeometry
      */
     double const diameter;
 
-    double area() const
-    {
-        constexpr double pi = boost::math::constants::pi<double>();
-        return pi * diameter * diameter / 4;
-    }
+    double area() const { return std::numbers::pi * diameter * diameter / 4; }
 };
 
 BoreholeGeometry createBoreholeGeometry(BaseLib::ConfigTree const& config);
diff --git a/ProcessLib/HeatTransportBHE/BHE/Pipe.h b/ProcessLib/HeatTransportBHE/BHE/Pipe.h
index 54959fb151c..9b0c9824d06 100644
--- a/ProcessLib/HeatTransportBHE/BHE/Pipe.h
+++ b/ProcessLib/HeatTransportBHE/BHE/Pipe.h
@@ -10,7 +10,8 @@
 
 #pragma once
 
-#include <boost/math/constants/constants.hpp>
+#include <cmath>
+#include <numbers>
 
 namespace BaseLib
 {
@@ -38,19 +39,16 @@ struct Pipe
 
     double wallThermalResistance() const
     {
-        constexpr double pi = boost::math::constants::pi<double>();
-
         double const outside_diameter = outsideDiameter();
 
         return std::log(outside_diameter / diameter) /
-               (2.0 * pi * wall_thermal_conductivity);
+               (2.0 * std::numbers::pi * wall_thermal_conductivity);
     }
 
 private:
     double circleArea(double const diameter) const
     {
-        constexpr double pi = boost::math::constants::pi<double>();
-        return pi * diameter * diameter / 4;
+        return std::numbers::pi * diameter * diameter / 4;
     }
 };
 
diff --git a/ProcessLib/HeatTransportBHE/BHE/ThermalResistancesCoaxial.h b/ProcessLib/HeatTransportBHE/BHE/ThermalResistancesCoaxial.h
index e5e6d3f15a3..bb755156934 100644
--- a/ProcessLib/HeatTransportBHE/BHE/ThermalResistancesCoaxial.h
+++ b/ProcessLib/HeatTransportBHE/BHE/ThermalResistancesCoaxial.h
@@ -10,6 +10,8 @@
 
 #pragma once
 
+#include <numbers>
+
 #include "GroutParameters.h"
 #include "Pipe.h"
 #include "RefrigerantProperties.h"
@@ -47,9 +49,10 @@ inline AdvectiveThermalResistanceCoaxial calculateAdvectiveThermalResistance(
     double const hydraulic_diameter =
         coaxialPipesAnnulusDiameter(inner_pipe, outer_pipe);
 
-    auto advective_thermal_resistance = [&](double Nu, double diameter_ratio) {
-        constexpr double pi = boost::math::constants::pi<double>();
-        return 1.0 / (Nu * fluid.thermal_conductivity * pi) * diameter_ratio;
+    auto advective_thermal_resistance = [&](double Nu, double diameter_ratio)
+    {
+        return 1.0 / (Nu * fluid.thermal_conductivity * std::numbers::pi) *
+               diameter_ratio;
     };
     return {advective_thermal_resistance(Nu_inner_pipe, 1.),
             advective_thermal_resistance(
@@ -70,8 +73,6 @@ calculateGroutAndGroutSoilExchangeThermalResistance(
     Pipe const& outer_pipe, GroutParameters const& grout_parameters,
     double const borehole_diameter)
 {
-    constexpr double pi = boost::math::constants::pi<double>();
-
     double const outer_pipe_outside_diameter = outer_pipe.outsideDiameter();
     double const chi =
         std::log(std::sqrt(borehole_diameter * borehole_diameter +
@@ -81,7 +82,7 @@ calculateGroutAndGroutSoilExchangeThermalResistance(
         std::log(borehole_diameter / outer_pipe_outside_diameter);
     double const R_g =
         std::log(borehole_diameter / outer_pipe_outside_diameter) / 2 /
-        (pi * grout_parameters.lambda_g);
+        (std::numbers::pi * grout_parameters.lambda_g);
     double const conductive_b = chi * R_g;
     double const grout_soil = (1 - chi) * R_g;
     return {conductive_b, grout_soil};
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
index 56b62e26e96..79e6b0a4aa4 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -10,6 +10,8 @@
  */
 #pragma once
 
+#include <numbers>
+
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
 
@@ -238,7 +240,7 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
             case PhaseFieldModel::COHESIVE:
             {
                 auto const local_Jac_COHESIVE =
-                    (2.0 / boost::math::double_constants::pi * gc *
+                    (2.0 / std::numbers::pi * gc *
                      (-N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
                      w)
                         .eval();
@@ -381,7 +383,7 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::computeEnergy(
             case PhaseFieldModel::COHESIVE:
             {
                 element_surface_energy +=
-                    gc / boost::math::double_constants::pi *
+                    gc / std::numbers::pi *
                     ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
                     w;
                 break;
diff --git a/Tests/GeoLib/TestEquidistantPointGeneration.cpp b/Tests/GeoLib/TestEquidistantPointGeneration.cpp
index 20fabb3e6e0..49202b9514e 100644
--- a/Tests/GeoLib/TestEquidistantPointGeneration.cpp
+++ b/Tests/GeoLib/TestEquidistantPointGeneration.cpp
@@ -10,7 +10,7 @@
 
 #include <gtest/gtest.h>
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 #include <random>
 
 #include "GeoLib/Utils.h"
@@ -22,7 +22,7 @@ struct GeoLibGenerateEquidistantPoints : public testing::Test
         // enable random engine
         std::random_device rd;
         std::mt19937 random_engine_mt19937(rd());
-        constexpr auto pi = boost::math::double_constants::pi;
+        constexpr auto pi = std::numbers::pi;
         std::normal_distribution<> normal_dist_phi(-pi, pi);  // azimuthal angle
         std::normal_distribution<> normal_dist_theta(0, pi);  // polar angle
         // generate point on unit sphere
diff --git a/Tests/GeoLib/TestSurfaceIsPointInSurface.cpp b/Tests/GeoLib/TestSurfaceIsPointInSurface.cpp
index d1b81f509c6..f617335d878 100644
--- a/Tests/GeoLib/TestSurfaceIsPointInSurface.cpp
+++ b/Tests/GeoLib/TestSurfaceIsPointInSurface.cpp
@@ -12,8 +12,8 @@
 #include <gtest/gtest.h>
 
 #include <array>
-#include <boost/math/constants/constants.hpp>
 #include <memory>
+#include <numbers>
 #include <random>
 #include <vector>
 
@@ -114,8 +114,7 @@ TEST(GeoLib, SurfaceIsPointInSurface)
                                                            n_steps, f));
 
         // random rotation angles
-        std::normal_distribution<> normal_dist_angles(
-            0, boost::math::double_constants::two_pi);
+        std::normal_distribution<> normal_dist_angles(0, 2. * std::numbers::pi);
         std::array<double, 3> euler_angles = {
             {normal_dist_angles(random_engine_mt19937),
              normal_dist_angles(random_engine_mt19937),
diff --git a/Tests/MaterialLib/TestGasPressureDependentPermeability.cpp b/Tests/MaterialLib/TestGasPressureDependentPermeability.cpp
index 8ea5b381a7f..a35c3cc7d7f 100644
--- a/Tests/MaterialLib/TestGasPressureDependentPermeability.cpp
+++ b/Tests/MaterialLib/TestGasPressureDependentPermeability.cpp
@@ -12,7 +12,6 @@
 #include <gtest/gtest.h>
 
 #include <Eigen/Core>
-#include <boost/math/constants/constants.hpp>
 
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Properties/CreateGasPressureDependentPermeability.h"
diff --git a/Tests/MaterialLib/TestStrainDependentPermeability.cpp b/Tests/MaterialLib/TestStrainDependentPermeability.cpp
index afacff34782..529ec131c67 100644
--- a/Tests/MaterialLib/TestStrainDependentPermeability.cpp
+++ b/Tests/MaterialLib/TestStrainDependentPermeability.cpp
@@ -12,7 +12,6 @@
 #include <gtest/gtest.h>
 
 #include <Eigen/Core>
-#include <boost/math/constants/constants.hpp>
 
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Properties/CreateStrainDependentPermeability.h"
diff --git a/Tests/MeshLib/TestTetQualityCriteria.cpp b/Tests/MeshLib/TestTetQualityCriteria.cpp
index 4af1e19668f..a80bc0999cb 100644
--- a/Tests/MeshLib/TestTetQualityCriteria.cpp
+++ b/Tests/MeshLib/TestTetQualityCriteria.cpp
@@ -10,9 +10,9 @@
 #include <gtest/gtest.h>
 
 #include <Eigen/Core>
-#include <boost/math/constants/constants.hpp>
 #include <cmath>
 #include <memory>
+#include <numbers>
 #include <numeric>
 #include <random>
 
@@ -148,7 +148,7 @@ std::pair<double, double> getMinMaxAngleFromTriangle(Eigen::Vector3d n0,
 double computeCriterionForTet(Eigen::Vector3d n0, Eigen::Vector3d n1,
                               Eigen::Vector3d n2, Eigen::Vector3d n3)
 {
-    using namespace boost::math::double_constants;
+    using namespace std::numbers;
     auto const alpha0 = getMinMaxAngleFromTriangle(n0, n2, n1);
     auto const alpha1 = getMinMaxAngleFromTriangle(n0, n1, n3);
     auto const alpha2 = getMinMaxAngleFromTriangle(n1, n2, n3);
@@ -157,8 +157,7 @@ double computeCriterionForTet(Eigen::Vector3d n0, Eigen::Vector3d n1,
         std::min({alpha0.first, alpha1.first, alpha2.first, alpha3.first});
     auto const max =
         std::max({alpha0.second, alpha1.second, alpha2.second, alpha3.second});
-    return std::max((max - third_pi) / two_thirds_pi,
-                    (third_pi - min) / third_pi);
+    return std::max((max - pi / 3) * 3 / (2 * pi), (pi / 3 - min) * 3 / pi);
 }
 
 TEST_F(TetElementQuality, EquiAngleSkew)
diff --git a/Tests/MeshLib/TestTriQualityCriteria.cpp b/Tests/MeshLib/TestTriQualityCriteria.cpp
index 245db1d0e5a..7f1da3c848a 100644
--- a/Tests/MeshLib/TestTriQualityCriteria.cpp
+++ b/Tests/MeshLib/TestTriQualityCriteria.cpp
@@ -9,9 +9,9 @@
 
 #include <gtest/gtest.h>
 
-#include <boost/math/constants/constants.hpp>
 #include <cmath>
 #include <memory>
+#include <numbers>
 #include <numeric>
 #include <random>
 
@@ -103,7 +103,7 @@ TEST_F(TriElementQuality, EdgeRatio)
 
 TEST_F(TriElementQuality, EquiAngleSkew)
 {
-    using namespace boost::math::double_constants;
+    using namespace std::numbers;
     auto const type = MeshLib::MeshQualityType::EQUIANGLESKEW;
     auto const element_quality_vector =
         getElementQualityVectorFromRegularTriMesh(type);
@@ -112,11 +112,11 @@ TEST_F(TriElementQuality, EquiAngleSkew)
                                        std::pow(1.0 / n_subdivisions[1], 2));
     std::array const angles = {
         std::asin((1.0 / n_subdivisions[0]) / hypothenuse),
-        std::asin((1.0 / n_subdivisions[1]) / hypothenuse), half_pi};
+        std::asin((1.0 / n_subdivisions[1]) / hypothenuse), pi / 2};
     auto const& min_max = std::minmax_element(angles.begin(), angles.end());
     auto const expected_value =
-        std::max((*min_max.second - third_pi) / two_thirds_pi,
-                 (third_pi - *min_max.first) / third_pi);
+        std::max((*min_max.second - pi / 3) * 3 / (2 * pi),
+                 (pi / 3 - *min_max.first) * 3 / pi);
 
     for (auto const element_quality : element_quality_vector)
     {
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index 2453a4a66a7..98d7d61b869 100644
--- a/Tests/NumLib/ODEs.h
+++ b/Tests/NumLib/ODEs.h
@@ -10,7 +10,7 @@
 
 #pragma once
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MathLib/LinAlg/UnifiedMatrixSetters.h"
@@ -142,8 +142,7 @@ public:
 
 const double ODETraits<ODE1>::t0 = 0.0;
 
-const double ODETraits<ODE1>::t_end =
-    2.0 * 3.141592653589793238462643383279502884197169399375105820974944;
+const double ODETraits<ODE1>::t_end = 2. * std::numbers::pi;
 // ODE 1 end //////////////////////////////////////////////////////
 
 // ODE 2 //////////////////////////////////////////////////////////
@@ -364,6 +363,5 @@ public:
 
 const double ODETraits<ODE3>::t0 = 0.0;
 
-const double ODETraits<ODE3>::t_end =
-    0.5 * boost::math::constants::pi<double>();
+const double ODETraits<ODE3>::t_end = std::numbers::pi / 2.;
 // ODE 3 end //////////////////////////////////////////////////////
-- 
GitLab