diff --git a/src/ExternalCurrent.cpp b/src/ExternalCurrent.cpp
index 7ddd0ca0caea198f5e2acc0e952fa0b26e999c3b..361f829ef214beb3af99960acf2e1dd0c23fa6e9 100644
--- a/src/ExternalCurrent.cpp
+++ b/src/ExternalCurrent.cpp
@@ -1,4 +1,7 @@
 #include "ExternalCurrent.hpp"
+
+#include <numbers>
+
 #include "GlobalDefinitions.hpp"
 
 std::function<double(double,double,double,double)> excitationFunc=IextStepfct;
@@ -17,7 +20,7 @@ double IextTanh(double t, double amp, double exc, double dur) {
 double IextArcTan(double t, double amp, double exc, double dur) {
   double start_t = 4 * (t -  exc);
   double end_t = 4 * (t - (exc + dur));
-  return (atan(start_t) - atan(end_t)) / Pi * amp;
+  return (atan(start_t) - atan(end_t)) / std::numbers::pi * amp;
 }
 
 void SelectExternalCurrentFunction(const std::string &externalCurrentFunction) {
diff --git a/src/elasticity/problem/DynamicBoundaryProblems.cpp b/src/elasticity/problem/DynamicBoundaryProblems.cpp
index e10228d3eb5de7d3d20e56c05bb98f1712a41815..768e619eb8b389ebc45e9cf2fc726dac37a6b88c 100644
--- a/src/elasticity/problem/DynamicBoundaryProblems.cpp
+++ b/src/elasticity/problem/DynamicBoundaryProblems.cpp
@@ -1,5 +1,7 @@
 #include "DynamicBoundaryProblems.hpp"
 
+#include <numbers>
+
 VectorField DExponentialStiffnessProblem::Deformation(double time, const Point &x) const {
     VectorField  u(
             exp(x[0]), exp(x[1]), spaceDim > 2 ? exp(x[2]) : 0.0
@@ -211,46 +213,59 @@ double DPolynomialNeumannProblem::pressure(double t, const Point &x, int bndCond
   }
 }
 
-
-VectorField DExponentialNeumannProblem::Deformation(double time, const Point &x) const {
-  VectorField u(cos(0.5 * Pi * x[0]), cos(0.5 * Pi * x[1]),cos(0.5 * Pi * x[2]));
-  return -2.0 / Pi * sin(0.5 * Pi * time) * u;
+VectorField DExponentialNeumannProblem::Deformation(double time,
+                                                    const Point &x) const {
+  using std::numbers::pi;
+  VectorField u(cos(0.5 * pi * x[0]), cos(0.5 * pi * x[1]),
+                cos(0.5 * pi * x[2]));
+  return -2.0 / pi * sin(0.5 * pi * time) * u;
 }
 
-VectorField DExponentialNeumannProblem::Velocity(double time, const Point &x) const {
-  VectorField u(cos(0.5 * Pi * x[0]), cos(0.5 * Pi * x[1]), cos(0.5 * Pi * x[2]));
-  return -cos(0.5 * Pi * time) * u;
+VectorField DExponentialNeumannProblem::Velocity(double time,
+                                                 const Point &x) const {
+  using std::numbers::pi;
+  VectorField u(cos(0.5 * pi * x[0]), cos(0.5 * pi * x[1]),
+                cos(0.5 * pi * x[2]));
+  return -cos(0.5 * pi * time) * u;
 }
 
-VectorField DExponentialNeumannProblem::Acceleration(double time, const Point &x) const {
-  VectorField u(cos(0.5 * Pi * x[0]), cos(0.5 * Pi * x[1]), cos(0.5 * Pi * x[2]) );
-  return 0.5 * Pi * sin(0.5 * Pi * time) * u;
+VectorField DExponentialNeumannProblem::Acceleration(double time,
+                                                     const Point &x) const {
+  using std::numbers::pi;
+  VectorField u(cos(0.5 * pi * x[0]), cos(0.5 * pi * x[1]),
+                cos(0.5 * pi * x[2]));
+  return 0.5 * pi * sin(0.5 * pi * time) * u;
 }
 
-
-Tensor DExponentialNeumannProblem::DeformationGradient(double time, const Point &x) const {
-  Tensor Du(- 0.5 * Pi * sin(0.5 * Pi * x[0]), 0.0, 0.0,
-            0.0, - 0.5 * Pi * sin(0.5 * Pi * x[1]), 0.0,
-            0.0, 0.0, - 0.5 * Pi * sin(0.5 * Pi * x[2])
-  );
-  return -2.0 / Pi * sin(0.5 * Pi * time) * Du;
+Tensor DExponentialNeumannProblem::DeformationGradient(double time,
+                                                       const Point &x) const {
+  using std::numbers::pi;
+  Tensor Du(-0.5 * pi * sin(0.5 * pi * x[0]), 0.0, 0.0, 0.0,
+            -0.5 * pi * sin(0.5 * pi * x[1]), 0.0, 0.0, 0.0,
+            -0.5 * pi * sin(0.5 * pi * x[2]));
+  return -2.0 / pi * sin(0.5 * pi * time) * Du;
 }
 
-VectorField DExponentialNeumannProblem::Load(double time, const Point &x, const Cell &c) const {
+VectorField DExponentialNeumannProblem::Load(double time, const Point &x,
+                                             const Cell &c) const {
+  using std::numbers::pi;
   VectorField a = Acceleration(time, x);
-  VectorField u(-0.25 * Pi * Pi * cos(0.5 * Pi * x[0]), -0.25 * Pi * Pi * cos(0.5 * Pi * x[1]),
-                 -0.25 * Pi * Pi * cos(0.5 * Pi * x[2]));
+  VectorField u(-0.25 * pi * pi * cos(0.5 * pi * x[0]),
+                -0.25 * pi * pi * cos(0.5 * pi * x[1]),
+                -0.25 * pi * pi * cos(0.5 * pi * x[2]));
 
-  return a - ((-2.0 / Pi * sin(0.5 * Pi * time)) * u);
+  return a - ((-2.0 / pi * sin(0.5 * pi * time)) * u);
 }
 
-double DExponentialNeumannProblem::pressure(double t, const Point &x, int bndCondition) const {
+double DExponentialNeumannProblem::pressure(double t, const Point &x,
+                                            int bndCondition) const {
+  using std::numbers::pi;
   if (bndCondition == 230) {
-    return  - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[0]) ;
+    return -sin(0.5 * pi * t) * sin(0.5 * pi * x[0]);
   } else if (bndCondition == 231) {
-    return  - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[1]);
+    return -sin(0.5 * pi * t) * sin(0.5 * pi * x[1]);
   } else if (bndCondition == 232) {
-    return  - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[2]) ;
+    return -sin(0.5 * pi * t) * sin(0.5 * pi * x[2]);
   } else {
     return 0;
   }
diff --git a/src/electrophysiology/problem/DiffusionProblem.cpp b/src/electrophysiology/problem/DiffusionProblem.cpp
index b66cddef6f4ba0ce1e89c3fc7909da4c568fb945..1df6419a78338d87c2e4c3150ef9272a5a994db0 100644
--- a/src/electrophysiology/problem/DiffusionProblem.cpp
+++ b/src/electrophysiology/problem/DiffusionProblem.cpp
@@ -1,5 +1,6 @@
 #include "DiffusionProblem.hpp"
 
+#include <numbers>
 #include <utility>
 
 DiffusionProblem::DiffusionProblem(std::array<double, 9> entries) :
@@ -38,7 +39,9 @@ DiffusionExampleOne::DiffusionExampleOne() : DiffusionProblem(
     {1, 0.5, 0.5, 0.5, 1.0, 0.5, 0.5, 0.5, 1.0}) {}
 
 double DiffusionExampleOne::Potential(double time, const Point &x) const {
-  return 0.25 * (1 - sin(time) * sin(time)) * (1 - cos(2 * Pi * x[0])) * (1 - cos(2 * Pi * x[1]));
+  using std::numbers::pi;
+  return 0.25 * (1 - sin(time) * sin(time)) * (1 - cos(2 * pi * x[0])) *
+         (1 - cos(2 * pi * x[1]));
 }
 
 double DiffusionExampleOne::IonCurrent(double time, const Point &x) const {
@@ -48,20 +51,23 @@ double DiffusionExampleOne::IonCurrent(double time, const Point &x) const {
 }
 
 double DiffusionExampleOne::timeDerivative(double time, const Point &x) const {
-  return -0.5 * (sin(time) * cos(time)) * (1 - cos(2 * Pi * x[0])) * (1 - cos(2 * Pi * x[1]));
+  using std::numbers::pi;
+  return -0.5 * (sin(time) * cos(time)) * (1 - cos(2 * pi * x[0])) *
+         (1 - cos(2 * pi * x[1]));
 }
 
 double DiffusionExampleOne::divergence(int row, const Point &x) const {
+  using std::numbers::pi;
   int i = row;
   int j = (row + 1) % SpaceDimension;
   int k = (row + 2) % SpaceDimension;
-  return Pi * Pi *
-         (conductivity[i][i] * cos(2 * Pi * x[i]) * (1 - cos(2 * Pi * x[j])) *
-          (1 - cos(2 * Pi * x[k]))) *
-         (conductivity[i][j] * sin(2 * Pi * x[j]) * (1 + sin(2 * Pi * x[i])) *
-          (1 - cos(2 * Pi * x[k]))) *
-         (conductivity[i][k] * sin(2 * Pi * x[k]) * (1 + sin(2 * Pi * x[i])) *
-          (1 - cos(2 * Pi * x[j])));
+  return pi * pi *
+         (conductivity[i][i] * cos(2 * pi * x[i]) * (1 - cos(2 * pi * x[j])) *
+          (1 - cos(2 * pi * x[k]))) *
+         (conductivity[i][j] * sin(2 * pi * x[j]) * (1 + sin(2 * pi * x[i])) *
+          (1 - cos(2 * pi * x[k]))) *
+         (conductivity[i][k] * sin(2 * pi * x[k]) * (1 + sin(2 * pi * x[i])) *
+          (1 - cos(2 * pi * x[j])));
 }
 
 DiffusionExampleTwo::DiffusionExampleTwo() : DiffusionProblem(
@@ -71,7 +77,7 @@ double DiffusionExampleTwo::Potential(double time, const Point &x) const {
   double v =
       256 * (1 - sin(time) * sin(time)) * x[0] * x[0] * (x[0] - 1) * (x[0] - 1) * x[1] * x[1] *
       (x[1] - 1) * (x[1] - 1) * x[2] * x[2] * (x[2] - 1) * (x[2] - 1);
-  double g = 0.5 * (1 + cos(Pi * (x[1] - 0.5)));
+  double g = 0.5 * (1 + cos(std::numbers::pi * (x[1] - 0.5)));
   return v * g;
 }
 
diff --git a/test/cellmodels/MTest.hpp b/test/cellmodels/MTest.hpp
index 84d0ecb4862a128500c781b973d5f098af2ef824..19783b10ac8ebec52db969259b6f22827b1da0da 100644
--- a/test/cellmodels/MTest.hpp
+++ b/test/cellmodels/MTest.hpp
@@ -1,6 +1,8 @@
 #ifndef MTEST_HPP
 #define MTEST_HPP
 
+#include <numbers>
+
 #include "MCellModel.hpp"
 
 class MTest : public MElphyModel {
@@ -56,11 +58,11 @@ public:
   MExact() : MTest() {}
 
   double Value(double t, const std::vector<double> &vw) const override {
-    return sin(2 * Pi * t);
+    return sin(2 * std::numbers::pi * t);
   }
 
   double Diff(double t, const std::vector<double> &vw) const override {
-    return 2 * Pi * cos(2 * Pi * t);
+    return 2 * std::numbers::pi * cos(2 * std::numbers::pi * t);
   }
 
   void InitialValues(double t, std::vector<double> &vw) override {
diff --git a/test/elasticity/TestCirculation.cpp b/test/elasticity/TestCirculation.cpp
index b28d92a1007d5318c3d74e52d17affe28e421246..82061aa107142079f206f313c42f1bfc167013f4 100644
--- a/test/elasticity/TestCirculation.cpp
+++ b/test/elasticity/TestCirculation.cpp
@@ -1,3 +1,5 @@
+#include <numbers>
+
 #include <GlobalDefinitions.hpp>
 #include "pressure_bnd/MCircModel.hpp"
 #include "TestEnvironmentCardMech.hpp"
@@ -10,6 +12,8 @@ TEST(CirculationTest, SingleStep) {
 
 
 TEST(CirculationTest, Periodic){
+  using std::numbers::pi;
+
   MCircModel circModel;
   TimeSeries TS(0.0, 1.0, 0.001);
 
@@ -25,23 +29,23 @@ TEST(CirculationTest, Periodic){
     if (periodicTime >= 1.0) periodicTime -= 1.0;
 
     if (periodicTime < 0.025) {
-      pressure[2] = pressure[3] = 5.0 * (-cos(0.5 * Pi / 0.025 * periodicTime) + 1);
+      pressure[2] = pressure[3] = 5.0 * (-cos(0.5 * pi / 0.025 * periodicTime) + 1);
     } else if (periodicTime < 0.05) {
-      pressure[2] = pressure[3] = 5.0 * (cos(0.5 * Pi / 0.025 * (periodicTime - 0.025)) + 1);
+      pressure[2] = pressure[3] = 5.0 * (cos(0.5 * pi / 0.025 * (periodicTime - 0.025)) + 1);
     } else if (periodicTime > 0.35) {
       pressure[2] = pressure[3] = 0.0;
     } else if (periodicTime > 0.3) {
-      pressure[2] = pressure[3] = 5.0 * (cos(0.5 * Pi / 0.05 * (periodicTime - 0.3)) + 1);
+      pressure[2] = pressure[3] = 5.0 * (cos(0.5 * pi / 0.05 * (periodicTime - 0.3)) + 1);
     } else if (periodicTime > 0.15) {
-      pressure[2] = pressure[3] = 5.0 * (-cos(0.5 * Pi / 0.15 * (periodicTime - 0.15)) + 1);
+      pressure[2] = pressure[3] = 5.0 * (-cos(0.5 * pi / 0.15 * (periodicTime - 0.15)) + 1);
     } else {
       pressure[2] = pressure[3] = 0.0;
     }
 
     if (periodicTime < 0.15) {
-      pressure[0] = pressure[1] = 60 * (-cos(0.5 * Pi / 0.15 * periodicTime) + 1);
+      pressure[0] = pressure[1] = 60 * (-cos(0.5 * pi / 0.15 * periodicTime) + 1);
     } else if (periodicTime < 0.3) {
-      pressure[0] = pressure[1] = 60.0 * (cos(0.5 * Pi / 0.15 * (periodicTime - 0.15)) + 1);
+      pressure[0] = pressure[1] = 60.0 * (cos(0.5 * pi / 0.15 * (periodicTime - 0.15)) + 1);
     } else pressure[0] = pressure[1] = 0.0;
 
     std::cout << "Volumes at t=" << t << " : " << std::endl;