diff --git a/Code/Source/svFSI/mat_models.cpp b/Code/Source/svFSI/mat_models.cpp index 4b706443..b6cdf6a1 100644 --- a/Code/Source/svFSI/mat_models.cpp +++ b/Code/Source/svFSI/mat_models.cpp @@ -137,11 +137,23 @@ void get_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibStrsT } } -/// @brief Compute 2nd Piola-Kirchhoff stress and material stiffness tensors -/// including both dilational and isochoric components. -/// -/// Reproduces the Fortran 'GETPK2CC' subroutine. -// +/** + * @brief Compute 2nd Piola-Kirchhoff stress and material stiffness tensors + * including both dilational and isochoric components. + * + * Reproduces the Fortran 'GETPK2CC' subroutine. + * + * @param[in] com_mod Object containing global common variables. + * @param[in] cep_mod Object containing electrophysiology-specific common variables. + * @param[in] lDmn Domain object. + * @param[in] F Deformation gradient tensor. + * @param[in] nfd Number of fiber directions. + * @param[in] fl Fiber directions. + * @param[in] ya Electrophysiology active stress. + * @param[out] S 2nd Piola-Kirchhoff stress tensor (modified in place). + * @param[out] Dm Material stiffness tensor (modified in place). + * @return None, but modifies S and Dm in place. + */ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Array& F, const int nfd, const Array& fl, const double ya, Array& S, Array& Dm) { @@ -479,10 +491,23 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn cc_to_voigt(nsd, CC, Dm); } -/// @brief Compute isochoric (deviatoric) component of 2nd Piola-Kirchhoff stress and material stiffness tensors. -/// -/// Reproduces 'SUBROUTINE GETPK2CCdev(lDmn, F, nfd, fl, ya, S, Dm, Ja)'. -// +/** + * @brief Compute isochoric (deviatoric) component of 2nd Piola-Kirchhoff stress and material stiffness tensors. + * + * Reproduces 'SUBROUTINE GETPK2CCdev(lDmn, F, nfd, fl, ya, S, Dm, Ja)'. + * + * @param[in] com_mod Object containing global common variables. + * @param[in] cep_mod Object containing electrophysiology-specific common variables. + * @param[in] lDmn Domain object. + * @param[in] F Deformation gradient tensor. + * @param[in] nfd Number of fiber directions. + * @param[in] fl Fiber directions. + * @param[in] ya Electrophysiology active stress. + * @param[out] S 2nd Piola-Kirchhoff stress tensor (isochoric part). + * @param[out] Dm Material stiffness tensor (isochoric part). + * @param[out] Ja Jacobian for active strain. + * @return None, but modifies S, Dm, and Ja in place. + */ void get_pk2cc_dev(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Array& F, const int nfd, const Array& fl, const double ya, Array& S, Array& Dm, double& Ja) { @@ -765,7 +790,7 @@ void get_pk2cc_dev(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& } double r1 = J2d*mat_ddot(C, Sb, nsd) / nd; - auto S = J2d*Sb - r1*Ci; + S = J2d*Sb - r1*Ci; auto PP = ten_ids(nsd) - (1.0/nd) * ten_dyad_prod(Ci, C, nsd); CC = ten_ddot(CCb, PP, nsd); @@ -1409,6 +1434,22 @@ void get_tau(const ComMod& com_mod, const dmnType& lDmn, const double detF, cons tauC = ctC * (he*c) * (rho0/detF); } +/** + * @brief Compute rho, beta, drho/dp, dbeta/dp for volumetric penalty terms in + * the ustruct formulation. + * + * See ustruct paper (https://doi.org/10.1016/j.cma.2018.03.045) Section 2.4. + * + * @param[in] com_mod Object containing global common variables + * @param[in] lDmn Domain object + * @param[in] p Pressure. + * @param[out] ro Solid density, rho. + * @param[out] bt Isothermal compressibility coefficient, beta. + * @param[out] dro Derivative of rho with respect to p. + * @param[out] dbt Derivative of beta with respect to p. + * @param[out] Ja Active strain Jacobian. + * @return None, but updates ro, bt, dro, and dbt. + */ void g_vol_pen(const ComMod& com_mod, const dmnType& lDmn, const double p, double& ro, double& bt, double& dro, double& dbt, const double Ja) { diff --git a/Code/Source/svFSI/mat_models_carray.h b/Code/Source/svFSI/mat_models_carray.h index a5010679..8946fcd6 100644 --- a/Code/Source/svFSI/mat_models_carray.h +++ b/Code/Source/svFSI/mat_models_carray.h @@ -160,10 +160,24 @@ void voigt_to_cc_carray(const double Dm[2*N][2*N], double CC[N][N][N][N]) { } } -//----------- -// get_pk2cc -//----------- -// +/** + * @brief Compute 2nd Piola-Kirchhoff stress and material stiffness tensors + * including both dilational and isochoric components. + * + * Reproduces the Fortran 'GETPK2CC' subroutine. + * + * @tparam N Template parameter for number of spatial dimensions. + * @param[in] com_mod Object containing global common variables. + * @param[in] cep_mod Object containing electrophysiology-specific common variables. + * @param[in] lDmn Domain object. + * @param[in] F Deformation gradient tensor. + * @param[in] nfd Number of fiber directions. + * @param[in] fl Fiber directions. + * @param[in] ya Electrophysiology active stress. + * @param[out] S 2nd Piola-Kirchhoff stress tensor (modified in place). + * @param[out] Dm Material stiffness tensor (modified in place). + * @return None, but modifies S and Dm in place. + */ template void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const double F[N][N], const int nfd, const Array& fl, const double ya, double S[N][N], double Dm[2*N][2*N]) diff --git a/Code/Source/svFSI/ustruct.cpp b/Code/Source/svFSI/ustruct.cpp index 3ddd9744..25ab4162 100644 --- a/Code/Source/svFSI/ustruct.cpp +++ b/Code/Source/svFSI/ustruct.cpp @@ -28,6 +28,24 @@ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ +/** + * @file ustruct.cpp + * @brief Structural mechanics implementation based on the following reference: + * + * Ju Liu, Alison L. Marsden, + * A unified continuum and variational multiscale formulation for fluids, solids, and fluid–structure interaction, + * Computer Methods in Applied Mechanics and Engineering, + * Volume 337, + * 2018, + * Pages 549-597, + * ISSN 0045-7825, + * https://doi.org/10.1016/j.cma.2018.03.045. + * + * This paper describes a unified framework for fluid, solids, and FSI. The code + * in this file is based on the solid mechanics portion of the paper, which can be + * found in Section 4. + */ + #include "ustruct.h" #include "all_fun.h" diff --git a/tests/unitTests/test.cpp b/tests/unitTests/test.cpp index 83945553..661a4680 100644 --- a/tests/unitTests/test.cpp +++ b/tests/unitTests/test.cpp @@ -37,7 +37,11 @@ using namespace std; // --------------------------- Neo-Hookean Material --------------------------- // ---------------------------------------------------------------------------- -// Test fixture class for Neo-Hookean material model +/** + * @brief Test fixture class for the Neo-Hookean material model. + * + * This class sets up the necessary parameters and objects for testing the Neo-Hookean material model. + */ class NeoHookeanTest : public ::testing::Test { protected: // Variables common across tests @@ -70,8 +74,24 @@ class NeoHookeanTest : public ::testing::Test { } }; +// ------------------------------ STRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for STRUCT Neo-Hookean material model. + * + * This class sets up the necessary parameters and objects for testing the STRUCT Neo-Hookean material model. + */ +class STRUCT_NeoHookeanTest : public NeoHookeanTest { +protected: + void SetUp() override { + NeoHookeanTest::SetUp(); + + // Use struct + TestNH->ustruct = false; + } +}; + // Test PK2 stress zero for F = I -TEST_F(NeoHookeanTest, TestPK2StressIdentityF) { +TEST_F(STRUCT_NeoHookeanTest, TestPK2StressIdentityF) { //verbose = true; // Show values of S and S_ref // Check identity F produces zero PK2 stress @@ -83,7 +103,7 @@ TEST_F(NeoHookeanTest, TestPK2StressIdentityF) { } // Test PK2 stress with finite difference for random F -TEST_F(NeoHookeanTest, TestPK2StressRandomF) { +TEST_F(STRUCT_NeoHookeanTest, TestPK2StressRandomF) { //verbose = true; // Show values of S, dE, SdE and dPsi // Compute reference PK2 stress with finite difference for random F @@ -96,7 +116,7 @@ TEST_F(NeoHookeanTest, TestPK2StressRandomF) { } // Test PK2 stress consistent with strain energy for random F -TEST_F(NeoHookeanTest, TestPK2StressConsistentRandomF) { +TEST_F(STRUCT_NeoHookeanTest, TestPK2StressConsistentRandomF) { //verbose = true; // Show values of S, dE, SdE and dPsi // Check random F produces consistent PK2 stress @@ -105,7 +125,66 @@ TEST_F(NeoHookeanTest, TestPK2StressConsistentRandomF) { } // Test material elasticity consistent with PK2 stress -TEST_F(NeoHookeanTest, TestMaterialElasticityConsistentRandomF) { +TEST_F(STRUCT_NeoHookeanTest, TestMaterialElasticityConsistentRandomF) { + //verbose = true; // Show values of CC, dE, CCdE and dS + + // Check with random F + create_random_F(F); + TestNH->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// ------------------------------ USTRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for USTRUCT Neo-Hookean material model. + * + * This class sets up the necessary parameters and objects for testing the USTRUCT Neo-Hookean material model. + */ +class USTRUCT_NeoHookeanTest : public NeoHookeanTest { +protected: + void SetUp() override { + NeoHookeanTest::SetUp(); + + // Use ustruct + TestNH->ustruct = true; + } +}; + +// Test PK2 stress zero for F = I +TEST_F(USTRUCT_NeoHookeanTest, TestPK2StressIdentityF) { + //verbose = true; // Show values of S and S_ref + + // Check identity F produces zero PK2 stress + double F[3][3] = {{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}; + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestNH->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress with finite difference for random F +TEST_F(USTRUCT_NeoHookeanTest, TestPK2StressRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Compute reference PK2 stress with finite difference for random F + create_random_F(F); + double S_ref[3][3]; // PK2 stress + TestNH->calcPK2StressFiniteDifference(F, delta, S_ref); + + // Check PK2 stress against reference value + TestNH->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress consistent with strain energy for random F +TEST_F(USTRUCT_NeoHookeanTest, TestPK2StressConsistentRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Check random F produces consistent PK2 stress + create_random_F(F); + TestNH->testPK2StressConsistentWithStrainEnergy(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// Test material elasticity consistent with PK2 stress +TEST_F(USTRUCT_NeoHookeanTest, TestMaterialElasticityConsistentRandomF) { //verbose = true; // Show values of CC, dE, CCdE and dS // Check with random F @@ -118,7 +197,11 @@ TEST_F(NeoHookeanTest, TestMaterialElasticityConsistentRandomF) { // --------------------------- Mooney-Rivlin Material ------------------------- // ---------------------------------------------------------------------------- -// Test fixture class for Mooney-Rivlin material model +/** + * @brief Test fixture class for the Mooney-Rivlin material model. + * + * This class sets up the necessary parameters and objects for testing the Mooney-Rivlin material model. + */ class MooneyRivlinTest : public ::testing::Test { protected: // Variables common across tests @@ -153,8 +236,24 @@ class MooneyRivlinTest : public ::testing::Test { } }; +// ------------------------------ STRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for STRUCT Mooney-Rivlin material model. + * + * This class sets up the necessary parameters and objects for testing the STRUCT Mooney-Rivlin material model. + */ +class STRUCT_MooneyRivlinTest : public MooneyRivlinTest { +protected: + void SetUp() override { + MooneyRivlinTest::SetUp(); + + // Use struct + TestMR->ustruct = false; + } +}; + // Test PK2 stress zero for F = I -TEST_F(MooneyRivlinTest, TestPK2StressIdentityF) { +TEST_F(STRUCT_MooneyRivlinTest, TestPK2StressIdentityF) { //verbose = true; // Show values of S and S_ref // Check identity F produces zero PK2 stress @@ -166,7 +265,7 @@ TEST_F(MooneyRivlinTest, TestPK2StressIdentityF) { } // Test PK2 stress with finite difference for random F -TEST_F(MooneyRivlinTest, TestPK2StressRandomF) { +TEST_F(STRUCT_MooneyRivlinTest, TestPK2StressRandomF) { //verbose = true; // Show values of S, dE, SdE and dPsi // Compute reference PK2 stress with finite difference for random F @@ -179,7 +278,7 @@ TEST_F(MooneyRivlinTest, TestPK2StressRandomF) { } // Test PK2 stress consistent with strain energy for random F -TEST_F(MooneyRivlinTest, TestPK2StressConsistentRandomF) { +TEST_F(STRUCT_MooneyRivlinTest, TestPK2StressConsistentRandomF) { //verbose = true; // Show values of S, dE, SdE and dPsi // Check random F produces consistent PK2 stress @@ -188,7 +287,66 @@ TEST_F(MooneyRivlinTest, TestPK2StressConsistentRandomF) { } // Test material elasticity consistent with PK2 stress -TEST_F(MooneyRivlinTest, TestMaterialElasticityConsistentRandomF) { +TEST_F(STRUCT_MooneyRivlinTest, TestMaterialElasticityConsistentRandomF) { + //verbose = true; // Show values of CC, dE, CCdE and dS + + // Check with random F + create_random_F(F); + TestMR->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// ------------------------------ USTRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for USTRUCT Mooney-Rivlin material model. + * + * This class sets up the necessary parameters and objects for testing the USTRUCT Mooney-Rivlin material model. + */ +class USTRUCT_MooneyRivlinTest : public MooneyRivlinTest { +protected: + void SetUp() override { + MooneyRivlinTest::SetUp(); + + // Use ustruct + TestMR->ustruct = true; + } +}; + +// Test PK2 stress zero for F = I +TEST_F(USTRUCT_MooneyRivlinTest, TestPK2StressIdentityF) { + //verbose = true; // Show values of S and S_ref + + // Check identity F produces zero PK2 stress + double F[3][3] = {1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0}; + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestMR->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress with finite difference for random F +TEST_F(USTRUCT_MooneyRivlinTest, TestPK2StressRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Compute reference PK2 stress with finite difference for random F + create_random_F(F); + double S_ref[3][3]; // PK2 stress + TestMR->calcPK2StressFiniteDifference(F, delta, S_ref); + + // Check PK2 stress against reference value + TestMR->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress consistent with strain energy for random F +TEST_F(USTRUCT_MooneyRivlinTest, TestPK2StressConsistentRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Check random F produces consistent PK2 stress + create_random_F(F); + TestMR->testPK2StressConsistentWithStrainEnergy(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// Test material elasticity consistent with PK2 stress +TEST_F(USTRUCT_MooneyRivlinTest, TestMaterialElasticityConsistentRandomF) { //verbose = true; // Show values of CC, dE, CCdE and dS // Check with random F @@ -206,7 +364,11 @@ TEST_F(MooneyRivlinTest, TestMaterialElasticityConsistentRandomF) { // ----------------------- Holzapfel-Ogden Material --------------------------- // ---------------------------------------------------------------------------- -// Test fixture class for Holzapfel-Ogden material model +/** + * @brief Test fixture class for the Holzapfel-Ogden material model. + * + * This class sets up the necessary parameters and objects for testing the Holzapfel-Ogden material model. +*/ class HolzapfelOgdenTest : public ::testing::Test { protected: // Variables common across tests @@ -233,7 +395,6 @@ class HolzapfelOgdenTest : public ::testing::Test { params.b_f = 16.026; // Pa params.b_s = 11.12; // Pa params.b_fs = 11.436; // Pa - params.kappa = 1.0e6; // Pa params.k = 100000.0; // Pa // Set random values for f between 0 and 1 and normalize @@ -283,8 +444,24 @@ class HolzapfelOgdenTest : public ::testing::Test { } }; +// ------------------------------ STRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for STRUCT Holzapfel-Ogden material model. + * + * This class sets up the necessary parameters and objects for testing the STRUCT Holzapfel-Ogden material model. + */ +class STRUCT_HolzapfelOgdenTest : public HolzapfelOgdenTest { +protected: + void SetUp() override { + HolzapfelOgdenTest::SetUp(); + + // Use struct + TestHO->ustruct = false; + } +}; + // Test PK2 stress zero for F = I -TEST_F(HolzapfelOgdenTest, TestPK2StressIdentityF) { +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressIdentityF) { //verbose = true; // Show values of S and S_ref // Check identity F produces zero PK2 stress @@ -296,7 +473,7 @@ TEST_F(HolzapfelOgdenTest, TestPK2StressIdentityF) { } // Test PK2 stress for triaxial stretch -TEST_F(HolzapfelOgdenTest, TestPK2StressTriaxialStretch) { +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressTriaxialStretch) { //verbose = true; // Show values of S and S_ref // Check triaxial stretch produces PK2 stress consistent with svFSI @@ -313,7 +490,7 @@ TEST_F(HolzapfelOgdenTest, TestPK2StressTriaxialStretch) { } // Test PK2 stress for triaxial compression -TEST_F(HolzapfelOgdenTest, TestPK2StressTriaxialCompression) { +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressTriaxialCompression) { //verbose = true; // Show values of S and S_ref // Check triaxial compression produces PK2 stress consistent with svFSI @@ -330,7 +507,7 @@ TEST_F(HolzapfelOgdenTest, TestPK2StressTriaxialCompression) { } // Test PK2 stress with finite difference for random F -TEST_F(HolzapfelOgdenTest, TestPK2StressRandomF) { +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressRandomF) { //verbose = true; // Show values of S, dE, SdE and dPsi // Compute reference PK2 stress with finite difference for F = I + random perturbations @@ -343,7 +520,7 @@ TEST_F(HolzapfelOgdenTest, TestPK2StressRandomF) { } // Test PK2 stress consistent with strain energy for random F -TEST_F(HolzapfelOgdenTest, TestPK2StressConsistentRandomF) { +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressConsistentRandomF) { //verbose = true; // Show values of S, dE, SdE and dPsi // Check F = I + random perturbations produces consistent PK2 stress @@ -352,7 +529,7 @@ TEST_F(HolzapfelOgdenTest, TestPK2StressConsistentRandomF) { } // Test material elasticity consistent with PK2 stress -TEST_F(HolzapfelOgdenTest, TestMaterialElasticityConsistentRandomF) { +TEST_F(STRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistentRandomF) { //verbose = true; // Show values of CC, dE, CCdE and dS // Check F = I + random perturbations produces consistent PK2 stress @@ -360,4 +537,620 @@ TEST_F(HolzapfelOgdenTest, TestMaterialElasticityConsistentRandomF) { TestHO->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); } +// ------------------------------ USTRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for USTRUCT Holzapfel-Ogden material model. + * + * This class sets up the necessary parameters and objects for testing the USTRUCT Holzapfel-Ogden material model. + */ +class USTRUCT_HolzapfelOgdenTest : public HolzapfelOgdenTest { +protected: + void SetUp() override { + HolzapfelOgdenTest::SetUp(); + + // Use ustruct + TestHO->ustruct = true; + } +}; + +// Test PK2 stress zero for F = I +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressIdentityF) { + //verbose = true; // Show values of S and S_ref + + // Check identity F produces zero PK2 stress + double F[3][3] = {{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}; + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestHO->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress for triaxial stretch +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressTriaxialStretch) { + //verbose = true; // Show values of S and S_ref + + // Check triaxial stretch produces PK2 stress consistent with svFSI + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.3}}; + + // Compute reference PK2 stress with finite difference for triaxial stretch + double S_ref[3][3]; // PK2 stress + TestHO->calcPK2StressFiniteDifference(F, delta, S_ref); + + // Check PK2 stress against reference value + TestHO->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress for triaxial compression +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressTriaxialCompression) { + //verbose = true; // Show values of S and S_ref + + // Check triaxial compression produces PK2 stress consistent with svFSI + double F[3][3] = {{0.9, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 0.7}}; + + // Compute reference PK2 stress with finite difference for triaxial compression + double S_ref[3][3]; // PK2 stress + TestHO->calcPK2StressFiniteDifference(F, delta, S_ref); + + // Check PK2 stress against reference value + TestHO->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress with finite difference for random F +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Compute reference PK2 stress with finite difference for F = I + random perturbations + create_random_perturbed_identity_F(F, 0.5); + double S_ref[3][3]; // PK2 stress + TestHO->calcPK2StressFiniteDifference(F, delta, S_ref); + + // Check PK2 stress against reference value + TestHO->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress consistent with strain energy for random F +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressConsistentRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Check F = I + random perturbations produces consistent PK2 stress + create_random_perturbed_identity_F(F, 0.5); + TestHO->testPK2StressConsistentWithStrainEnergy(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// Test material elasticity consistent with PK2 stress +TEST_F(USTRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistentRandomF) { + //verbose = true; // Show values of CC, dE, CCdE and dS + + // Check F = I + random perturbations produces consistent PK2 stress + create_random_perturbed_identity_F(F, 0.5); + TestHO->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + + + + +// ---------------------------------------------------------------------------- +// ---------------------------------------------------------------------------- +// ----------------------- Volumetric penalty models ------------------------- +// ---------------------------------------------------------------------------- +// ---------------------------------------------------------------------------- + +// ---------------------------------------------------------------------------- +// --------------------------- Quadratic Volumetric Penalty Model ------------------------ +// ---------------------------------------------------------------------------- +/** + * @brief Test fixture class for the Quadratic Volumetric penalty model. + * + * This class sets up the necessary parameters and objects for testing the Quadratic Volumetric penalty model. + */ +class QuadraticVolumetricPenaltyTest : public ::testing::Test { +protected: + // Variables common across tests + VolumetricPenaltyParams params; + double F[3][3] = {}; // Deformation gradient + 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-9; // absolute tolerance for comparing values + double delta = 1e-7; // perturbation scaling factor + bool verbose = false; // Show values of S, dE, SdE and dPsi + + // Add the test object + TestQuadraticVolumetricPenalty* TestQVP; + + // Setup method to initialize variables before each test + void SetUp() override { + + // Set random values for the Quadratic penalty parameters between 1000 and 10000 + params.kappa = getRandomDouble(1000.0, 10000.0); + + // Initialize the test object + TestQVP = new TestQuadraticVolumetricPenalty(params); + } + + // TearDown method to clean up after each test, if needed + void TearDown() override { + // Clean up the test object + delete TestQVP; + TestQVP = nullptr; + } +}; + +// ------------------------------ STRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for STRUCT Quadratic penalty model. + * + * This class sets up the necessary parameters and objects for testing the STRUCT Quadratic penalty model. + */ +class STRUCT_QuadraticVolumetricPenaltyTest : public QuadraticVolumetricPenaltyTest { +protected: + void SetUp() override { + QuadraticVolumetricPenaltyTest::SetUp(); + + // Use struct + //TestQVP->ustruct = false; + } +}; + +// Test PK2 stress zero for F = I +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressIdentityF) { + //verbose = true; // Show values of S and S_ref + + // Check identity F produces zero PK2 stress + double F[3][3] = {{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}; + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestQVP->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress zero for prescribed isochoric deformation +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation) { + //verbose = true; // Show values of S and S_ref + + // Check isochoric deformation produces zero PK2 stress + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.0/(1.1*1.2)}}; + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestQVP->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress zero for random isochoric deformation +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) { + //verbose = true; // Show values of S and S_ref + + // Create random deformation gradient tensor + create_random_F(F); + + // Make det(F) = 1 + double J = mat_fun_carray::mat_det(F); + double J13 = pow(J, 1.0/3.0); + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + F[i][j] /= J13; + } + } + + // Check random isochoric deformation produces zero PK2 stress + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestQVP->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress with finite difference for random F +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Compute reference PK2 stress with finite difference for random F + create_random_F(F); + double S_ref[3][3]; // PK2 stress + TestQVP->calcPK2StressFiniteDifference(F, delta, S_ref); + + // Check PK2 stress against reference value + TestQVP->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress consistent with strain energy for random F +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressConsistentRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Check random F produces consistent PK2 stress + create_random_F(F); + TestQVP->testPK2StressConsistentWithStrainEnergy(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// Test material elasticity consistent with PK2 stress +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) { + //verbose = true; // Show values of CC, dE, CCdE and dS + + // Check with random F + create_random_F(F); + TestQVP->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// ------------------------------ USTRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for USTRUCT Quadratic penalty model. + * + * This class sets up the necessary parameters and objects for testing the USTRUCT Quadratic penalty model. + */ +class USTRUCT_QuadraticVolumetricPenaltyTest : public QuadraticVolumetricPenaltyTest { +protected: + void SetUp() override { + QuadraticVolumetricPenaltyTest::SetUp(); + + // Use ustruct + //TestQVP->ustruct = true; + } +}; + +// Test rho, beta, drho/dp and dbeta/dp for random pressure +TEST_F(USTRUCT_QuadraticVolumetricPenaltyTest, TestRhoBeta) { + //verbose = true; // Show values of rho, beta, drho/dp and dbeta/dp and reference values + + // Generate random pressure p between 100 and 1000 + double p = getRandomDouble(100.0, 1000.0); + + // Other parameters + double rho0 = 1000.0; // Reference density + double kappa = TestQVP->params.kappa; // Volumetric penalty parameter + + // Compute reference values for rho, beta, drho/dp and dbeta/dp + // See ustruct paper (https://doi.org/10.1016/j.cma.2018.03.045) Section 2.4 + double rho_ref = rho0 / (1.0 - p/kappa); + double beta_ref = 1.0 / (kappa - p); + double drhodp_ref = -rho0 / pow(1.0 - p/kappa, 2) * (-1.0/kappa); // Derivative of rho with respect to p + double dbetadp_ref = -1.0 / pow(kappa - p, 2) * (-1.0); // Derivative of beta with respect to p + + // Check rho, beta, drho/dp and dbeta/dp against reference values + TestQVP->testRhoBetaAgainstReference(p, rho0, rho_ref, beta_ref, drhodp_ref, dbetadp_ref, rel_tol, abs_tol, verbose); +} + +// ---------------------------------------------------------------------------- +// --------------------------- Simo-Taylor91 Volumetric Penalty Model --------- +// ---------------------------------------------------------------------------- + +/** + * @brief Test fixture class for the Simo-Taylor91 Volumetric penalty model. + * + * This class sets up the necessary parameters and objects for testing the Simo-Taylor91 Volumetric penalty model. + */ +class SimoTaylor91VolumetricPenaltyTest : public ::testing::Test { +protected: + // Variables common across tests + VolumetricPenaltyParams params; + double F[3][3] = {}; // Deformation gradient + 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-9; // absolute tolerance for comparing values + double delta = 1e-7; // perturbation scaling factor + bool verbose = false; // Show values of S, dE, SdE and dPsi + + // Add the test object + TestSimoTaylor91VolumetricPenalty* TestST91; + + // Setup method to initialize variables before each test + void SetUp() override { + + // Set random values for the Simo-Taylor91 penalty parameters between 1000 and 10000 + params.kappa = getRandomDouble(1000.0, 10000.0); + + // Initialize the test object + TestST91 = new TestSimoTaylor91VolumetricPenalty(params); + } + + // TearDown method to clean up after each test, if needed + void TearDown() override { + // Clean up the test object + delete TestST91; + TestST91 = nullptr; + } +}; + +// ------------------------------ STRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for STRUCT Simo-Taylor91 penalty model. + * + * This class sets up the necessary parameters and objects for testing the STRUCT Simo-Taylor91 penalty model. + */ +class STRUCT_SimoTaylor91VolumetricPenaltyTest : public SimoTaylor91VolumetricPenaltyTest { +protected: + void SetUp() override { + SimoTaylor91VolumetricPenaltyTest::SetUp(); + + // Use struct + //TestST91->ustruct = false; + } +}; + +// Test PK2 stress zero for F = I +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressIdentityF) { + //verbose = true; // Show values of S and S_ref + + // Check identity F produces zero PK2 stress + double F[3][3] = {{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}; + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestST91->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress zero for prescribed isochoric deformation +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation) { + //verbose = true; // Show values of S and S_ref + + // Check isochoric deformation produces zero PK2 stress + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.0/(1.1*1.2)}}; + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestST91->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress zero for random isochoric deformation +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) { + //verbose = true; // Show values of S and S_ref + + // Create random deformation gradient tensor + create_random_F(F); + + // Make det(F) = 1 + double J = mat_fun_carray::mat_det(F); + double J13 = pow(J, 1.0/3.0); + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + F[i][j] /= J13; + } + } + + // Check random isochoric deformation produces zero PK2 stress + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestST91->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress with finite difference for random F +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Compute reference PK2 stress with finite difference for random F + create_random_F(F); + double S_ref[3][3]; // PK2 stress + TestST91->calcPK2StressFiniteDifference(F, delta, S_ref); + + // Check PK2 stress against reference value + TestST91->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress consistent with strain energy for random F +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressConsistentRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Check random F produces consistent PK2 stress + create_random_F(F); + TestST91->testPK2StressConsistentWithStrainEnergy(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// Test material elasticity consistent with PK2 stress +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) { + //verbose = true; // Show values of CC, dE, CCdE and dS + + // Check with random F + create_random_F(F); + TestST91->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// ------------------------------ USTRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for USTRUCT Simo-Taylor91 penalty model. + * + * This class sets up the necessary parameters and objects for testing the USTRUCT Simo-Taylor91 penalty model. + */ +class USTRUCT_SimoTaylor91VolumetricPenaltyTest : public SimoTaylor91VolumetricPenaltyTest { +protected: + void SetUp() override { + SimoTaylor91VolumetricPenaltyTest::SetUp(); + + // Use ustruct + //TestST91->ustruct = true; + } +}; + +// Test rho and beta values for random p +TEST_F(USTRUCT_SimoTaylor91VolumetricPenaltyTest, TestRhoBeta) { + //verbose = true; // Show values of rho, beta and rho_ref + + // Generate random pressure p between 100 and 1000 + double p = getRandomDouble(100.0, 1000.0); + + // Other parameters + double rho0 = 1000.0; // Reference density + double kappa = TestST91->params.kappa; // Volumetric penalty parameter + + // Compute reference values for rho, beta, drho/dp and dbeta/dp + // See ustruct paper (https://doi.org/10.1016/j.cma.2018.03.045) Section 2.4 + double rho_ref = rho0 / kappa * (sqrt(p*p + kappa*kappa) + p); + double beta_ref = 1.0 / sqrt(p*p + kappa*kappa); + double drhodp_ref = rho0 / kappa * (1.0/2.0 * (1.0/sqrt(p*p + kappa*kappa)) * 2.0*p + 1.0); // Derivative of rho with respect to p + double dbetadp_ref = -1.0/2.0 * 1.0/pow(p*p + kappa*kappa, 3.0/2.0) * 2.0*p; // Derivative of beta with respect to p + + // Check rho, beta, drho/dp and dbeta/dp against reference values + TestST91->testRhoBetaAgainstReference(p, rho0, rho_ref, beta_ref, drhodp_ref, dbetadp_ref, rel_tol, abs_tol, verbose); +} + +// ---------------------------------------------------------------------------- +// --------------------------- Miehe94 Volumetric Penalty Model --------------- +// ---------------------------------------------------------------------------- +/** + * @brief Test fixture class for the Miehe94 Volumetric penalty model. + * + * This class sets up the necessary parameters and objects for testing the Miehe94 Volumetric penalty model. + */ +class Miehe94VolumetricPenaltyTest : public ::testing::Test { +protected: + // Variables common across tests + VolumetricPenaltyParams params; + double F[3][3] = {}; // Deformation gradient + 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-9; // absolute tolerance for comparing values + double delta = 1e-7; // perturbation scaling factor + bool verbose = false; // Show values of S, dE, SdE and dPsi + + // Add the test object + TestMiehe94VolumetricPenalty* TestM94; + + // Setup method to initialize variables before each test + void SetUp() override { + + // Set random values for the Miehe94 penalty parameters between 1000 and 10000 + params.kappa = getRandomDouble(1000.0, 10000.0); + + // Initialize the test object + TestM94 = new TestMiehe94VolumetricPenalty(params); + } + + // TearDown method to clean up after each test, if needed + void TearDown() override { + // Clean up the test object + delete TestM94; + TestM94 = nullptr; + } +}; + +// ------------------------------ STRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for STRUCT Miehe94 penalty model. + * + * This class sets up the necessary parameters and objects for testing the STRUCT Miehe94 penalty model. + */ +class STRUCT_Miehe94VolumetricPenaltyTest : public Miehe94VolumetricPenaltyTest { +protected: + void SetUp() override { + Miehe94VolumetricPenaltyTest::SetUp(); + + // Use struct + //TestM94->ustruct = false; + } +}; + +// Test PK2 stress zero for F = I +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressIdentityF) { + //verbose = true; // Show values of S and S_ref + + // Check identity F produces zero PK2 stress + double F[3][3] = {{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}; + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestM94->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress zero for prescribed isochoric deformation +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation) { + //verbose = true; // Show values of S and S_ref + + // Check isochoric deformation produces zero PK2 stress + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.0/(1.1*1.2)}}; + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestM94->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress zero for random isochoric deformation +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) { + //verbose = true; // Show values of S and S_ref + + // Create random deformation gradient tensor + create_random_F(F); + + // Make det(F) = 1 + double J = mat_fun_carray::mat_det(F); + double J13 = pow(J, 1.0/3.0); + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + F[i][j] /= J13; + } + } + + // Check random isochoric deformation produces zero PK2 stress + double S_ref[3][3] = {}; // PK2 stress initialized to zero + TestM94->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress with finite difference for random F +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Compute reference PK2 stress with finite difference for random F + create_random_F(F); + double S_ref[3][3]; // PK2 stress + TestM94->calcPK2StressFiniteDifference(F, delta, S_ref); + + // Check PK2 stress against reference value + TestM94->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test PK2 stress consistent with strain energy for random F +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressConsistentRandomF) { + //verbose = true; // Show values of S, dE, SdE and dPsi + + // Check random F produces consistent PK2 stress + create_random_F(F); + TestM94->testPK2StressConsistentWithStrainEnergy(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// Test material elasticity consistent with PK2 stress +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) { + //verbose = true; // Show values of CC, dE, CCdE and dS + + // Check with random F + create_random_F(F); + TestM94->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); +} + +// ------------------------------ USTRUCT TESTS -------------------------------- +/** + * @brief Test fixture class for USTRUCT Miehe94 penalty model. + * + * This class sets up the necessary parameters and objects for testing the USTRUCT Miehe94 penalty model. + */ +class USTRUCT_Miehe94VolumetricPenaltyTest : public Miehe94VolumetricPenaltyTest { +protected: + void SetUp() override { + Miehe94VolumetricPenaltyTest::SetUp(); + + // Use ustruct + //TestM94->ustruct = true; + } +}; + +// Test rho and beta values for random p +TEST_F(USTRUCT_Miehe94VolumetricPenaltyTest, TestRhoBeta) { + //verbose = true; // Show values of rho, beta and rho_ref + + // Generate random pressure p between 100 and 1000 + double p = getRandomDouble(100.0, 1000.0); + + // Other parameters + double rho0 = 1000.0; // Reference density + double kappa = TestM94->params.kappa; // Volumetric penalty parameter + + // Compute reference values for rho, beta, drho/dp and dbeta/dp + // See ustruct paper (https://doi.org/10.1016/j.cma.2018.03.045) Section 2.4 + double rho_ref = rho0 * (1.0 + p/kappa); + double beta_ref = 1.0 / (p + kappa); + double drhodp_ref = rho0 / kappa; // Derivative of rho with respect to p + double dbetadp_ref = -1.0 / pow(p + kappa, 2); // Derivative of beta with respect to p + + // Check rho, beta, drho/dp and dbeta/dp against reference values + TestM94->testRhoBetaAgainstReference(p, rho0, rho_ref, beta_ref, drhodp_ref, dbetadp_ref, rel_tol, abs_tol, verbose); +} + +// ---------------------------------------------------------------------------- +// ---------------------------------------------------------------------------- + diff --git a/tests/unitTests/test.h b/tests/unitTests/test.h index 2c409a9e..0f95e7b2 100644 --- a/tests/unitTests/test.h +++ b/tests/unitTests/test.h @@ -28,6 +28,12 @@ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ +// -------------------------------------------------------------- +// To run the tests in test.cpp +// 1. Navigate to /build/svFSI-build/Source/svFSI +// 2. Run `make` to build the tests +// 3. Run `ctest --verbose` to run the tests + // -------------------------------------------------------------- // To add a new material model to test: @@ -99,7 +105,6 @@ class HolzapfelOgdenParams : public MatParams { double b_s; double a_fs; double b_fs; - double kappa; double f[3]; // Fiber direction double s[3]; // Sheet direction @@ -108,7 +113,7 @@ class HolzapfelOgdenParams : public MatParams { bool full_anisotropic_invariants; // Flag to use full anisotropic invariants in strain energy density function // Default constructor - HolzapfelOgdenParams() : a(0.0), b(0.0), a_f(0.0), b_f(0.0), a_s(0.0), b_s(0.0), a_fs(0.0), b_fs(0.0), kappa(0.0), k(0.0) { + HolzapfelOgdenParams() : a(0.0), b(0.0), a_f(0.0), b_f(0.0), a_s(0.0), b_s(0.0), a_fs(0.0), b_fs(0.0), k(0.0) { for (int i = 0; i < 3; i++) { f[i] = 0.0; s[i] = 0.0; @@ -116,7 +121,7 @@ class HolzapfelOgdenParams : public MatParams { } // Constructor with parameters - HolzapfelOgdenParams(double a, double b, double a_f, double b_f, double a_s, double b_s, double a_fs, double b_fs, double kappa, double k, double f[3], double s[3]) : a(a), b(b), a_f(a_f), b_f(b_f), a_s(a_s), b_s(b_s), a_fs(a_fs), b_fs(b_fs), kappa(kappa), k(k) { + HolzapfelOgdenParams(double a, double b, double a_f, double b_f, double a_s, double b_s, double a_fs, double b_fs, double k, double f[3], double s[3]) : a(a), b(b), a_f(a_f), b_f(b_f), a_s(a_s), b_s(b_s), a_fs(a_fs), b_fs(b_fs), k(k) { for (int i = 0; i < 3; i++) { this->f[i] = f[i]; this->s[i] = s[i]; @@ -124,11 +129,28 @@ class HolzapfelOgdenParams : public MatParams { } }; +// Class to contain volumetric penalty parameters (just the penalty parameter) +class VolumetricPenaltyParams : public MatParams { +public: + double kappa; + + // Default constructor + VolumetricPenaltyParams() : kappa(0.0) {} + + // Constructor with parameters + VolumetricPenaltyParams(double kappa) : kappa(kappa) {} +}; + // -------------------------------------------------------------- // ---------------------- Helper functions ---------------------- // -------------------------------------------------------------- -// Creates an identity deformation gradient F +/** + * @brief Creates an identity deformation gradient F. + * + * @param[out] F The deformation gradient tensor to be set to the identity matrix. + * @return None, but fills F with the identity matrix. + */ template void create_identity_F(double F[N][N]) { for (int i = 0; i < N; i++) { @@ -138,7 +160,15 @@ void create_identity_F(double F[N][N]) { } } -// Inline getRandomDouble function +/** + * @brief Generates a random double value. + * + * This function generates a random double value within a specified range. + * + * @param[in] min The minimum value of the range. + * @param[in] max The maximum value of the range. + * @return A random double value between min and max. + */ 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(); @@ -148,8 +178,17 @@ inline double getRandomDouble(double min, double max) { return distribution(engine); } -// Creates a random deformation gradient F with values between min and max, -// and det(F) > 0 +/** + * @brief Creates a random deformation gradient F with values between min and max, and det(F) > 0. + * + * This function generates a random deformation gradient tensor F such that the determinant of F is greater than 0. + * + * @tparam N The size of the deformation gradient tensor (NxN). + * @param[out] F The generated deformation gradient tensor. + * @param[in] min The minimum value for the elements of the deformation gradient tensor (default is 0.1). + * @param[in] max The maximum value for the elements of the deformation gradient tensor (default is 10.0). + * @return None. + */ template void create_random_F(double F[N][N], double min=0.1, double max=10.0) { // Create a random deformation gradient with values between min and max, @@ -165,7 +204,17 @@ void create_random_F(double F[N][N], double min=0.1, double max=10.0) { } } -// Creates a deformation matrix F of random deviations from the identity matrix +/** + * @brief Creates a deformation matrix F of random deviations from the identity matrix. + * + * This function generates a random deformation gradient tensor F with values perturbed from the identity matrix, + * such that the determinant of F is greater than 0. + * + * @tparam N The size of the deformation gradient tensor (NxN). + * @param[out] F The generated deformation gradient tensor. + * @param[in] max_deviation The maximum deviation from the identity matrix elements. + * @return None. + */ template void create_random_perturbed_identity_F(double F[N][N], double max_deviation) { // Create a random deformation gradient with values perturbed from the identity matrix, @@ -181,8 +230,18 @@ void create_random_perturbed_identity_F(double F[N][N], double max_deviation) { } } -// Perturb the deformation gradient F by delta times a random number between -1 and 1 -// and store the perturbed deformation gradient in F_tilde +/** + * @brief Perturbs the deformation gradient F by delta times a random number between -1 and 1. + * + * This function perturbs the given deformation gradient tensor F by adding delta times a random number + * between -1 and 1 to each element, and stores the perturbed deformation gradient in F_tilde. + * + * @tparam N The size of the deformation gradient tensor (NxN). + * @param[in] F The original deformation gradient tensor. + * @param[in] delta The perturbation factor. + * @param[out] F_tilde The perturbed deformation gradient tensor. + * @return None. + */ template void perturb_random_F(const double F[N][N], const double delta, double F_tilde[N][N]) { @@ -196,8 +255,19 @@ void perturb_random_F(const double F[N][N], const double delta, double F_tilde[N } } -// Compute Jacobian J, right Cauchy-Green deformation tensor C, and Green-Lagrange -// strain tensor E from the deformation gradient F +/** + * @brief Computes the Jacobian J, right Cauchy-Green deformation tensor C, and Green-Lagrange strain tensor E from the deformation gradient F. + * + * This function computes the Jacobian of the deformation gradient tensor F, the right Cauchy-Green deformation tensor C, + * and the Green-Lagrange strain tensor E. + * + * @tparam N The size of the deformation gradient tensor (NxN). + * @param[in] F The deformation gradient tensor. + * @param[out] J The computed Jacobian of F. + * @param[out] C The computed right Cauchy-Green deformation tensor. + * @param[out] E The computed Green-Lagrange strain tensor. + * @return None. + */ template void calc_JCE(const double F[N][N], double &J, double C[N][N], double E[N][N]) { // Compute Jacobian of F @@ -218,20 +288,34 @@ void calc_JCE(const double F[N][N], double &J, double C[N][N], double E[N][N]) { } } -// Structure to store solid mechanics terms used to compute strain energy density -// functions. +/** + * @brief Structure to store solid mechanics terms used to compute strain energy density functions. + * + * @tparam N The size of the deformation gradient tensor (NxN). + */ template struct solidMechanicsTerms { - double J; - double C[N][N]; - double E[N][N]; - double E2[N][N]; - double C_bar[N][N]; - double I1, I2, Ib1, Ib2; + double J; /**< Jacobian of the deformation gradient tensor. */ + double C[N][N]; /**< Right Cauchy-Green deformation tensor. */ + double E[N][N]; /**< Green-Lagrange strain tensor. */ + double E2[N][N]; /**< Second-order Green-Lagrange strain tensor. */ + double C_bar[N][N];/**< Modified right Cauchy-Green deformation tensor. */ + double I1; /**< First invariant of the right Cauchy-Green deformation tensor. */ + double I2; /**< Second invariant of the right Cauchy-Green deformation tensor. */ + double Ib1; /**< First invariant of the modified right Cauchy-Green deformation tensor. */ + double Ib2; /**< Second invariant of the modified right Cauchy-Green deformation tensor. */ }; -// Function to compute the solid mechanics terms used to compute strain energy density -// functions +/** + * @brief Computes the solid mechanics terms used to compute strain energy density functions. + * + * This function computes various solid mechanics terms such as the Jacobian, right Cauchy-Green deformation tensor, + * Green-Lagrange strain tensor, and their invariants from the given deformation gradient tensor F. + * + * @tparam N The size of the deformation gradient tensor (NxN). + * @param[in] F The deformation gradient tensor. + * @return A structure containing the computed solid mechanics terms. + */ template solidMechanicsTerms calcSolidMechanicsTerms(const double F[N][N]) { solidMechanicsTerms out; @@ -348,10 +432,12 @@ class TestMaterialModel { int nFn; Array fN; double ya_g; + bool ustruct; TestMaterialModel(const consts::ConstitutiveModelType matType, const consts::ConstitutiveModelType penType) { int nsd = com_mod.nsd; - mat_fun_carray::ten_init(nsd); // initialize tensor index pointer + mat_fun_carray::ten_init(nsd); // initialize tensor index pointer for mat_fun_carray + mat_fun::ten_init(nsd); // initialize tensor index pointer for mat_fun // Set material and penalty models auto &dmn = com_mod.mockEq.mockDmn; @@ -363,6 +449,11 @@ class TestMaterialModel { fN = Array(nsd, nFn); // Fiber directions array (initialized to zeros) ya_g = 0.0; // ? + // Flag to use struct or ustruct material models + // If struct, calls get_pk2cc() and uses strain energy composed of isochoric and volumetric parts + // If ustruct, calls get_pk2cc_dev() and uses strain energy composed of isochoric part only + ustruct = false; + // Material parameters are set in each derived class } @@ -372,15 +463,96 @@ class TestMaterialModel { // Pure virtual method for computing Strain Energy virtual double computeStrainEnergy(const double F[3][3]) = 0; - // Function to get S and Dm for a given F, from the get_pk2cc function in mat_models_carray.h - void get_pk2cc(double F[3][3], double S[3][3], double Dm[6][6]) { + /** + * @brief Computes the PK2 stress tensor S and material elasticity tensor Dm for a given deformation gradient F. + * + * This function computes the PK2 stress tensor S and the material elasticity tensor Dm from the deformation gradient F. + * If `ustruct` is true, the deviatoric part of the PK2 stress tensor is returned using the `get_pk2cc_dev` function. + * + * @param[in] F The deformation gradient tensor. + * @param[out] S The computed PK2 stress tensor. + * @param[out] Dm The computed material elasticity tensor. + * @return None, but fills S and Dm with the computed values. + */ + void get_pk2cc(const double F[3][3], double S[3][3], double Dm[6][6]) { auto &dmn = com_mod.mockEq.mockDmn; - mat_models_carray::get_pk2cc<3>(com_mod, cep_mod, dmn, F, nFn, fN, ya_g, S, Dm); + if (ustruct) { + double J = 0; // Jacobian (not used in this function) + + // Cast F, S, and Dm to Array for use in get_pk2cc_dev + Array F_arr(3,3); + Array S_arr(3,3); + Array Dm_arr(6,6); + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + F_arr(i, J) = F[i][J]; + } + } + + mat_models::get_pk2cc_dev(com_mod, cep_mod, dmn, F_arr, nFn, fN, ya_g, S_arr, Dm_arr, J); + + // Copy data from S_arr and Dm_arr to S and Dm + for (int I = 0; I < 3; I++) { + for (int J = 0; J < 3; J++) { + S[I][J] = S_arr(I, J); + } + } + for (int I = 0; I < 6; I++) { + for (int J = 0; J < 6; J++) { + Dm[I][J] = Dm_arr(I, J); + } + } + + } else { + mat_models_carray::get_pk2cc<3>(com_mod, cep_mod, dmn, F, nFn, fN, ya_g, S, Dm); + } } - // Function to compute the PK2 stress tensor S(F) from the strain energy density Psi(F) - // using finite differences + /** + * @brief Computes the solid density, isothermal compressibility coefficient, and their derivatives for a given pressure. + * + * This function computes the solid density (rho), isothermal compressibility coefficient (beta), + * and their derivatives with respect to pressure (drho and dbeta) for a given pressure (p) using the g_vol_pen() function + * from mat_models.h. + * + * @param[in] p Pressure. + * @param[in] rho0 Initial solid density. + * @param[out] rho Computed solid density. + * @param[out] beta Computed isothermal compressibility coefficient. + * @param[out] drho Computed Derivative of solid density with respect to pressure. + * @param[out] dbeta Computed Derivative of beta with respect to pressure. + * @param[in] Ja Jacobian (not used in this function). + * @return None. + */ + void g_vol_pen(const double p, const double rho0, double &rho, double &beta, double &drho, double &dbeta, const double Ja) { + auto &dmn = com_mod.mockEq.mockDmn; + dmn.prop[consts::PhysicalProperyType::solid_density] = rho0; // Set initial solid density + + mat_models::g_vol_pen(com_mod, dmn, p, rho, beta, drho, dbeta, Ja); + } + + /** + * @brief Computes the PK2 stress tensor S(F) from the strain energy density Psi(F) using finite differences. + * + * Analytically, we should have S = dPsi/dE. Since we have Psi(F), we cannot directly compute S. + * Instead, we compute S = F^-1 * P, where P = dPsi/dF is computed using finite differences in each component of F. + * + * Pseudocode: + * - Compute strain energy density Psi(F) + * - For each component of F, F[i][J] + * - Perturb F[i][J] by delta to get F_tilde + * - Compute Psi(F_tilde) + * - Compute dPsi = Psi(F_tilde) - Psi(F) + * - Compute P[i][J] = dPsi / delta + * - Compute S = F^-1 * P + * + * @tparam N The size of the deformation gradient tensor (NxN). + * @param[in] F The deformation gradient tensor. + * @param[in] delta The perturbation scaling factor. + * @param[out] S The computed PK2 stress tensor. + * @return None, but fills S with the computed values. + */ template void calcPK2StressFiniteDifference(const double F[N][N], double delta, double (&S)[N][N]) { // Compute strain energy density given F @@ -414,29 +586,30 @@ class TestMaterialModel { double F_inv[N][N]; mat_fun_carray::mat_inv(F, F_inv); mat_fun_carray::mat_mul(F_inv, P, S); - - } - // Test the consistency of the PK2 stress tensor S(F) from get_pk2cc() with the strain - // energy density Psi(F) provided by the user. - // Analytically, we should have S = dPsi/dE. We are checking whether - // S:dE = dPsi, where dE and dPsi are computed using finite differences in F. - // - // ARGS: - // - F: Deformation gradient - // - n_iter: Number of random perturbations to test - // - rel_tol: Relative tolerance for comparing dPsi and S:dE - // - delta: Perturbation scaling factor - // - verbose: Show values of S, dE, SdE and dPsi - // - // Psuedocode: - // - Compute Psi(F) - // - Compute S(F) from get_pk2cc() - // - For many random dF - // - Compute dPsi = Psi(F + dF) - Psi(F) - // - Compute dE from dF - // - Check that S:dE = dPsi + /** + * @brief Tests the consistency of the PK2 stress tensor S(F) from get_pk2cc() with the strain energy density Psi(F) provided by the user. + * + * Analytically, we should have S = dPsi/dE. This function checks whether S:dE = dPsi, where dE and dPsi are computed using finite differences in F. + * + * Pseudocode: + * - Compute Psi(F) + * - Compute S(F) from get_pk2cc() + * - For many random dF + * - Compute dPsi = Psi(F + dF) - Psi(F) + * - Compute dE from dF + * - Check that S:dE = dPsi + * + * @param[in] F Deformation gradient. + * @param[in] n_iter Number of random perturbations to test. + * @param[in] rel_tol Relative tolerance for comparing dPsi and S:dE. + * @param[in] abs_tol Absolute tolerance for comparing dPsi and S:dE. + * @param[in] delta Perturbation scaling factor. + * @param[in] verbose Show values of S, dE, SdE, and dPsi if true. + * @return None. + * + */ void testPK2StressConsistentWithStrainEnergy(double F[3][3], int n_iter, double rel_tol, double abs_tol, double delta, bool verbose = false) { // Compute E from F double J, C[3][3], E[3][3]; @@ -473,7 +646,6 @@ class TestMaterialModel { } } - // Compute S:dE double SdE = mat_fun_carray::mat_ddot<3>(S, dE); @@ -516,27 +688,29 @@ class TestMaterialModel { } } - // Test the consistency of the material elasticity tensor CC(F) from get_pk2cc() with the - // PK2 stress tensor S(F) from get_pk2cc() - // Analytically, we should have CC:dE = dS. We are checking whether - // CC:dE = dS, where dE and dS are computed using finite differences in F. - // - // ARGS: - // - F: Deformation gradient - // - n_iter: Number of random perturbations to test - // - rel_tol: Relative tolerance for comparing dS and CC:dE - // - delta: Perturbation scaling factor - // - verbose: Show values of CC, dE, CCdE and dS - // - // Psuedocode: - // - Compute S(F) and CC(F) from get_pk2cc() - // - For many random dF - // - Compute S(F + dF) from get_pk2cc() - // - Compute dS = S(F + dF) - S(F) - // - Compute dE from dF - // - Check that CC:dE = dS + /** + * @brief Tests the consistency of the material elasticity tensor CC(F) from get_pk2cc() with the PK2 stress tensor S(F) from get_pk2cc(). + * + * Analytically, we should have CC:dE = dS. This function checks whether CC:dE = dS, where dE and dS are computed using finite differences in F. + * + * Pseudocode: + * - Compute S(F) and CC(F) from get_pk2cc() + * - For many random dF + * - Compute S(F + dF) from get_pk2cc() + * - Compute dS = S(F + dF) - S(F) + * - Compute dE from dF + * - Check that CC:dE = dS + * + * @param[in] F Deformation gradient. + * @param[in] n_iter Number of random perturbations to test. + * @param[in] rel_tol Relative tolerance for comparing dS and CC:dE. + * @param[in] abs_tol Absolute tolerance for comparing dS and CC:dE. + * @param[in] delta Perturbation scaling factor. + * @param[in] verbose Show values of CC, dE, CCdE, and dS if true. + * @return None. + */ void testMaterialElasticityConsistentWithPK2Stress(double F[3][3], int n_iter, double rel_tol, double abs_tol, double delta, bool verbose = false) { - + // Compute E from F double J, C[3][3], E[3][3]; calc_JCE(F, J, C, E); @@ -550,12 +724,12 @@ class TestMaterialModel { // Calculate CC from Dm double CC[3][3][3][3]; mat_models_carray::voigt_to_cc_carray<3>(Dm, CC); - + // ------- Ancillary test --------- // Calculate Dm_check from CC double Dm_check[6][6]; mat_models_carray::cc_to_voigt_carray<3>(CC, Dm_check); - + // Check that Dm_check = Dm, for sanity for (int i = 0; i < 6; i++) { for (int j = 0; j < 6; j++) { @@ -563,7 +737,7 @@ class TestMaterialModel { } } // ------------------------------- - + // Generate many random dF and check that CC:dE = dS // CC was obtained from get_pk2cc(), and dS = S(F + dF) - S(F), // where S is also obtained from get_pk2cc() @@ -575,20 +749,20 @@ class TestMaterialModel { for (int i = 0; i < n_iter; i++) { // Perturb the deformation gradient perturb_random_F<3>(F, delta, F_tilde); - + // Compute perturbed E_tilde from F_tilde calc_JCE(F_tilde, J_tilde, C_tilde, E_tilde); - + // Compute dE for (int i = 0; i < 3; i++) { for (int J = 0; J < 3; J++) { dE[i][J] = E_tilde[i][J] - E[i][J]; } } - + // Compute perturbed S_tilde with perturbed deformation gradient get_pk2cc(F_tilde, S_tilde, Dm_tilde); - + // Compute dS double dS[3][3]; for (int i = 0; i < 3; i++) { @@ -596,7 +770,7 @@ class TestMaterialModel { dS[i][j] = S_tilde[i][j] - S[i][j]; } } - + // Check that CC_ijkl dE_kl = dS_ij mat_fun_carray::ten_mat_ddot<3>(CC, dE, CCdE); for (int i = 0; i < 3; i++) { @@ -604,13 +778,13 @@ class TestMaterialModel { EXPECT_NEAR(CCdE[i][j], dS[i][j], fmax(abs_tol, rel_tol * fabs(dS[i][j]))); } } - + // Print results if verbose if (verbose) { std::cout << "Iteration " << i << ":" << std::endl; - + printMaterialParameters(); - + std::cout << "F =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -618,7 +792,7 @@ class TestMaterialModel { } std::cout << std::endl; } - + std::cout << "CC =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -631,7 +805,7 @@ class TestMaterialModel { } std::cout << std::endl; } - + std::cout << "dE =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -639,7 +813,7 @@ class TestMaterialModel { } std::cout << std::endl; } - + std::cout << "dS =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -647,7 +821,7 @@ class TestMaterialModel { } std::cout << std::endl; } - + std::cout << "CCdE =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -660,27 +834,35 @@ class TestMaterialModel { } } - // Function to compare PK2 stress with reference solution - // ARGS: - // - F: Deformation gradient - // - S_ref: Reference solution for PK2 stress - // - rel_tol: Relative tolerance for comparing S with S_ref + /** + * @brief Compares the PK2 stress tensor S(F) with a reference solution. + * + * This function computes the PK2 stress tensor S(F) from the deformation gradient F using get_pk2cc() + * and compares it with a reference solution S_ref. The comparison is done using relative and absolute tolerances. + * + * @param[in] F Deformation gradient. + * @param[in] S_ref Reference solution for PK2 stress. + * @param[in] rel_tol Relative tolerance for comparing S with S_ref. + * @param[in] abs_tol Absolute tolerance for comparing S with S_ref. + * @param[in] verbose Show values of F, S, and S_ref if true. + * @return None. + */ void testPK2StressAgainstReference(double F[3][3], double S_ref[3][3], double rel_tol, double abs_tol, bool verbose = false) { // Compute S(F) from get_pk2cc() double S[3][3], Dm[6][6]; get_pk2cc(F, S, Dm); - + // Compare S with reference solution - for (int i = 0; i < 3; i++){ - for (int j = 0; j < 3; j++){ + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { EXPECT_NEAR(S[i][j], S_ref[i][j], fmax(abs_tol, rel_tol * fabs(S_ref[i][j]))); } } - + // Print results if verbose if (verbose) { printMaterialParameters(); - + std::cout << "F =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -688,7 +870,7 @@ class TestMaterialModel { } std::cout << std::endl; } - + std::cout << "S =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -696,7 +878,7 @@ class TestMaterialModel { } std::cout << std::endl; } - + std::cout << "S_ref =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -708,36 +890,44 @@ class TestMaterialModel { } } - // Function to compare material elasticity tensor with reference solution - // ARGS: - // - F: Deformation gradient - // - CC_ref: Reference solution for material elasticity tensor - // - rel_tol: Relative tolerance for comparing CC with CC_ref + /** + * @brief Compares the material elasticity tensor CC(F) with a reference solution. + * + * This function computes the material elasticity tensor CC(F) from the deformation gradient F using get_pk2cc() + * and compares it with a reference solution CC_ref. The comparison is done using relative and absolute tolerances. + * + * @param[in] F Deformation gradient. + * @param[in] CC_ref Reference solution for material elasticity tensor. + * @param[in] rel_tol Relative tolerance for comparing CC with CC_ref. + * @param[in] abs_tol Absolute tolerance for comparing CC with CC_ref. + * @param[in] verbose Show values of F, CC, and CC_ref if true. + * @return None. + */ void testMaterialElasticityAgainstReference(double F[3][3], double CC_ref[3][3][3][3], double rel_tol, double abs_tol, bool verbose = false) { // Compute CC(F) from get_pk2cc() double S[3][3], Dm[6][6]; get_pk2cc(F, S, Dm); - + // Calculate CC from Dm double CC[3][3][3][3]; mat_models_carray::voigt_to_cc_carray<3>(Dm, CC); - + // Compare CC with reference solution - for (int i = 0; i < 3; i++){ - for (int j = 0; j < 3; j++){ - for (int k = 0; k < 3; k++){ - for (int l = 0; l < 3; l++){ + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + for (int k = 0; k < 3; k++) { + for (int l = 0; l < 3; l++) { EXPECT_NEAR(CC[i][j][k][l], CC_ref[i][j][k][l], fmax(abs_tol, rel_tol * fabs(CC_ref[i][j][k][l]))); } } } } - + // Print results if verbose if (verbose) { printMaterialParameters(); - + std::cout << "F =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -745,7 +935,7 @@ class TestMaterialModel { } std::cout << std::endl; } - + std::cout << "CC =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -758,7 +948,7 @@ class TestMaterialModel { } std::cout << std::endl; } - + std::cout << "CC_ref =" << std::endl; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { @@ -774,33 +964,101 @@ class TestMaterialModel { std::cout << std::endl; } } + + /** + * @brief Tests rho, beta, drho/dp, and dbeta/dp from g_vol_pen() against reference solutions. + * + * This function computes rho, beta, drho/dp, and dbeta/dp from the pressure p using g_vol_pen() + * and compares them with reference solutions. + * These values are required for treating a volumetric penalty term in the ustruct formulation. + * + * @param[in] p Pressure. + * @param[in] rho0 Initial solid density. + * @param[in] rho_ref Reference solution for rho. + * @param[in] beta_ref Reference solution for beta. + * @param[in] drhodp_ref Reference solution for drho/dp. + * @param[in] dbetadp_ref Reference solution for dbeta/dp. + * @param[in] rel_tol Relative tolerance for comparing rho and beta with reference solutions. + * @param[in] abs_tol Absolute tolerance for comparing rho and beta with reference solutions. + * @param[in] verbose Show values of p, rho, beta, rho_ref, beta_ref if true. + * @return None. + */ + void testRhoBetaAgainstReference(double p, double rho0, double rho_ref, double beta_ref, double drhodp_ref, double dbetadp_ref, double rel_tol, double abs_tol, bool verbose = false) { + double rho, beta, drhodp, dbetadp; + double Ja = 1.0; // Active strain Jacobian (not used in this function) + + // Compute rho, beta, drhodp, dbetadp from g_vol_pen() + g_vol_pen(p, rho0, rho, beta, drhodp, dbetadp, Ja); + + // Compare rho, beta, drho, dbeta with reference solutions + EXPECT_NEAR(rho, rho_ref, fmax(abs_tol, rel_tol * fabs(rho_ref))); + EXPECT_NEAR(beta, beta_ref, fmax(abs_tol, rel_tol * fabs(beta_ref))); + EXPECT_NEAR(drhodp, drhodp_ref, fmax(abs_tol, rel_tol * fabs(drhodp_ref))); + EXPECT_NEAR(dbetadp, dbetadp_ref, fmax(abs_tol, rel_tol * fabs(dbetadp_ref))); + + // Print results if verbose + if (verbose) { + printMaterialParameters(); + + std::cout << "p = " << p << std::endl; + std::cout << "rho0 = " << rho0 << std::endl; + std::cout << "rho = " << rho << ", rho_ref = " << rho_ref << std::endl; + std::cout << "beta = " << beta << ", beta_ref = " << beta_ref << std::endl; + std::cout << "drhodp = " << drhodp << ", drhodp_ref = " << drhodp_ref << std::endl; + std::cout << "dbetadp = " << dbetadp << ", dbetadp_ref = " << dbetadp_ref << std::endl; + std::cout << std::endl; + } + } }; // -------------------------------------------------------------- -// Class for test of Neo-Hookean material model +// --------------------- Material Model Classes ----------------- +// -------------------------------------------------------------- + +/** + * @brief Class for testing the Neo-Hookean material model. + * + * This class provides methods to set up and test the Neo-Hookean material model, including + * computing the strain energy and printing material parameters. + */ class TestNeoHookean : public TestMaterialModel { public: + /** + * @brief Parameters for the Neo-Hookean material model. + */ NeoHookeanParams params; + /** + * @brief Constructor for the TestNeoHookean class. + * + * Initializes the Neo-Hookean material parameters for svFSIplus. + * + * @param[in] params_ Parameters for the Neo-Hookean material model. + */ TestNeoHookean(const NeoHookeanParams ¶ms_) : TestMaterialModel( consts::ConstitutiveModelType::stIso_nHook, consts::ConstitutiveModelType::stVol_ST91), params(params_) { - // Set Neo-Hookean material parameters for svFSIplus auto &dmn = com_mod.mockEq.mockDmn; dmn.stM.C10 = params.C10; dmn.stM.Kpen = 0.0; // Zero volumetric penalty parameter } - // Print Neo-Hookean material parameters + /** + * @brief Prints the Neo-Hookean material parameters. + */ void printMaterialParameters() { std::cout << "C10 = " << params.C10 << std::endl; } - // Compute strain energy for Neo-Hookean material model + /** + * @brief Computes the strain energy for the Neo-Hookean material model. + * + * @param[in] F Deformation gradient. + * @return Strain energy density for the Neo-Hookean material model. + */ double computeStrainEnergy(const double F[3][3]) { - // Compute solid mechanics terms solidMechanicsTerms smTerms = calcSolidMechanicsTerms(F); @@ -812,17 +1070,30 @@ class TestNeoHookean : public TestMaterialModel { } }; -// -------------------------------------------------------------- -// Class for test of Mooney-Rivlin material model +/** + * @brief Class for testing the Mooney-Rivlin material model. + * + * This class provides methods to set up and test the Mooney-Rivlin material model, including + * computing the strain energy and printing material parameters. + */ class TestMooneyRivlin : public TestMaterialModel { public: + /** + * @brief Parameters for the Mooney-Rivlin material model. + */ MooneyRivlinParams params; + /** + * @brief Constructor for the TestMooneyRivlin class. + * + * Initializes the Mooney-Rivlin material parameters for svFSIplus. + * + * @param[in] params_ Parameters for the Mooney-Rivlin material model. + */ TestMooneyRivlin(const MooneyRivlinParams ¶ms_) : TestMaterialModel( consts::ConstitutiveModelType::stIso_MR, consts::ConstitutiveModelType::stVol_ST91), params(params_) { - // Set Mooney-Rivlin material parameters for svFSIplus auto &dmn = com_mod.mockEq.mockDmn; dmn.stM.C01 = params.C01; @@ -830,14 +1101,20 @@ class TestMooneyRivlin : public TestMaterialModel { dmn.stM.Kpen = 0.0; // Zero volumetric penalty parameter } - // Print Mooney-Rivlin material parameters + /** + * @brief Prints the Mooney-Rivlin material parameters. + */ void printMaterialParameters() { std::cout << "C01 = " << params.C01 << ", C10 = " << params.C10 << std::endl; } - // Compute strain energy for Mooney-Rivlin material model + /** + * @brief Computes the strain energy for the Mooney-Rivlin material model. + * + * @param[in] F Deformation gradient. + * @return Strain energy density for the Mooney-Rivlin material model. + */ double computeStrainEnergy(const double F[3][3]) { - // Compute solid mechanics terms solidMechanicsTerms smTerms = calcSolidMechanicsTerms(F); @@ -849,18 +1126,36 @@ class TestMooneyRivlin : public TestMaterialModel { } }; -// ---------------------------------------------------------------------------- -// Class for test of Holzapfel-Ogden material model + +/** + * @brief Class for testing the Holzapfel-Ogden material model. + * + * This class provides methods to set up and test the Holzapfel-Ogden material + * model, including computing the strain energy and printing material parameters. + * + * Implements two versions of the Holzapfel-Ogden model: + * - Full anisotropic invariants (I1_bar, I4_f, I4_s, I8_fs) used by cardiac mechanics benchmark paper (Arostica et al., 2024) + * - Modified anisotropic invariants (I1_bar, I4_bar_f, I4_bar_s, I8_bar_fs) used by svFSIplus + */ class TestHolzapfelOgden : public TestMaterialModel { public: + /** + * @brief Parameters for the Holzapfel-Ogden material model. + */ HolzapfelOgdenParams params; + /** + * @brief Constructor for the TestHolzapfelOgden class. + * + * Initializes the Holzapfel-Ogden material parameters for svFSIplus. + * + * @param[in] params_ Parameters for the Holzapfel-Ogden material model. + */ TestHolzapfelOgden(const HolzapfelOgdenParams ¶ms_) : TestMaterialModel( consts::ConstitutiveModelType::stIso_HO, consts::ConstitutiveModelType::stVol_ST91), - params(params_) + params(params_) { - - // Set Holzapfel-Ogden material parameters + // Set Holzapfel-Ogden material parameters for svFSIplus auto &dmn = com_mod.mockEq.mockDmn; dmn.stM.a = params.a; dmn.stM.b = params.b; @@ -870,7 +1165,7 @@ class TestHolzapfelOgden : public TestMaterialModel { dmn.stM.bss = params.b_s; dmn.stM.afs = params.a_fs; dmn.stM.bfs = params.b_fs; - dmn.stM.Kpen = params.kappa; + dmn.stM.Kpen = 0.0; // Zero volumetric penalty parameter // Set number of fiber directions and fiber directions nFn = 2; @@ -880,7 +1175,9 @@ class TestHolzapfelOgden : public TestMaterialModel { fN.set_col(1, s); } - // Print Holzapfel-Ogden material parameters + /** + * @brief Prints the Holzapfel-Ogden material parameters. + */ void printMaterialParameters() { std::cout << "a = " << params.a << std::endl; std::cout << "b = " << params.b << std::endl; @@ -890,28 +1187,33 @@ class TestHolzapfelOgden : public TestMaterialModel { std::cout << "b_s = " << params.b_s << std::endl; std::cout << "a_fs = " << params.a_fs << std::endl; std::cout << "b_fs = " << params.b_fs << std::endl; - std::cout << "kappa = " << params.kappa << std::endl; std::cout << "k = " << params.k << std::endl; std::cout << "f = " << "[" << params.f[0] << " " << params.f[1] << " " << params.f[2] << "]" << std::endl; std::cout << "s = " << "[" << params.s[0] << " " << params.s[1] << " " << params.s[2] << "]" << std::endl; std::cout << "full_anisotropic_invariants = " << params.full_anisotropic_invariants << std::endl; } - // Smooth Heaviside function centered at 1 + /** + * @brief Smoothed Heaviside function centered at x = 1. + * + * @param[in] x Input value. + * @param[in] k Smoothing parameter. + * @return Smoothed Heaviside function. + */ double chi(const double x, const double k=100) const { return 1. / (1. + exp(-k * (x - 1.))); } - // Compute strain energy for Holzapfel-Ogden material model + /** + * @brief Computes the strain energy for the Holzapfel-Ogden material model. + * + * @param[in] F Deformation gradient. + * @return Strain energy density for the Holzapfel-Ogden material model. + */ double computeStrainEnergy(const double F[3][3]) { - // Compute solid mechanics terms solidMechanicsTerms smTerms = calcSolidMechanicsTerms(F); - // Strain energy density for Holzapfel-Ogden material model - // Psi = a/2b * exp{b(I1_bar - 3)} + Sum_{i=f,s} [a_i/2b_i * chi(I4_i) * (exp{b_i(I4_i - 1)^2} - 1)] - // + a_fs/2b_fs * (exp{b_fs*I8_fs^2} - 1) + kappa/4 * (J^2 - 1 - 2*ln(J)) - // Material parameters double a = params.a; double b = params.b; @@ -921,7 +1223,6 @@ class TestHolzapfelOgden : public TestMaterialModel { double b_s = params.b_s; double a_fs = params.a_fs; double b_fs = params.b_fs; - double kappa = params.kappa; // Smoothed Heaviside parameter double k = params.k; @@ -930,8 +1231,14 @@ class TestHolzapfelOgden : public TestMaterialModel { double f[3] = {params.f[0], params.f[1], params.f[2]}; double s[3] = {params.s[0], params.s[1], params.s[2]}; + // Strain energy density for Holzapfel-Ogden material model + // Formulation used by cardiac mechanics benchmark paper (Arostica et al., 2024) // Uses I1_bar (bar = isochoric), but I4_f, I4_s, I8_fs (not bar) + // Psi = a/2b * exp{b(I1_bar - 3)} + // + a_f/2b_f * chi(I4_f) * (exp{b_f(I4_f - 1)^2} - 1 + // + a_s/2b_s * chi(I4_s) * (exp{b_s(I4_s - 1)^2} - 1 + // + a_fs/2b_fs * (exp{b_fs*I8_fs^2} - 1) // This corresponds to the HO-ma (modified anisotropy) implementation in svFSIplus if (params.full_anisotropic_invariants) { // Invariants @@ -951,12 +1258,15 @@ class TestHolzapfelOgden : public TestMaterialModel { Psi += a_f / (2.0 * b_f) * chi(I4_f, k) * (exp(b_f * pow(I4_f - 1.0, 2)) - 1.0); // Fiber term Psi += a_s / (2.0 * b_s) * chi(I4_s, k) * (exp(b_s * pow(I4_s - 1.0, 2)) - 1.0); // Sheet term Psi += a_fs / (2.0 * b_fs) * (exp(b_fs * pow(I8_fs, 2)) - 1.0); // Cross-fiber term - Psi += kappa / 4.0 * (pow(smTerms.J, 2) - 1.0 - 2.0 * log(smTerms.J)); // Simo-Taylor 91 volumetric penalty term return Psi; } // Formulation with fully decoupled isochoric-volumetric split // Uses I1_bar, I4_bar_f, I4_bar_s, I8_bar_fs (bar = isochoric) + // Psi = a/2b * exp{b(I1_bar - 3)} + // + a_f/2b_f * chi(I4_bar_f) * (exp{b_f(I4_bar_f - 1)^2} - 1 + // + a_s/2b_s * chi(I4_bar_s) * (exp{b_s(I4_bar_s - 1)^2} - 1 + // + a_fs/2b_fs * (exp{b_fs*I8_bar_fs^2} - 1) // This corresponds to the HO-d (decoupled) implementation in svFSIplus else { // Invariants @@ -976,9 +1286,186 @@ class TestHolzapfelOgden : public TestMaterialModel { Psi += a_f / (2.0 * b_f) * chi(I4_bar_f, k) * (exp(b_f * pow(I4_bar_f - 1.0, 2)) - 1.0); // Fiber term Psi += a_s / (2.0 * b_s) * chi(I4_bar_s, k) * (exp(b_s * pow(I4_bar_s - 1.0, 2)) - 1.0); // Sheet term Psi += a_fs / (2.0 * b_fs) * (exp(b_fs * pow(I8_bar_fs, 2)) - 1.0); // Cross-fiber term - Psi += kappa / 4.0 * (pow(smTerms.J, 2) - 1.0 - 2.0 * log(smTerms.J)); // Volumetric penalty term return Psi; } } }; + + +/** + * @brief Class for testing the quadratic volumetric penalty model. + * + * This class provides methods to set up and test the quadratic volumetric penalty model, including + * computing the strain energy and printing material parameters. + */ +class TestQuadraticVolumetricPenalty : public TestMaterialModel { +public: + + /** + * @brief Parameters for the volumetric penalty model. + */ + VolumetricPenaltyParams params; + + /** + * @brief Constructor for the TestQuadraticVolumetricPenalty class. + * + * Initializes the volumetric penalty parameters for svFSIplus. + * + * @param[in] params_ Parameters for the volumetric penalty model. + */ + TestQuadraticVolumetricPenalty(const VolumetricPenaltyParams ¶ms_) : TestMaterialModel( consts::ConstitutiveModelType::stIso_nHook, consts::ConstitutiveModelType::stVol_Quad), + params(params_) + { + + // Set volumetric penalty parameter for svFSIplus + auto &dmn = com_mod.mockEq.mockDmn; + dmn.stM.Kpen = params.kappa; // Volumetric penalty parameter + + // Note: Use Neo-Hookean material model for isochoric part, but set parameters to zero + dmn.stM.C10 = 0.0; // Zero Neo-Hookean parameter + } + + /** + * @brief Prints the volumetric penalty parameters. + */ + void printMaterialParameters() { + std::cout << "kappa = " << params.kappa << std::endl; + } + + /** + * @brief Computes the strain energy for the quadratic volumetric penalty model. + * + * @param[in] F Deformation gradient. + * @return Strain energy density for the quadratic volumetric penalty model. + */ + double computeStrainEnergy(const double F[3][3]) { + + // Compute solid mechanics terms + solidMechanicsTerms smTerms = calcSolidMechanicsTerms(F); + + // Strain energy density for quadratic volumetric penalty model + // Psi = kappa/2 * (J - 1)^2 + double Psi = params.kappa/2.0 * pow(smTerms.J - 1.0, 2); + + return Psi; + } +}; + +/** + * @brief Class for testing the Simo-Taylor91 volumetric penalty model. + * + * This class provides methods to set up and test the Simo-Taylor91 volumetric penalty model, including + * computing the strain energy and printing material parameters. + */ +class TestSimoTaylor91VolumetricPenalty : public TestMaterialModel { +public: + + /** + * @brief Parameters for the volumetric penalty model. + */ + VolumetricPenaltyParams params; + + /** + * @brief Constructor for the TestSimoTaylor91VolumetricPenalty class. + * + * Initializes the volumetric penalty parameters for svFSIplus. + * + * @param[in] params_ Parameters for the volumetric penalty model. + */ + TestSimoTaylor91VolumetricPenalty(const VolumetricPenaltyParams ¶ms_) : TestMaterialModel( consts::ConstitutiveModelType::stIso_nHook, consts::ConstitutiveModelType::stVol_ST91), + params(params_) + { + + // Set volumetric penalty parameter for svFSIplus + auto &dmn = com_mod.mockEq.mockDmn; + dmn.stM.Kpen = params.kappa; // Volumetric penalty parameter + + // Note: Use Neo-Hookean material model for isochoric part, but set parameters to zero + dmn.stM.C10 = 0.0; // Zero Neo-Hookean parameter + } + + /** + * @brief Prints the volumetric penalty parameters. + */ + void printMaterialParameters() { + std::cout << "kappa = " << params.kappa << std::endl; + } + + /** + * @brief Computes the strain energy for the Simo-Taylor91 volumetric penalty model. + * + * @param[in] F Deformation gradient. + * @return Strain energy density for the Simo-Taylor91 volumetric penalty model. + */ + double computeStrainEnergy(const double F[3][3]) { + + // Compute solid mechanics terms + solidMechanicsTerms smTerms = calcSolidMechanicsTerms(F); + + // Strain energy density for Simo-Taylor91 volumetric penalty model + // Psi = kappa/4 * (J^2 - 1 - 2*ln(J)) + double Psi = params.kappa/4.0 * (pow(smTerms.J, 2) - 1.0 - 2.0 * log(smTerms.J)); + + return Psi; + } +}; + +/** + * @brief Class for testing the Miehe94 volumetric penalty model. + * + * This class provides methods to set up and test the Miehe94 volumetric penalty model, including + * computing the strain energy and printing material parameters. + */ +class TestMiehe94VolumetricPenalty : public TestMaterialModel { +public: + + /** + * @brief Parameters for the volumetric penalty model. + */ + VolumetricPenaltyParams params; + + /** + * @brief Constructor for the TestMiehe94VolumetricPenalty class. + * + * Initializes the volumetric penalty parameters for svFSIplus. + * + * @param[in] params_ Parameters for the volumetric penalty model. + */ + TestMiehe94VolumetricPenalty(const VolumetricPenaltyParams ¶ms_) : TestMaterialModel( consts::ConstitutiveModelType::stIso_nHook, consts::ConstitutiveModelType::stVol_M94), + params(params_) + { + + // Set volumetric penalty parameter for svFSIplus + auto &dmn = com_mod.mockEq.mockDmn; + dmn.stM.Kpen = params.kappa; // Volumetric penalty parameter + + // Note: Use Neo-Hookean material model for isochoric part, but set parameters to zero + dmn.stM.C10 = 0.0; // Zero Neo-Hookean parameter + } + + /** + * @brief Prints the volumetric penalty parameters. + */ + void printMaterialParameters() { + std::cout << "kappa = " << params.kappa << std::endl; + } + + /** + * @brief Computes the strain energy for the Miehe94 volumetric penalty model. + * + * @param[in] F Deformation gradient. + * @return Strain energy density for the Miehe94 volumetric penalty model. + */ + double computeStrainEnergy(const double F[3][3]) { + + // Compute solid mechanics terms + solidMechanicsTerms smTerms = calcSolidMechanicsTerms(F); + + // Strain energy density for Miehe94 volumetric penalty model + // Psi = kappa * (J - ln(J) - 1) + double Psi = params.kappa * (smTerms.J - log(smTerms.J) - 1.0); + + return Psi; + } +}; \ No newline at end of file