diff --git a/tests/unitTests/test.cpp b/tests/unitTests/test.cpp index 1de72885..83945553 100644 --- a/tests/unitTests/test.cpp +++ b/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(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(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(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]; diff --git a/tests/unitTests/test.h b/tests/unitTests/test.h index 990b8cab..2c409a9e 100644 --- a/tests/unitTests/test.h +++ b/tests/unitTests/test.h @@ -47,6 +47,8 @@ #include #include +#include +#include #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 distribution(min, max); + return distribution(engine); +} + // Creates a random deformation gradient F with values between min and max, // and det(F) > 0 template @@ -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(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(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 } }