Switch to C++ random library, reduce finite difference deltas to 1e-7

tests/unitTests/test.cpp

@@ -46,7 +46,7 @@ class NeoHookeanTest : public ::testing::Test {
  int n_iter = 10; // Number of random perturbations to test
  double rel_tol = 1e-3; // relative tolerance for comparing values
  double abs_tol = 1e-11; // absolute tolerance for comparing values
- double delta = 1e-6; // perturbation scaling factor
+ double delta = 1e-7; // perturbation scaling factor
  bool verbose = false; // Show values of S, dE, SdE and dPsi
 
  // Add the test object
@@ -55,11 +55,8 @@ class NeoHookeanTest : public ::testing::Test {
  // Setup method to initialize variables before each test
  void SetUp() override {
 
- // Seed random number generator
- //srand(static_cast<unsigned int>(time(0)));
-
- // Set random values for the Neo-Hookean parameter between 1000 and 10000
- params.C10 = 1000 + 9000 * (double)rand() / RAND_MAX;
+ // Set random values for the Neo-Hookean parameters between 1000 and 10000
+ params.C10 = getRandomDouble(1000.0, 10000.0);
 
  // Initialize the test object
  TestNH = new TestNeoHookean(params);
@@ -130,7 +127,7 @@ class MooneyRivlinTest : public ::testing::Test {
  int n_iter = 10; // Number of random perturbations to test
  double rel_tol = 1e-3; // relative tolerance for comparing dPsi and dS with values from svFSI
  double abs_tol = 1e-11; // absolute tolerance for comparing values
- double delta = 1e-6; // perturbation scaling factor
+ double delta = 1e-7; // perturbation scaling factor
  bool verbose = false; // Show values of S, dE, SdE and dPsi
 
 
@@ -140,12 +137,9 @@ class MooneyRivlinTest : public ::testing::Test {
  // Setup method to initialize variables before each test
  void SetUp() override {
 
- // Seed random number generator
- //srand(static_cast<unsigned int>(time(0)));
-
  // Set random values for the Mooney-Rivlin parameters between 1000 and 10000
- params.C01 = 1000 + 9000 * (double)rand() / RAND_MAX;
- params.C10 = 1000 + 9000 * (double)rand() / RAND_MAX;
+ params.C01 = getRandomDouble(1000.0, 10000.0);
+ params.C10 = getRandomDouble(1000.0, 10000.0);
 
  // Initialize the test object
  TestMR = new TestMooneyRivlin(params);
@@ -221,7 +215,7 @@ class HolzapfelOgdenTest : public ::testing::Test {
  int n_iter = 10; // Number of random perturbations to test
  double rel_tol = 1e-3; // relative tolerance for comparing dPsi and dS with values from svFSI
  double abs_tol = 1e-11; // absolute tolerance for comparing values
- double delta = 1e-6; // perturbation scaling factor
+ double delta = 1e-7; // perturbation scaling factor
  bool verbose = false; // Show values of S, dE, SdE and dPsi
 
  // Add the test object
@@ -230,9 +224,6 @@ class HolzapfelOgdenTest : public ::testing::Test {
  // Setup method to initialize variables before each test
  void SetUp() override {
 
- // Seed random number generator
- //srand(static_cast<unsigned int>(time(0)));
-
  // Set Holzapfel-Ogden parameters from cardiac benchmark paper
  params.a = 59.0; // Pa
  params.a_f = 18472.0; // Pa
@@ -245,14 +236,14 @@ class HolzapfelOgdenTest : public ::testing::Test {
  params.kappa = 1.0e6; // Pa
  params.k = 100000.0; // Pa 
 
- // Set random values for f and normalize
- params.f[0] = (double)rand() / RAND_MAX;
- params.f[1] = (double)rand() / RAND_MAX;
- params.f[2] = (double)rand() / RAND_MAX;
+ // Set random values for f between 0 and 1 and normalize
+ params.f[0] = getRandomDouble(0.0, 1.0);
+ params.f[1] = getRandomDouble(0.0, 1.0);
+ params.f[2] = getRandomDouble(0.0, 1.0);
  double norm_f = sqrt(params.f[0]*params.f[0] + params.f[1]*params.f[1] + params.f[2]*params.f[2]);
  params.f[0] /= norm_f; params.f[1] /= norm_f; params.f[2] /= norm_f;
 
- // Create random orthogonal s
+ // Create s orthogonal to f
  if (fabs(params.f[0]) < 0.9) { // Check if f[0] is not the dominant component
  params.s[0] = 0;
  params.s[1] = params.f[2];

tests/unitTests/test.h

@@ -47,6 +47,8 @@
 
 #include <stdlib.h>
 #include <iostream>
+#include <random>
+#include <chrono>
 #include "gtest/gtest.h" // include GoogleTest
 #include "mat_fun.h"
 #include "mat_fun_carray.h"
@@ -136,6 +138,16 @@ void create_identity_F(double F[N][N]) {
  }
 }
 
+// Inline getRandomDouble function
+inline double getRandomDouble(double min, double max) {
+ // Uncomment to use a random seed
+ //unsigned int seed = std::chrono::system_clock::now().time_since_epoch().count();
+ unsigned int seed = 42;
+ static std::default_random_engine engine(seed);
+ std::uniform_real_distribution<double> distribution(min, max);
+ return distribution(engine);
+}
+
 // Creates a random deformation gradient F with values between min and max,
 // and det(F) > 0
 template<int N>
@@ -146,7 +158,7 @@ void create_random_F(double F[N][N], double min=0.1, double max=10.0) {
  while (J < 0) {
  for (int i = 0; i < N; i++) {
  for (int J = 0; J < N; J++) {
- F[i][J] = min + (max - min) * rand() / RAND_MAX;
+ F[i][J] = getRandomDouble(min, max);
  }
  }
  J = mat_fun_carray::mat_det<N>(F);
@@ -162,7 +174,7 @@ void create_random_perturbed_identity_F(double F[N][N], double max_deviation) {
  while (J < 0) {
  for (int i = 0; i < N; i++) {
  for (int J = 0; J < N; J++) {
- F[i][J] = (i == J) + max_deviation * (2.0 * rand() / RAND_MAX - 1.0);
+ F[i][J] = (i == J) + max_deviation * getRandomDouble(-1.0, 1.0);
  }
  }
  J = mat_fun_carray::mat_det<N>(F);
@@ -178,7 +190,7 @@ void perturb_random_F(const double F[N][N], const double delta, double F_tilde[N
  double dF_iJ;
  for (int i = 0; i < N; i++) {
  for (int J = 0; J < N; J++) {
- dF_iJ = delta * (2.0 * std::rand() / RAND_MAX - 1.0); // (random number between -1 and 1) * delta
+ dF_iJ = delta * getRandomDouble(-1.0, 1.0);
  F_tilde[i][J] = F[i][J] + dF_iJ; // perturbed deformation gradient
  }
  }