From f9fe6fe4b5025c54188aad4c1e0859d578230d70 Mon Sep 17 00:00:00 2001 From: divyaadil23 <146143076+divyaadil23@users.noreply.github.com> Date: Wed, 16 Oct 2024 10:04:34 -0700 Subject: [PATCH] Holzapfel Ogden Modified Anisotropy (HO-ma) Constitutive Model (#272) * Added the Holzapfel Ogden Modified Anisotropy (HO-ma) Constitutive Model to material models. Updated Holzapfel Ogden Model. Both HO and HO-ma models are taken from https://github.com/vvedula22/svFSI/blob/master/Code/Source/svFSI/MATMODELS.f and modified. Added additional terms containing the second derivative of smoothed heaviside funciton to the elasticity tensor to make it exact. Changed the approximate derivatives of heaviside function to exact. Removed fds (dot product of fiber directions in 8th invariant). Tested it with the new framework of unit tests - the unit tests now test the order of convergence for finite difference estimates. Added some integrated tests for struct and ustruct for both models. * Changed input file names to svFSIplus.xml and responded to changes in comments. * Added the Holzapfel Ogden Modified Anisotropy (HO-ma) Constitutive Model to material models. Updated Holzapfel Ogden Model. Both HO and HO-ma models are taken from https://github.com/vvedula22/svFSI/blob/master/Code/Source/svFSI/MATMODELS.f and modified. Added additional terms containing the second derivative of smoothed heaviside funciton to the elasticity tensor to make it exact. Changed the approximate derivatives of heaviside function to exact. Removed fds (dot product of fiber directions in 8th invariant). Tested it with the new framework of unit tests - the unit tests now test the order of convergence for finite difference estimates. Added some integrated tests for struct and ustruct for both models. * Changed input file names to svFSIplus.xml and responded to changes in comments. * Removed redundant norm function and updated result file without solid viscosity. * Removed commented solid viscosity from input files * Removing extra files and updating tests/README.md with the current version * Removing extra files and updating README.md with the current version * Undoing meaningless changes to README --------- Co-authored-by: Martin R. Pfaller Co-authored-by: Aaron Brown <71424733+aabrown100-git@users.noreply.github.com> --- Code/Source/svFSI/ComMod.h | 3 + Code/Source/svFSI/Parameters.cpp | 12 +- Code/Source/svFSI/Parameters.h | 3 +- Code/Source/svFSI/consts.cpp | 5 +- Code/Source/svFSI/distribute.cpp | 2 +- Code/Source/svFSI/mat_fun_carray.h | 1 + Code/Source/svFSI/mat_models.cpp | 345 ++- Code/Source/svFSI/mat_models_carray.h | 359 ++- Code/Source/svFSI/post.cpp | 2 + Code/Source/svFSI/read_files.cpp | 2 +- Code/Source/svFSI/set_material_props.h | 20 + .../endo_pressure.dat | 0 .../mesh/fibersLongCells.vtu | 0 .../mesh/fibersSheetCells.vtu | 0 .../mesh/mesh-complete.exterior.vtp | 0 .../mesh/mesh-complete.mesh.vtu | 0 .../mesh/mesh-surfaces/endo.vtp | 0 .../mesh/mesh-surfaces/epi.vtp | 0 .../mesh/mesh-surfaces/top.vtp | 0 .../result_001.vtu | 3 + .../svFSIplus.xml | 21 +- .../endo_pressure.dat | 3 + .../mesh/fibersLongCells.vtu | 3 + .../mesh/fibersSheetCells.vtu | 3 + .../mesh/mesh-complete.exterior.vtp | 3 + .../mesh/mesh-complete.mesh.vtu | 3 + .../mesh/mesh-surfaces/endo.vtp | 3 + .../mesh/mesh-surfaces/epi.vtp | 3 + .../mesh/mesh-surfaces/top.vtp | 3 + .../LV_HolzapfelOgden_passive/result_001.vtu | 3 + .../LV_HolzapfelOgden_passive/svFSIplus.xml | 112 + .../LV_Holzapfel_passive/result_001.vtu | 3 - .../endo_pressure.dat | 3 + .../mesh/fibersLongCells.vtu | 3 + .../mesh/fibersSheetCells.vtu | 3 + .../mesh/mesh-complete.exterior.vtp | 3 + .../mesh/mesh-complete.mesh.vtu | 3 + .../mesh/mesh-surfaces/endo.vtp | 3 + .../mesh/mesh-surfaces/epi.vtp | 3 + .../mesh/mesh-surfaces/top.vtp | 3 + .../result_001.vtu | 3 + .../svFSIplus.xml | 113 + .../endo_pressure.dat | 3 + .../mesh/fibersLongCells.vtu | 3 + .../mesh/fibersSheetCells.vtu | 3 + .../mesh/mesh-complete.exterior.vtp | 3 + .../mesh/mesh-complete.mesh.vtu | 3 + .../mesh/mesh-surfaces/endo.vtp | 3 + .../mesh/mesh-surfaces/epi.vtp | 3 + .../mesh/mesh-surfaces/top.vtp | 3 + .../LV_HolzapfelOgden_passive/result_001.vtu | 3 + .../LV_HolzapfelOgden_passive/svFSIplus.xml | 113 + tests/test_struct.py | 7 +- tests/test_ustruct.py | 8 + tests/unitTests/test.cpp | 2373 ++++++++++++----- tests/unitTests/test.h | 1336 +++++++--- 56 files changed, 3805 insertions(+), 1119 deletions(-) rename tests/cases/struct/{LV_Holzapfel_passive => LV_HolzapfelOgdenModifiedAnisotropy_passive}/endo_pressure.dat (100%) rename tests/cases/struct/{LV_Holzapfel_passive => LV_HolzapfelOgdenModifiedAnisotropy_passive}/mesh/fibersLongCells.vtu (100%) rename tests/cases/struct/{LV_Holzapfel_passive => LV_HolzapfelOgdenModifiedAnisotropy_passive}/mesh/fibersSheetCells.vtu (100%) rename tests/cases/struct/{LV_Holzapfel_passive => LV_HolzapfelOgdenModifiedAnisotropy_passive}/mesh/mesh-complete.exterior.vtp (100%) rename tests/cases/struct/{LV_Holzapfel_passive => LV_HolzapfelOgdenModifiedAnisotropy_passive}/mesh/mesh-complete.mesh.vtu (100%) rename tests/cases/struct/{LV_Holzapfel_passive => LV_HolzapfelOgdenModifiedAnisotropy_passive}/mesh/mesh-surfaces/endo.vtp (100%) rename tests/cases/struct/{LV_Holzapfel_passive => LV_HolzapfelOgdenModifiedAnisotropy_passive}/mesh/mesh-surfaces/epi.vtp (100%) rename tests/cases/struct/{LV_Holzapfel_passive => LV_HolzapfelOgdenModifiedAnisotropy_passive}/mesh/mesh-surfaces/top.vtp (100%) create mode 100644 tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/result_001.vtu rename tests/cases/struct/{LV_Holzapfel_passive => LV_HolzapfelOgdenModifiedAnisotropy_passive}/svFSIplus.xml (86%) create mode 100755 tests/cases/struct/LV_HolzapfelOgden_passive/endo_pressure.dat create mode 100644 tests/cases/struct/LV_HolzapfelOgden_passive/mesh/fibersLongCells.vtu create mode 100644 tests/cases/struct/LV_HolzapfelOgden_passive/mesh/fibersSheetCells.vtu create mode 100644 tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-complete.exterior.vtp create mode 100644 tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-complete.mesh.vtu create mode 100644 tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/endo.vtp create mode 100644 tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/epi.vtp create mode 100644 tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/top.vtp create mode 100644 tests/cases/struct/LV_HolzapfelOgden_passive/result_001.vtu create mode 100644 tests/cases/struct/LV_HolzapfelOgden_passive/svFSIplus.xml delete mode 100644 tests/cases/struct/LV_Holzapfel_passive/result_001.vtu create mode 100755 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/endo_pressure.dat create mode 100644 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersLongCells.vtu create mode 100644 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersSheetCells.vtu create mode 100644 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.exterior.vtp create mode 100644 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.mesh.vtu create mode 100644 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/endo.vtp create mode 100644 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/epi.vtp create mode 100644 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/top.vtp create mode 100644 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/result_001.vtu create mode 100644 tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/svFSIplus.xml create mode 100755 tests/cases/ustruct/LV_HolzapfelOgden_passive/endo_pressure.dat create mode 100644 tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/fibersLongCells.vtu create mode 100644 tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/fibersSheetCells.vtu create mode 100644 tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-complete.exterior.vtp create mode 100644 tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-complete.mesh.vtu create mode 100644 tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/endo.vtp create mode 100644 tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/epi.vtp create mode 100644 tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/top.vtp create mode 100644 tests/cases/ustruct/LV_HolzapfelOgden_passive/result_001.vtu create mode 100644 tests/cases/ustruct/LV_HolzapfelOgden_passive/svFSIplus.xml diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h index f5bb6403..9f4d1924 100644 --- a/Code/Source/svFSI/ComMod.h +++ b/Code/Source/svFSI/ComMod.h @@ -333,6 +333,9 @@ class fibStrsType // Constant steady value double g = 0.0; + // Cross fiber stress parameter + double eta_s = 0.0; + // Unsteady time-dependent values fcType gt; }; diff --git a/Code/Source/svFSI/Parameters.cpp b/Code/Source/svFSI/Parameters.cpp index bcea3dde..8d741c52 100644 --- a/Code/Source/svFSI/Parameters.cpp +++ b/Code/Source/svFSI/Parameters.cpp @@ -467,7 +467,8 @@ const std::string ConstitutiveModelParameters::xml_element_name_ = "Constitutive // [TODO] Should use the types defined in consts.h. const std::string ConstitutiveModelParameters::GUCCIONE_MODEL = "Guccione"; const std::string ConstitutiveModelParameters::HGO_MODEL = "HGO"; -const std::string ConstitutiveModelParameters::HOLZAPFEL_MODEL = "Holzapfel"; +const std::string ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL = "HolzapfelOgden"; +const std::string ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL = "HolzapfelOgden-ModifiedAnisotropy"; const std::string ConstitutiveModelParameters::LEE_SACKS = "Lee-Sacks"; const std::string ConstitutiveModelParameters::NEOHOOKEAN_MODEL = "neoHookean"; const std::string ConstitutiveModelParameters::STVENANT_KIRCHHOFF_MODEL = "stVenantKirchhoff"; @@ -479,7 +480,9 @@ const std::map ConstitutiveModelParameters::constituti {ConstitutiveModelParameters::HGO_MODEL, ConstitutiveModelParameters::HGO_MODEL}, - {ConstitutiveModelParameters::HOLZAPFEL_MODEL, ConstitutiveModelParameters::HOLZAPFEL_MODEL}, + {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL, ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL}, + + {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL, ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL}, {ConstitutiveModelParameters::LEE_SACKS, ConstitutiveModelParameters::LEE_SACKS}, @@ -498,7 +501,8 @@ using SetConstitutiveModelParamMapType = std::map void {cp->guccione.set_values(params);}}, {ConstitutiveModelParameters::HGO_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->holzapfel_gasser_ogden.set_values(params);}}, - {ConstitutiveModelParameters::HOLZAPFEL_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->holzapfel.set_values(params);}}, + {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->holzapfel.set_values(params);}}, + {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->holzapfel.set_values(params);}}, {ConstitutiveModelParameters::LEE_SACKS, [](CmpType cp, CmpXmlType params) -> void {cp->lee_sacks.set_values(params);}}, {ConstitutiveModelParameters::NEOHOOKEAN_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->neo_hookean.set_values(params);}}, {ConstitutiveModelParameters::STVENANT_KIRCHHOFF_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->stvenant_kirchhoff.set_values(params);}}, @@ -510,6 +514,8 @@ using PrintConstitutiveModelParamMapType = std::map void {cp->guccione.print_parameters();}}, {ConstitutiveModelParameters::HGO_MODEL, [](CmpType cp) -> void {cp->holzapfel_gasser_ogden.print_parameters();}}, + {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL, [](CmpType cp) -> void {cp->holzapfel.print_parameters();}}, + {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL, [](CmpType cp) -> void {cp->holzapfel.print_parameters();}}, {ConstitutiveModelParameters::LEE_SACKS, [](CmpType cp) -> void {cp->lee_sacks.print_parameters();}}, {ConstitutiveModelParameters::NEOHOOKEAN_MODEL, [](CmpType cp) -> void {cp->neo_hookean.print_parameters();}}, {ConstitutiveModelParameters::STVENANT_KIRCHHOFF_MODEL, [](CmpType cp) -> void {cp->stvenant_kirchhoff.print_parameters();}}, diff --git a/Code/Source/svFSI/Parameters.h b/Code/Source/svFSI/Parameters.h index 4e1a7830..64de4804 100644 --- a/Code/Source/svFSI/Parameters.h +++ b/Code/Source/svFSI/Parameters.h @@ -517,7 +517,8 @@ class ConstitutiveModelParameters : public ParameterLists // Model types supported. static const std::string GUCCIONE_MODEL; static const std::string HGO_MODEL; - static const std::string HOLZAPFEL_MODEL; + static const std::string HOLZAPFEL_OGDEN_MODEL; + static const std::string HOLZAPFEL_OGDEN_MA_MODEL; static const std::string LEE_SACKS; static const std::string NEOHOOKEAN_MODEL; static const std::string STVENANT_KIRCHHOFF_MODEL; diff --git a/Code/Source/svFSI/consts.cpp b/Code/Source/svFSI/consts.cpp index bd2e521c..698850ef 100644 --- a/Code/Source/svFSI/consts.cpp +++ b/Code/Source/svFSI/consts.cpp @@ -86,7 +86,10 @@ const std::map constitutive_model_name_to_typ {"Gucci", ConstitutiveModelType::stIso_Gucci}, {"HO", ConstitutiveModelType::stIso_HO}, - {"Holzapfel", ConstitutiveModelType::stIso_HO}, + {"HolzapfelOgden", ConstitutiveModelType::stIso_HO}, + + {"HO_ma", ConstitutiveModelType::stIso_HO_ma}, + {"HolzapfelOgden-ModifiedAnisotropy", ConstitutiveModelType::stIso_HO_ma}, {"quad", ConstitutiveModelType::stVol_Quad}, {"Quad", ConstitutiveModelType::stVol_Quad}, diff --git a/Code/Source/svFSI/distribute.cpp b/Code/Source/svFSI/distribute.cpp index 93532359..09bc061a 100644 --- a/Code/Source/svFSI/distribute.cpp +++ b/Code/Source/svFSI/distribute.cpp @@ -1221,7 +1221,7 @@ void dist_mat_consts(const ComMod& com_mod, const CmMod& cm_mod, const cmType& c cm.bcast(cm_mod, lStM.Tf.gt.r, "lStM.Tf.gt.r"); cm.bcast(cm_mod, lStM.Tf.gt.i, "lStM.Tf.gt.i"); } - + cm.bcast(cm_mod, &lStM.Tf.eta_s); } diff --git a/Code/Source/svFSI/mat_fun_carray.h b/Code/Source/svFSI/mat_fun_carray.h index 5069199d..212f876f 100644 --- a/Code/Source/svFSI/mat_fun_carray.h +++ b/Code/Source/svFSI/mat_fun_carray.h @@ -459,6 +459,7 @@ double norm(const Vector& u, const double v[N]) } + template void mat_symm(const double A[N][N], double S[N][N]) { diff --git a/Code/Source/svFSI/mat_models.cpp b/Code/Source/svFSI/mat_models.cpp index b6cdf6a1..30b060aa 100644 --- a/Code/Source/svFSI/mat_models.cpp +++ b/Code/Source/svFSI/mat_models.cpp @@ -34,6 +34,7 @@ #include "fft.h" #include "mat_fun.h" +#include "mat_fun_carray.h" #include "utils.h" #include @@ -179,6 +180,7 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn // Fiber-reinforced stress double Tfa = 0.0; get_fib_stress(com_mod, cep_mod, stM.Tf, Tfa); + double Tsa = Tfa*stM.Tf.eta_s; // Electromechanics coupling - active stress if (cep_mod.cem.aStress) { @@ -428,43 +430,67 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn double Ess = Inv6 - 1.0; double Efs = Inv8; + // Smoothed Heaviside function: 1 / (1 + exp(-kx)) = 1 - 1 / (1 + exp(kx)) + double k = stM.khs; + double one_over_exp_plus_one_f = 1.0 / (exp(k * Eff) + 1.0); + double one_over_exp_plus_one_s = 1.0 / (exp(k * Ess) + 1.0); + double c4f = 1.0 - one_over_exp_plus_one_f; + double c4s = 1.0 - one_over_exp_plus_one_s; + + // Exact first derivative of smoothed heaviside function (from Wolfram Alpha) + double dc4f = k * (one_over_exp_plus_one_f - pow(one_over_exp_plus_one_f,2)); + double dc4s = k * (one_over_exp_plus_one_s - pow(one_over_exp_plus_one_s,2)); + + // Exact second derivative of smoothed heaviside function (from Wolfram Alpha) + double ddc4f = pow(k,2) * (-one_over_exp_plus_one_f + 3.0*pow(one_over_exp_plus_one_f,2) - 2.0*pow(one_over_exp_plus_one_f,3)); + double ddc4s = pow(k,2) * (-one_over_exp_plus_one_s + 3.0*pow(one_over_exp_plus_one_s,2) - 2.0*pow(one_over_exp_plus_one_s,3)); + + // Isotropic + fiber-sheet interaction stress double g1 = stM.a * exp(stM.b*(Inv1-3.0)); double g2 = 2.0 * stM.afs * Efs * exp(stM.bfs*Efs*Efs); auto Hfs = mat_symm_prod(fl.col(0), fl.col(1), nsd); auto Sb = g1*Idm + g2*Hfs; - Efs = Efs * Efs; + // Isotropic + fiber-sheet interaction stiffness g1 = 2.0*J4d*stM.b*g1; - g2 = 4.0*J4d*stM.afs*(1.0 + 2.0*stM.bfs*Efs)* exp(stM.bfs*Efs); + g2 = 4.0*J4d*stM.afs*(1.0 + 2.0*stM.bfs*Efs*Efs)* exp(stM.bfs*Efs*Efs); auto CCb = g1 * ten_dyad_prod(Idm, Idm, nsd) + g2 * ten_dyad_prod(Hfs, Hfs, nsd); - // Fiber reinforcement/active stress - if (Eff > 0.0) { - g1 = Tfa; - - g1 = g1 + 2.0 * stM.aff * Eff * exp(stM.bff*Eff*Eff); - auto Hff = mat_dyad_prod(fl.col(0), fl.col(0), nsd); - Sb = Sb + g1*Hff; - - Eff = Eff * Eff; - g1 = 4.0*J4d*stM.aff*(1.0 + 2.0*stM.bff*Eff)*exp(stM.bff*Eff); - CCb = CCb + g1*ten_dyad_prod(Hff, Hff, nsd); - } - - if (Ess > 0.0) { - g2 = 2.0 * stM.ass * Ess * exp(stM.bss*Ess*Ess); - auto Hss = mat_dyad_prod(fl.col(1), fl.col(1), nsd); - Sb = Sb + g2*Hss; + // Fiber-fiber interaction stress + additional reinforcement (Tfa) + double rexp = exp(stM.bff*Eff*Eff); + g1 = c4f * Eff * rexp; + g1 = g1 + (0.5*dc4f/stM.bff) * (rexp - 1.0); + g1 = 2.0 * stM.aff * g1 + Tfa; + auto Hff = mat_dyad_prod(fl.col(0), fl.col(0), nsd); + Sb = Sb + g1*Hff; + + // Fiber-fiber interaction stiffness + g1 = c4f * (1.0 + 2.0*stM.bff*Eff*Eff); + g1 = (g1 + 2.0*dc4f*Eff) * rexp; + g1 = g1 + (0.5*ddc4f/stM.bff)*(rexp - 1.0); + g1 = 4.0 * J4d * stM.aff * g1; + CCb = CCb + g1*ten_dyad_prod(Hff, Hff, nsd); + + // Sheet-sheet interaction stress + additional cross-fiber stress + rexp = exp(stM.bss*Ess*Ess); + g2 = c4s * Ess * rexp; + g2 = g2 + (0.5*dc4s/stM.bss) * (rexp - 1.0); + g2 = 2.0 * stM.ass * g2 + Tsa; + auto Hss = mat_dyad_prod(fl.col(1), fl.col(1), nsd); + Sb = Sb + g2 * Hss; - Ess = Ess * Ess; - g2 = 4.0*J4d*stM.ass*(1.0 + 2.0*stM.bss*Ess)*exp(stM.bss*Ess); - CCb = CCb + g2*ten_dyad_prod(Hss, Hss, nsd); - } + // Sheet-sheet interaction stiffness + g2 = c4s * (1.0 + 2.0 * stM.bss * Ess * Ess); + g2 = (g2 + 2.0*dc4s*Ess) * rexp; + g2 = g2 + (0.5*ddc4s/stM.bss)*(rexp - 1.0); + g2 = 4.0 * J4d * stM.ass * g2; + CCb = CCb + g2*ten_dyad_prod(Hss, Hss, nsd); + // Isochoric 2nd-Piola-Kirchhoff stress and stiffness tensors double r1 = J2d*mat_ddot(C, Sb, nsd) / nd; 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); CC = ten_transpose(CC, nsd); CC = ten_ddot(PP, CC, nsd); @@ -483,6 +509,99 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn } } break; + // HO (Holzapfel-Ogden)-MA model for myocardium with full invariants for the anisotropy terms (modified-anisotropy) + case ConstitutiveModelType::stIso_HO_ma: { + if (nfd != 2) { + //err = "Min fiber directions not defined for Holzapfel material model (2)" + } + double Inv4 = norm(fl.col(0), mat_mul(C, fl.col(0))); + double Inv6 = norm(fl.col(1), mat_mul(C, fl.col(1))); + double Inv8 = norm(fl.col(0), mat_mul(C, fl.col(1))); + + //double fds = norm(fl.col(0),fl.col(1)); + double Eff = Inv4 - 1.0; + double Ess = Inv6 - 1.0; + double Efs = Inv8; + + // Smoothed Heaviside function: 1 / (1 + exp(-kx)) = 1 - 1 / (1 + exp(kx)) + double k = stM.khs; + double one_over_exp_plus_one_f = 1.0 / (exp(k * Eff) + 1.0); + double one_over_exp_plus_one_s = 1.0 / (exp(k * Ess) + 1.0); + double c4f = 1.0 - one_over_exp_plus_one_f; + double c4s = 1.0 - one_over_exp_plus_one_s; + + // Approx. derivative of smoothed heaviside function + //double dc4f = 0.25*stM.khs*exp(-stM.khs*abs(Eff)); + //double dc4s = 0.25*stM.khs*exp(-stM.khs*abs(Ess)); + + // Exact first derivative of smoothed heaviside function (from Wolfram Alpha) + double dc4f = k * (one_over_exp_plus_one_f - pow(one_over_exp_plus_one_f,2)); + double dc4s = k * (one_over_exp_plus_one_s - pow(one_over_exp_plus_one_s,2)); + + // Exact second derivative of smoothed heaviside function (from Wolfram Alpha) + double ddc4f = pow(k,2) * (-one_over_exp_plus_one_f + 3.0*pow(one_over_exp_plus_one_f,2) - 2.0*pow(one_over_exp_plus_one_f,3)); + double ddc4s = pow(k,2) * (-one_over_exp_plus_one_s + 3.0*pow(one_over_exp_plus_one_s,2) - 2.0*pow(one_over_exp_plus_one_s,3)); + + // Isochoric stress and stiffness + double g1 = stM.a * exp(stM.b*(Inv1-3.0)); + auto Sb = g1*Idm; + double r1 = J2d/nd*mat_fun::mat_ddot(C, Sb, nsd); + + g1 = g1*2.0*J4d*stM.b; + auto CCb = g1 * ten_dyad_prod(Idm, Idm, nsd); + + // Add isochoric stress and stiffness contribution + 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); + CC = ten_transpose(CC, nsd); + CC = ten_ddot(PP, CC, nsd); + CC = CC + 2.0*r1 * (ten_symm_prod(Ci, Ci, nsd) - + 1.0/nd * ten_dyad_prod(Ci, Ci, nsd)) + - 2.0/nd * (ten_dyad_prod(Ci, S, nsd) + ten_dyad_prod(S, Ci, nsd)); + + // Now add aniostropic components + // Fiber-sheet interaction terms + g1 = 2.0 * stM.afs * exp(stM.bfs*Efs*Efs); + auto Hfs = mat_symm_prod(fl.col(0), fl.col(1), nsd); + S = S + g1*Efs*Hfs; + + g1 = g1 * 2.0*(1.0 + 2.0*stM.bfs*Efs*Efs); + + CC = CC + g1*ten_dyad_prod(Hfs, Hfs, nsd); + + // Fiber-fiber interaction stress + additional reinforcement (Tfa) + double rexp = exp(stM.bff * Eff * Eff); + g1 = c4f*Eff*rexp; + g1 = g1 + (0.5*dc4f/stM.bff)*(rexp - 1.0); + g1 = (2.0*stM.aff*g1) + Tfa; + auto Hff = mat_dyad_prod(fl.col(0), fl.col(0), nsd); + S = S + g1*Hff; + + // Fiber-fiber interaction stiffness + g1 = c4f*(1.0 + (2.0*stM.bff*Eff*Eff)); + g1 = (g1 + (2.0*dc4f*Eff))*rexp; + g1 = g1 + (0.5*ddc4f/stM.bff)*(rexp - 1.0); + g1 = 4.0*stM.aff*g1; + CC = CC + g1*ten_dyad_prod(Hff, Hff, nsd); + + // Sheet-sheet interaction stress + additional cross-fiber stress (Tsa) + rexp = exp(stM.bss * Ess * Ess); + double g2 = c4s*Ess*rexp; + g2 = g2 + (0.5*dc4s/stM.bss)*(rexp - 1.0); + g2 = 2.0*stM.ass*g2 + Tsa; + auto Hss = mat_dyad_prod(fl.col(1), fl.col(1), nsd); + S = S + g2*Hss; + + // Sheet-sheet interaction stiffness + g2 = c4s*(1.0 + (2.0*stM.bss*Ess*Ess)); + g2 = (g2 + (2.0*dc4s*Ess))*rexp; + g2 = g2 + (0.5*ddc4s/stM.bss)*(rexp - 1.0); + g2 = 4.0*stM.ass*g2; + CC = CC + g2*ten_dyad_prod(Hss, Hss, nsd); + } break; + default: throw std::runtime_error("Undefined material constitutive model."); } @@ -532,6 +651,7 @@ void get_pk2cc_dev(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& // Fiber-reinforced stress double Tfa = 0.0; get_fib_stress(com_mod, cep_mod, stM.Tf, Tfa); + double Tsa = Tfa*stM.Tf.eta_s; // Electromechanics coupling - active stress if (cep_mod.cem.aStress) { @@ -756,39 +876,68 @@ void get_pk2cc_dev(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& double Ess = Inv6 - 1.0; double Efs = Inv8; + // Smoothed Heaviside function: 1 / (1 + exp(-kx)) = 1 - 1 / (1 + exp(kx)) + double k = stM.khs; + double one_over_exp_plus_one_f = 1.0 / (exp(k * Eff) + 1.0); + double one_over_exp_plus_one_s = 1.0 / (exp(k * Ess) + 1.0); + double c4f = 1.0 - one_over_exp_plus_one_f; + double c4s = 1.0 - one_over_exp_plus_one_s; + + // Approx. derivative of smoothed heaviside function + //double dc4f = 0.25*stM.khs*exp(-stM.khs*abs(Eff)); + //double dc4s = 0.25*stM.khs*exp(-stM.khs*abs(Ess)); + + // Exact first derivative of smoothed heaviside function (from Wolfram Alpha) + double dc4f = k * (one_over_exp_plus_one_f - pow(one_over_exp_plus_one_f,2)); + double dc4s = k * (one_over_exp_plus_one_s - pow(one_over_exp_plus_one_s,2)); + + // Exact second derivative of smoothed heaviside function (from Wolfram Alpha) + double ddc4f = pow(k,2) * (-one_over_exp_plus_one_f + 3.0*pow(one_over_exp_plus_one_f,2) - 2.0*pow(one_over_exp_plus_one_f,3)); + double ddc4s = pow(k,2) * (-one_over_exp_plus_one_s + 3.0*pow(one_over_exp_plus_one_s,2) - 2.0*pow(one_over_exp_plus_one_s,3)); + + // Isotropic + fiber-sheet interaction stress double g1 = stM.a * exp(stM.b*(Inv1-3.0)); - double g2 = 2.0 * stM.afs * Efs * exp(stM.bfs*Efs*Efs); + double g2 = 2.0 * stM.afs * exp(stM.bfs*Efs*Efs); auto Hfs = mat_symm_prod(fl.col(0), fl.col(1), nsd); - auto Sb = g1*IDm + g2*Hfs; + auto Sb = g1*IDm + g2*Efs*Hfs; - Efs = Efs * Efs; - g1 = 2.0*J4d*stM.b*g1; - g2 = 4.0*J4d*stM.afs*(1.0 + 2.0*stM.bfs*Efs)* exp(stM.bfs*Efs); + // Isotropic + fiber-sheet interaction stiffness + g1 = g1 * 2.0 * J4d * stM.b; + g2 = g2 * 2.0 * J4d * (1.0 + 2.0*stM.bfs*Efs*Efs); auto CCb = g1 * ten_dyad_prod(IDm, IDm, nsd) + g2 * ten_dyad_prod(Hfs, Hfs, nsd); - // Fiber reinforcement/active stress - // - if (Eff > 0.0) { - g1 = Tfa; - g1 = g1 + 2.0 * stM.aff * Eff * exp(stM.bff*Eff*Eff); - auto Hff = mat_dyad_prod(fl.col(0), fl.col(0), nsd); - Sb = Sb + g1*Hff; - - Eff = Eff * Eff; - g1 = 4.0*J4d*stM.aff*(1.0 + 2.0*stM.bff*Eff)*exp(stM.bff*Eff); - CCb = CCb + g1*ten_dyad_prod(Hff, Hff, nsd); - } - - if (Ess > 0.0) { - g2 = 2.0 * stM.ass * Ess * exp(stM.bss*Ess*Ess); - auto Hss = mat_dyad_prod(fl.col(1), fl.col(1), nsd); - Sb = Sb + g2*Hss; - - Ess = Ess * Ess; - g2 = 4.0*J4d*stM.ass*(1.0 + 2.0*stM.bss*Ess)*exp(stM.bss*Ess); - CCb = CCb + g2*ten_dyad_prod(Hss, Hss, nsd); - } + // Fiber-fiber interaction stress + additional reinforcement (Tfa) + double rexp = exp(stM.bff*Eff*Eff); + g1 = c4f * Eff * rexp; + g1 = g1 + (0.5*dc4f/stM.bff) * (rexp - 1.0); + g1 = 2.0 * stM.aff * g1 + Tfa; + auto Hff = mat_dyad_prod(fl.col(0), fl.col(0), nsd); + Sb = Sb + g1*Hff; + + // Fiber-fiber interaction stiffness + g1 = c4f * (1.0 + 2.0*stM.bff*Eff*Eff); + g1 = (g1 + 2.0*dc4f*Eff) * rexp; + g1 = g1 + (0.5*ddc4f/stM.bff)*(rexp - 1.0); + g1 = 4.0 * J4d * stM.aff * g1; + CCb = CCb + g1*ten_dyad_prod(Hff, Hff, nsd); + + // Sheet-sheet interaction stress + additional cross-fiber stress (Tsa) + rexp = exp(stM.bss*Ess*Ess); + g2 = c4s * Ess * rexp; + g2 = g2 + (0.5*dc4s/stM.bss) * (rexp - 1.0); + g2 = 2.0 * stM.ass * g2 + Tsa; + auto Hss = mat_dyad_prod(fl.col(1), fl.col(1), nsd); + Sb = Sb + g2*Hss; + + // Sheet-sheet interaction stiffness + g2 = c4s * (1.0 + 2.0*stM.bss*Ess*Ess); + g2 = (g2 + 2.0*dc4s*Ess) * rexp; + g2 = g2 + (0.5*ddc4s/stM.bss)*(rexp - 1.0); + g2 = 4.0 * J4d * stM.ass * g2; + CCb = CCb + g2*ten_dyad_prod(Hss, Hss, nsd); + + // Isochoric 2nd-Piola-Kirchoff stress and stiffness tensors double r1 = J2d*mat_ddot(C, Sb, nsd) / nd; S = J2d*Sb - r1*Ci; @@ -810,6 +959,100 @@ void get_pk2cc_dev(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& } break; + // HO (Holzapfel-Ogden)-MA model for myocardium with full invariants for the anisotropy terms (modified-anisotropy) + case ConstitutiveModelType::stIso_HO_ma: { + if (nfd != 2) { + //err = "Min fiber directions not defined for Holzapfel material model (2)" + } + double Inv4 = norm(fl.col(0), mat_mul(C, fl.col(0))); + double Inv6 = norm(fl.col(1), mat_mul(C, fl.col(1))); + double Inv8 = norm(fl.col(0), mat_mul(C, fl.col(1))); + + //double fds = norm(fl.col(0),fl.col(1)); + double Eff = Inv4 - 1.0; + double Ess = Inv6 - 1.0; + double Efs = Inv8; + + // Smoothed Heaviside function: 1 / (1 + exp(-kx)) = 1 - 1 / (1 + exp(kx)) + double k = stM.khs; + double one_over_exp_plus_one_f = 1.0 / (exp(k * Eff) + 1.0); + double one_over_exp_plus_one_s = 1.0 / (exp(k * Ess) + 1.0); + double c4f = 1.0 - one_over_exp_plus_one_f; + double c4s = 1.0 - one_over_exp_plus_one_s; + + // Approx. derivative of smoothed heaviside function + //double dc4f = 0.25*stM.khs*exp(-stM.khs*abs(Eff)); + //double dc4s = 0.25*stM.khs*exp(-stM.khs*abs(Ess)); + + // Exact first derivative of smoothed heaviside function (from Wolfram Alpha) + double dc4f = k * (one_over_exp_plus_one_f - pow(one_over_exp_plus_one_f,2)); + double dc4s = k * (one_over_exp_plus_one_s - pow(one_over_exp_plus_one_s,2)); + + // Exact second derivative of smoothed heaviside function (from Wolfram Alpha) + double ddc4f = pow(k,2) * (-one_over_exp_plus_one_f + 3.0*pow(one_over_exp_plus_one_f,2) - 2.0*pow(one_over_exp_plus_one_f,3)); + double ddc4s = pow(k,2) * (-one_over_exp_plus_one_s + 3.0*pow(one_over_exp_plus_one_s,2) - 2.0*pow(one_over_exp_plus_one_s,3)); + + // Isochoric stress and stiffness + double g1 = stM.a * exp(stM.b*(Inv1-3.0)); + auto Sb = g1*IDm; + double r1 = J2d/nd*mat_fun::mat_ddot(C, Sb, nsd); + + g1 = g1*2.0*J4d*stM.b; + auto CCb = g1 * ten_dyad_prod(IDm, IDm, nsd); + + // Add isochoric stress and stiffness contribution + 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); + CC = ten_transpose(CC, nsd); + CC = ten_ddot(PP, CC, nsd); + CC = CC + 2.0*r1 * (ten_symm_prod(Ci, Ci, nsd) - + 1.0/nd * ten_dyad_prod(Ci, Ci, nsd)) + - 2.0/nd * (ten_dyad_prod(Ci, S, nsd) + ten_dyad_prod(S, Ci, nsd)); + + // Now add aniostropic components + // Fiber-sheet interaction terms + g1 = 2.0 * stM.afs * exp(stM.bfs*Efs*Efs); + auto Hfs = mat_symm_prod(fl.col(0), fl.col(1), nsd); + S = S + g1*Efs*Hfs; + + g1 = g1 * 2.0*(1.0 + 2.0*stM.bfs*Efs*Efs); + + CC = CC + g1*ten_dyad_prod(Hfs, Hfs, nsd); + + // Fiber-fiber interaction stress + additional reinforcement (Tfa) + double rexp = exp(stM.bff * Eff * Eff); + g1 = c4f*Eff*rexp; + g1 = g1 + (0.5*dc4f/stM.bff)*(rexp - 1.0); + g1 = (2.0*stM.aff*g1) + Tfa; + auto Hff = mat_dyad_prod(fl.col(0), fl.col(0), nsd); + S = S + g1*Hff; + + // Fiber-fiber interaction stiffness + g1 = c4f*(1.0 + (2.0*stM.bff*Eff*Eff)); + g1 = (g1 + (2.0*dc4f*Eff))*rexp; + g1 = g1 + (0.5*ddc4f/stM.bff)*(rexp - 1.0); + g1 = 4.0*stM.aff*g1; + CC = CC + g1*ten_dyad_prod(Hff, Hff, nsd); + + // Sheet-sheet interaction stress + additional cross-fiber stress (Tsa) + rexp = exp(stM.bss * Ess * Ess); + double g2 = c4s*Ess*rexp; + g2 = g2 + (0.5*dc4s/stM.bss)*(rexp - 1.0); + g2 = 2.0*stM.ass*g2 + Tsa; + auto Hss = mat_dyad_prod(fl.col(1), fl.col(1), nsd); + S = S + g2*Hss; + + // Sheet-sheet interaction stiffness + g2 = c4s*(1.0 + (2.0*stM.bss*Ess*Ess)); + g2 = (g2 + (2.0*dc4s*Ess))*rexp; + g2 = g2 + (0.5*ddc4s/stM.bss)*(rexp - 1.0); + g2 = 4.0*stM.ass*g2; + CC = CC + g2*ten_dyad_prod(Hss, Hss, nsd); + } break; + + default: throw std::runtime_error("Undefined isochoric material constitutive model."); } diff --git a/Code/Source/svFSI/mat_models_carray.h b/Code/Source/svFSI/mat_models_carray.h index 85f80d31..3794e998 100644 --- a/Code/Source/svFSI/mat_models_carray.h +++ b/Code/Source/svFSI/mat_models_carray.h @@ -222,6 +222,7 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn // Fiber-reinforced stress double Tfa = 0.0; mat_models::get_fib_stress(com_mod, cep_mod, stM.Tf, Tfa); + double Tsa = Tfa*stM.Tf.eta_s; // Electromechanics coupling - active stress if (cep_mod.cem.aStress) { @@ -906,41 +907,57 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn //err = "Min fiber directions not defined for Holzapfel material model (2)" } + // Compute fiber-based isochoric invariants double C_fl[N]; mat_fun_carray::mat_mul(C, fl.rcol(0), C_fl); double Inv4 = J2d * mat_fun_carray::norm(fl.rcol(0), C_fl); - //double Inv4 = J2d*utils::norm(fl.col(0), mat_mul(C, fl.col(0))); mat_fun_carray::mat_mul(C, fl.rcol(1), C_fl); double Inv6 = J2d * mat_fun_carray::norm(fl.rcol(1), C_fl); - //double Inv6 = J2d*utils::norm(fl.col(1), mat_mul(C, fl.col(1))); - mat_fun_carray::mat_mul(C, fl.rcol(0), C_fl); - double Inv8 = J2d * mat_fun_carray::norm(fl.rcol(1), C_fl); - //double Inv8 = J2d*utils::norm(fl.col(0), mat_mul(C, fl.col(1))); + mat_fun_carray::mat_mul(C, fl.rcol(1), C_fl); + double Inv8 = J2d * mat_fun_carray::norm(fl.rcol(0), C_fl); double Eff = Inv4 - 1.0; double Ess = Inv6 - 1.0; double Efs = Inv8; + // Smoothed Heaviside function: 1 / (1 + exp(-kx)) = 1 - 1 / (1 + exp(kx)) + double k = stM.khs; + double one_over_exp_plus_one_f = 1.0 / (exp(k * Eff) + 1.0); + double one_over_exp_plus_one_s = 1.0 / (exp(k * Ess) + 1.0); + double c4f = 1.0 - one_over_exp_plus_one_f; + double c4s = 1.0 - one_over_exp_plus_one_s; + + // Approx. derivative of smoothed heaviside function + // double dc4f = 0.25*stM.khs*exp(-stM.khs*abs(Eff)); + // double dc4s = 0.25*stM.khs*exp(-stM.khs*abs(Ess)); + + // Exact first derivative of smoothed heaviside function (from Wolfram Alpha) + double dc4f = k * (one_over_exp_plus_one_f - pow(one_over_exp_plus_one_f,2)); + double dc4s = k * (one_over_exp_plus_one_s - pow(one_over_exp_plus_one_s,2)); + + // Exact second derivative of smoothed heaviside function (from Wolfram Alpha) + double ddc4f = pow(k,2) * (-one_over_exp_plus_one_f + 3.0*pow(one_over_exp_plus_one_f,2) - 2.0*pow(one_over_exp_plus_one_f,3)); + double ddc4s = pow(k,2) * (-one_over_exp_plus_one_s + 3.0*pow(one_over_exp_plus_one_s,2) - 2.0*pow(one_over_exp_plus_one_s,3)); + + // Isotropic + fiber-sheet interaction stress double g1 = stM.a * exp(stM.b*(Inv1-3.0)); - double g2 = 2.0 * stM.afs * Efs * exp(stM.bfs*Efs*Efs); + double g2 = 2.0 * stM.afs * exp(stM.bfs*Efs*Efs); double Hfs[N][N]; mat_fun_carray::mat_symm_prod(fl.rcol(0), fl.rcol(1), Hfs); - //auto Hfs = mat_symm_prod(fl.col(0), fl.col(1), nsd); double Sb[N][N]; for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { - Sb[i][j] = g1*Idm[i][j] + g2*Hfs[i][j]; + Sb[i][j] = g1*Idm[i][j] + g2*Efs*Hfs[i][j]; } } - //auto Sb = g1*Idm + g2*Hfs; - Efs = Efs * Efs; - g1 = 2.0*J4d*stM.b*g1; - g2 = 4.0*J4d*stM.afs*(1.0 + 2.0*stM.bfs*Efs)* exp(stM.bfs*Efs); + // Isotropic + fiber-sheet interaction stiffness + g1 = g1 * 2.0 * J4d * stM.b; + g2 = g2 * 2.0 * J4d * (1.0 + 2.0*stM.bfs*Efs*Efs); CArray4 Idm_prod; mat_fun_carray::ten_dyad_prod(Idm, Idm, Idm_prod); @@ -959,84 +976,79 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn } } } - //auto CCb = g1 * ten_dyad_prod(Idm, Idm, nsd) + g2 * ten_dyad_prod(Hfs, Hfs, nsd); - - // Fiber reinforcement/active stress - if (Eff > 0.0) { - g1 = Tfa; - g1 = g1 + 2.0 * stM.aff * Eff * exp(stM.bff*Eff*Eff); - double Hff[N][N]; - mat_fun_carray::mat_dyad_prod(fl.col(0), fl.col(0), Hff); - //auto Hff = mat_dyad_prod(fl.col(0), fl.col(0), nsd); - - for (int i = 0; i < N; i++) { + // Fiber-fiber interaction stress + additional reinforcement (Tfa) + double rexp = exp(stM.bff*Eff*Eff); + g1 = c4f * Eff * rexp; + g1 = g1 + (0.5*dc4f/stM.bff) * (rexp - 1.0); + g1 = 2.0 * stM.aff * g1 + Tfa; + double Hff[N][N]; + mat_fun_carray::mat_dyad_prod(fl.col(0), fl.col(0), Hff); + for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { Sb[i][j] += g1 * Hff[i][j]; } } - //Sb = Sb + g1*Hff; - - Eff = Eff * Eff; - g1 = 4.0*J4d*stM.aff*(1.0 + 2.0*stM.bff*Eff)*exp(stM.bff*Eff); - CArray4 Hff_prod; - mat_fun_carray::ten_dyad_prod(Hff, Hff, Hff_prod); + // Fiber-fiber interaction stiffness + g1 = c4f * (1.0 + 2.0*stM.bff*Eff*Eff); + g1 = (g1 + 2.0*dc4f*Eff) * rexp; + g1 = g1 + (0.5*ddc4f/stM.bff)*(rexp - 1.0); + g1 = 4.0 * J4d * stM.aff * g1; + CArray4 Hff_prod; + mat_fun_carray::ten_dyad_prod(Hff, Hff, Hff_prod); - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - for (int k = 0; k < N; k++) { - for (int l = 0; l < N; l++) { - CCb[i][j][k][l] += g1 * Hff_prod[i][j][k][l]; - } + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + CCb[i][j][k][l] += g1 * Hff_prod[i][j][k][l]; } } } - //CCb = CCb + g1*ten_dyad_prod(Hff, Hff, nsd); } - if (Ess > 0.0) { - g2 = 2.0 * stM.ass * Ess * exp(stM.bss*Ess*Ess); - - double Hss[N][N]; - mat_fun_carray::mat_dyad_prod(fl.col(1), fl.col(1), Hss); - //auto Hss = mat_dyad_prod(fl.col(1), fl.col(1), nsd); + // Sheet-sheet interaction stress + additional cross-fiber stress + rexp = exp(stM.bss*Ess*Ess); + g2 = c4s * Ess * rexp; + g2 = g2 + (0.5*dc4s/stM.bss) * (rexp - 1.0); + g2 = 2.0 * stM.ass * g2 + Tsa; + double Hss[N][N]; + mat_fun_carray::mat_dyad_prod(fl.col(1), fl.col(1), Hss); - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - Sb[i][j] += g2 * Hss[i][j]; - } + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + Sb[i][j] += g2 * Hss[i][j]; } - //Sb = Sb + g2*Hss; + } - Ess = Ess * Ess; - g2 = 4.0 * J4d * stM.ass * (1.0 + 2.0*stM.bss*Ess) * exp(stM.bss * Ess); + // Sheet-sheet interaction stiffness + g2 = c4s * (1.0 + 2.0 * stM.bss * Ess * Ess); + g2 = (g2 + 2.0*dc4s*Ess) * rexp; + g2 = g2 + (0.5*ddc4s/stM.bss)*(rexp - 1.0); + g2 = 4.0 * J4d * stM.ass * g2; - CArray4 Hss_prod; - mat_fun_carray::ten_dyad_prod(Hss, Hss, Hss_prod); + CArray4 Hss_prod; + mat_fun_carray::ten_dyad_prod(Hss, Hss, Hss_prod); - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - for (int k = 0; k < N; k++) { - for (int l = 0; l < N; l++) { - CCb[i][j][k][l] += g2 * Hss_prod[i][j][k][l]; - } + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + CCb[i][j][k][l] += g2 * Hss_prod[i][j][k][l]; } } } - //CCb = CCb + g2*ten_dyad_prod(Hss, Hss, nsd); - } + // Isochoric 2nd-Piola-Kirchhoff stress and stiffness tensors double r1 = J2d * mat_fun_carray::mat_ddot(C, Sb) / nd; - //double r1 = J2d*mat_ddot(C, Sb, nsd) / nd; for (int i = 0; i < nsd; i++) { for (int j = 0; j < nsd; j++) { S[i][j] = J2d*Sb[i][j] - r1*Ci[i][j]; } } - //S = J2d*Sb - r1*Ci; CArray4 Ci_C_prod; mat_fun_carray::ten_dyad_prod(Ci, C, Ci_C_prod); @@ -1051,17 +1063,13 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn } } } - //auto PP = ten_ids(nsd) - (1.0/nd) * ten_dyad_prod(Ci, C, nsd); mat_fun_carray::ten_ddot(CCb, PP, CC); - //CC = ten_ddot(CCb, PP, nsd); CArray4 CC_t; mat_fun_carray::ten_transpose(CC, CC_t); - //CC = ten_transpose(CC, nsd); mat_fun_carray::ten_ddot(PP, CC_t, CC); - //CC = ten_ddot(PP, CC, nsd); CArray4 Ci_S_prod; mat_fun_carray::ten_dyad_prod(Ci, S, Ci_S_prod); @@ -1078,14 +1086,13 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn } } } - // CC = CC - (2.0/nd) * ( ten_dyad_prod(Ci, S, nsd) + ten_dyad_prod(S, Ci, nsd) ); + // Add pressure contribution for (int i = 0; i < nsd; i++) { for (int j = 0; j < nsd; j++) { S[i][j] += p * J * Ci[i][j]; } } - //S = S + p*J*Ci; CArray4 Ci_sym_prod; mat_fun_carray::ten_symm_prod(Ci, Ci, Ci_sym_prod); @@ -1102,7 +1109,6 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn } } } - //CC = CC + 2.0*(r1 - p*J) * ten_symm_prod(Ci, Ci, nsd) + (pl*J - 2.0*r1/nd) * ten_dyad_prod(Ci, Ci, nsd); if (cep_mod.cem.aStrain) { double S_prod[N][N]; @@ -1127,6 +1133,223 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn } break; +// HO (Holzapfel-Ogden)-MA model for myocardium with full invariants for the anisotropy terms (modified-anisotropy) + case ConstitutiveModelType::stIso_HO_ma: { + if (nfd != 2) { + //err = "Min fiber directions not defined for Holzapfel material model (2)" + } + + // Compute fiber-based full invariants (not isochoric) + double C_fl[N]; + mat_fun_carray::mat_mul(C, fl.rcol(0), C_fl); + double Inv4 = mat_fun_carray::norm(fl.rcol(0), C_fl); + + mat_fun_carray::mat_mul(C, fl.rcol(1), C_fl); + double Inv6 = mat_fun_carray::norm(fl.rcol(1), C_fl); + + mat_fun_carray::mat_mul(C, fl.rcol(1), C_fl); + double Inv8 = mat_fun_carray::norm(fl.rcol(0), C_fl); + + double Eff = Inv4 - 1.0; + double Ess = Inv6 - 1.0; + double Efs = Inv8; + + // Smoothed Heaviside function: 1 / (1 + exp(-kx)) = 1 - 1 / (1 + exp(kx)) + double k = stM.khs; + double one_over_exp_plus_one_f = 1.0 / (exp(k * Eff) + 1.0); + double one_over_exp_plus_one_s = 1.0 / (exp(k * Ess) + 1.0); + double c4f = 1.0 - one_over_exp_plus_one_f; + double c4s = 1.0 - one_over_exp_plus_one_s; + + // Approx. derivative of smoothed heaviside function + // double dc4f = 0.25*stM.khs*exp(-stM.khs*abs(Eff)); + // double dc4s = 0.25*stM.khs*exp(-stM.khs*abs(Ess)); + + // Exact first derivative of smoothed heaviside function (from Wolfram Alpha) + double dc4f = k * (one_over_exp_plus_one_f - pow(one_over_exp_plus_one_f,2)); + double dc4s = k * (one_over_exp_plus_one_s - pow(one_over_exp_plus_one_s,2)); + + // Exact second derivative of smoothed heaviside function (from Wolfram Alpha) + double ddc4f = pow(k,2) * (-one_over_exp_plus_one_f + 3.0*pow(one_over_exp_plus_one_f,2) - 2.0*pow(one_over_exp_plus_one_f,3)); + double ddc4s = pow(k,2) * (-one_over_exp_plus_one_s + 3.0*pow(one_over_exp_plus_one_s,2) - 2.0*pow(one_over_exp_plus_one_s,3)); + + // Isochoric stress and stiffness + double g1 = stM.a * exp(stM.b*(Inv1-3.0)); + + double Sb[N][N]; + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + Sb[i][j] = g1*Idm[i][j]; + } + } + double r1 = J2d/nd*mat_fun_carray::mat_ddot(C,Sb); + g1 = g1*2.0*J4d*stM.b; + + CArray4 Idm_prod; + mat_fun_carray::ten_dyad_prod(Idm, Idm, Idm_prod); + + CArray4 CCb; + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + CCb[i][j][k][l] = g1 * Idm_prod[i][j][k][l]; + } + } + } + } + + // Add isochoric stress and stiffness contribution + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + S[i][j] = J2d*Sb[i][j] - r1*Ci[i][j]; + } + } + CArray4 Ci_C_prod; + mat_fun_carray::ten_dyad_prod(Ci, C, Ci_C_prod); + double PP[N][N][N][N]; + + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + PP[i][j][k][l] = Ids[i][j][k][l] - (1.0/nd) * Ci_C_prod[i][j][k][l]; + } + } + } + } + + mat_fun_carray::ten_ddot(CCb, PP, CC); + + CArray4 CC_t; + mat_fun_carray::ten_transpose(CC, CC_t); + + mat_fun_carray::ten_ddot(PP, CC_t, CC); + + CArray4 Ci_S_prod; + mat_fun_carray::ten_dyad_prod(Ci, S, Ci_S_prod); + + CArray4 S_Ci_prod; + mat_fun_carray::ten_dyad_prod(S, Ci, S_Ci_prod); + + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + CC[i][j][k][l] -= (2.0/nd) * (Ci_S_prod[i][j][k][l] + S_Ci_prod[i][j][k][l]); + } + } + } + } + + // Add pressure contribution to stress and stiffness + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + S[i][j] += p*J*Ci[i][j]; + } + } + + CArray4 Ci_Ci_prod; + mat_fun_carray::ten_dyad_prod(Ci, Ci, Ci_Ci_prod); + CArray4 Ci_Ci_symprod; + mat_fun_carray::ten_symm_prod(Ci, Ci, Ci_Ci_symprod); + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + CC[i][j][k][l] += 2.0*(r1 - p*J)*Ci_Ci_symprod[i][j][k][l] + (pl*J - 2.0*r1/nd)*Ci_Ci_prod[i][j][k][l]; + } + } + } + } + + // Now that both isochoric and volumetric components were added, anisotropic components need to be added + + // Fiber-sheet interaction terms + g1 = 2.0 * stM.afs * exp(stM.bfs*Efs*Efs); + double Hfs[N][N]; + mat_fun_carray::mat_symm_prod(fl.rcol(0), fl.rcol(1), Hfs); + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + S[i][j] += g1*Efs*Hfs[i][j]; + } + } + + g1 = g1 * 2.0*(1.0 + 2.0*stM.bfs*Efs*Efs); + + CArray4 Hfs_Hfs_prod; + mat_fun_carray::ten_dyad_prod(Hfs, Hfs, Hfs_Hfs_prod); + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + CC[i][j][k][l] += g1*Hfs_Hfs_prod[i][j][k][l]; + } + } + } + } + + // Fiber-fiber interaction stress + additional reinforcement (Tfa) + double rexp = exp(stM.bff * Eff * Eff); + g1 = c4f*Eff*rexp; + g1 = g1 + (0.5*dc4f/stM.bff)*(rexp - 1.0); + g1 = (2.0*stM.aff*g1) + Tfa; + double Hff[N][N]; + mat_fun_carray::mat_dyad_prod(fl.rcol(0), fl.rcol(0), Hff); + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + S[i][j] += g1*Hff[i][j]; + } + } + + // Fiber-fiber interaction stiffness + g1 = c4f*(1.0 + (2.0*stM.bff*Eff*Eff)); + g1 = (g1 + (2.0*dc4f*Eff))*rexp; + g1 = g1 + (0.5*ddc4f/stM.bff)*(rexp - 1.0); + g1 = 4.0*stM.aff*g1; + CArray4 Hff_Hff_prod; + mat_fun_carray::ten_dyad_prod(Hff, Hff, Hff_Hff_prod); + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + CC[i][j][k][l] += g1*Hff_Hff_prod[i][j][k][l]; + } + } + } + } + + // Sheet-sheet interaction stress + additional cross-fiber stress + rexp = exp(stM.bss * Ess * Ess); + double g2 = c4s*Ess*rexp; + g2 = g2 + (0.5*dc4s/stM.bss)*(rexp - 1.0); + g2 = 2.0*stM.ass*g2 + Tsa; + double Hss[N][N]; + mat_fun_carray::mat_dyad_prod(fl.rcol(1), fl.rcol(1), Hss); + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + S[i][j] += g2*Hss[i][j]; + } + } + + // Sheet-sheet interaction stiffness + g2 = c4s*(1.0 + (2.0*stM.bss*Ess*Ess)); + g2 = (g2 + (2.0*dc4s*Ess))*rexp; + g2 = g2 + (0.5*ddc4s/stM.bss)*(rexp - 1.0); + g2 = 4.0*stM.ass*g2; + CArray4 Hss_Hss_prod; + mat_fun_carray::ten_dyad_prod(Hss,Hss,Hss_Hss_prod); + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + CC[i][j][k][l] += g2*Hss_Hss_prod[i][j][k][l]; + } + } + } + } + } break; + default: throw std::runtime_error("Undefined material constitutive model."); } diff --git a/Code/Source/svFSI/post.cpp b/Code/Source/svFSI/post.cpp index ba0f1f7d..68f3f643 100644 --- a/Code/Source/svFSI/post.cpp +++ b/Code/Source/svFSI/post.cpp @@ -36,6 +36,8 @@ #include "initialize.h" #include "mat_fun.h" #include "mat_models.h" +#include "mat_fun_carray.h" +#include "mat_models_carray.h" #include "nn.h" #include "shells.h" #include "utils.h" diff --git a/Code/Source/svFSI/read_files.cpp b/Code/Source/svFSI/read_files.cpp index 52ae6809..a60d848b 100644 --- a/Code/Source/svFSI/read_files.cpp +++ b/Code/Source/svFSI/read_files.cpp @@ -1803,7 +1803,7 @@ void read_files(Simulation* simulation, const std::string& file_name) if ((dmn.phys != EquationType::phys_ustruct) && (dmn.phys != EquationType::phys_struct)) { continue; } - if (dmn.stM.isoType != ConstitutiveModelType::stIso_HO) { + if ((dmn.stM.isoType != ConstitutiveModelType::stIso_HO)) { throw std::runtime_error("Active strain is allowed with Holzapfel-Ogden passive constitutive model only"); } } diff --git a/Code/Source/svFSI/set_material_props.h b/Code/Source/svFSI/set_material_props.h index 70d8eb22..ccd5bab2 100644 --- a/Code/Source/svFSI/set_material_props.h +++ b/Code/Source/svFSI/set_material_props.h @@ -155,6 +155,26 @@ SeMaterialPropertiesMapType set_material_props = { lDmn.stM.khs = params.k.value(); } }, +//---------------------------// +// stIso_HO_ma // +//---------------------------// +// +{consts::ConstitutiveModelType::stIso_HO_ma, [](DomainParameters* domain_params, double mu, double kap, double lam, + dmnType& lDmn) -> void +{ + lDmn.stM.isoType = consts::ConstitutiveModelType::stIso_HO_ma; + auto& params = domain_params->constitutive_model.holzapfel; + lDmn.stM.a = params.a.value(); + lDmn.stM.b = params.b.value(); + lDmn.stM.aff = params.a4f.value(); + lDmn.stM.bff = params.b4f.value(); + lDmn.stM.ass = params.a4s.value(); + lDmn.stM.bss = params.b4s.value(); + lDmn.stM.afs = params.afs.value(); + lDmn.stM.bfs = params.bfs.value(); + lDmn.stM.khs = params.k.value(); +} }, + //---------------------------// // stIso_LS // //---------------------------// diff --git a/tests/cases/struct/LV_Holzapfel_passive/endo_pressure.dat b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/endo_pressure.dat similarity index 100% rename from tests/cases/struct/LV_Holzapfel_passive/endo_pressure.dat rename to tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/endo_pressure.dat diff --git a/tests/cases/struct/LV_Holzapfel_passive/mesh/fibersLongCells.vtu b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersLongCells.vtu similarity index 100% rename from tests/cases/struct/LV_Holzapfel_passive/mesh/fibersLongCells.vtu rename to tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersLongCells.vtu diff --git a/tests/cases/struct/LV_Holzapfel_passive/mesh/fibersSheetCells.vtu b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersSheetCells.vtu similarity index 100% rename from tests/cases/struct/LV_Holzapfel_passive/mesh/fibersSheetCells.vtu rename to tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersSheetCells.vtu diff --git a/tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-complete.exterior.vtp b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.exterior.vtp similarity index 100% rename from tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-complete.exterior.vtp rename to tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.exterior.vtp diff --git a/tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-complete.mesh.vtu b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.mesh.vtu similarity index 100% rename from tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-complete.mesh.vtu rename to tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.mesh.vtu diff --git a/tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-surfaces/endo.vtp b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/endo.vtp similarity index 100% rename from tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-surfaces/endo.vtp rename to tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/endo.vtp diff --git a/tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-surfaces/epi.vtp b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/epi.vtp similarity index 100% rename from tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-surfaces/epi.vtp rename to tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/epi.vtp diff --git a/tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-surfaces/top.vtp b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/top.vtp similarity index 100% rename from tests/cases/struct/LV_Holzapfel_passive/mesh/mesh-surfaces/top.vtp rename to tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/top.vtp diff --git a/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/result_001.vtu b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/result_001.vtu new file mode 100644 index 00000000..001d74ad --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:282b0f90edcdcd4b7e47032c2aa47c4ab47aaf66717050ed1001d6483e36c741 +size 249794 diff --git a/tests/cases/struct/LV_Holzapfel_passive/svFSIplus.xml b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/svFSIplus.xml similarity index 86% rename from tests/cases/struct/LV_Holzapfel_passive/svFSIplus.xml rename to tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/svFSIplus.xml index 88e6694d..8f1c7955 100644 --- a/tests/cases/struct/LV_Holzapfel_passive/svFSIplus.xml +++ b/tests/cases/struct/LV_HolzapfelOgdenModifiedAnisotropy_passive/svFSIplus.xml @@ -5,7 +5,7 @@ 0 3 1 - 1e-3 + 1e-2 0.50 STOP_SIM @@ -52,14 +52,14 @@ true 3 - 4 + 20 1e-10 1.0 1.0e6 0.483333 - + 590.0 8.023 184720.0 @@ -72,6 +72,7 @@ ST91 + 1.0e7 true @@ -88,22 +89,14 @@ fsils - 1e-12 + 1e-16 1000 50 - - Robin - 1.0e7 - 5.0e2 - 1 - - - Robin - 1.0e4 - 5.0e2 + Dirichlet + 0.0 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/endo_pressure.dat b/tests/cases/struct/LV_HolzapfelOgden_passive/endo_pressure.dat new file mode 100755 index 00000000..9762d51b --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/endo_pressure.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3cf8ccb1b43920269cf72e84d8a5332c053a1a24865ce845284987aa1a3872ad +size 20865 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/fibersLongCells.vtu b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/fibersLongCells.vtu new file mode 100644 index 00000000..61b9d063 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/fibersLongCells.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:de695cbe805a262926972c8d2e59be6d82f43bbb0e10e6e0bb4780832b32cc0d +size 93997 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/fibersSheetCells.vtu b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/fibersSheetCells.vtu new file mode 100644 index 00000000..00c92f85 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/fibersSheetCells.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:00ce8af58380bc6f682a5fbab0473818eeeb9f8024b419fbe928d38f5f8e34ec +size 99744 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-complete.exterior.vtp b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-complete.exterior.vtp new file mode 100644 index 00000000..37ab77b3 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:516ca5f67f51c432410c65a1f98d73baca653b0378664607a20eba3fd5e94db4 +size 19832 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-complete.mesh.vtu b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..cf145a1d --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9b7c72acb6c365fae7c33091f8d932c95863d95ffee7fed934a2f832aeba976d +size 27192 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/endo.vtp b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/endo.vtp new file mode 100644 index 00000000..e284a4c6 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/endo.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4028808552b03a2854ce33bf5016243c6b994e62a61928625691b72938f1480d +size 10038 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/epi.vtp b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/epi.vtp new file mode 100644 index 00000000..66272c04 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f125c7001fc4a5ad90313fcff8ef24749ccd054a6b045521dfc3d06169160a09 +size 11982 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/top.vtp b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/top.vtp new file mode 100644 index 00000000..2d90c222 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2050d191642b719f651de6587571a7021eff67801eddd97c0fec397f3e9b62cd +size 3940 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/result_001.vtu b/tests/cases/struct/LV_HolzapfelOgden_passive/result_001.vtu new file mode 100644 index 00000000..72d5b06b --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:98380148f70f0379a869619d69bab3c8916003185bb9ea5301eda2f88bb9aef9 +size 249906 diff --git a/tests/cases/struct/LV_HolzapfelOgden_passive/svFSIplus.xml b/tests/cases/struct/LV_HolzapfelOgden_passive/svFSIplus.xml new file mode 100644 index 00000000..932d005a --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_passive/svFSIplus.xml @@ -0,0 +1,112 @@ + + + + + 0 + 3 + 1 + 1e-2 + 0.50 + STOP_SIM + + 1 + result + + 1 + 1 + + 50 + 0 + + 1 + 1 + 1 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/endo.vtp + + + + mesh/mesh-surfaces/epi.vtp + + + + mesh/mesh-surfaces/top.vtp + + + mesh/fibersLongCells.vtu + mesh/fibersSheetCells.vtu + + 100.0 + + + + + + true + 3 + 20 + 1e-10 + + 1.0 + 1.0e6 + 0.483333 + + + 590.0 + 8.023 + 184720.0 + 16.026 + 24810.0 + 11.12 + 2160.0 + 11.436 + 100.0 + + + ST91 + 1.0e7 + + + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-16 + 1000 + 50 + + + + Dirichlet + 0.0 + + + + Neu + Unsteady + endo_pressure.dat + false + true + + + + + diff --git a/tests/cases/struct/LV_Holzapfel_passive/result_001.vtu b/tests/cases/struct/LV_Holzapfel_passive/result_001.vtu deleted file mode 100644 index 3988c33f..00000000 --- a/tests/cases/struct/LV_Holzapfel_passive/result_001.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:5480109844a69127597cf1ea8b16d6374e0b0fea5b6de06d1a17074df4dc0453 -size 249954 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/endo_pressure.dat b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/endo_pressure.dat new file mode 100755 index 00000000..9762d51b --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/endo_pressure.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3cf8ccb1b43920269cf72e84d8a5332c053a1a24865ce845284987aa1a3872ad +size 20865 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersLongCells.vtu b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersLongCells.vtu new file mode 100644 index 00000000..61b9d063 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersLongCells.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:de695cbe805a262926972c8d2e59be6d82f43bbb0e10e6e0bb4780832b32cc0d +size 93997 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersSheetCells.vtu b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersSheetCells.vtu new file mode 100644 index 00000000..00c92f85 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/fibersSheetCells.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:00ce8af58380bc6f682a5fbab0473818eeeb9f8024b419fbe928d38f5f8e34ec +size 99744 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.exterior.vtp b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.exterior.vtp new file mode 100644 index 00000000..37ab77b3 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:516ca5f67f51c432410c65a1f98d73baca653b0378664607a20eba3fd5e94db4 +size 19832 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.mesh.vtu b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..cf145a1d --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9b7c72acb6c365fae7c33091f8d932c95863d95ffee7fed934a2f832aeba976d +size 27192 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/endo.vtp b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/endo.vtp new file mode 100644 index 00000000..e284a4c6 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/endo.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4028808552b03a2854ce33bf5016243c6b994e62a61928625691b72938f1480d +size 10038 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/epi.vtp b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/epi.vtp new file mode 100644 index 00000000..66272c04 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f125c7001fc4a5ad90313fcff8ef24749ccd054a6b045521dfc3d06169160a09 +size 11982 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/top.vtp b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/top.vtp new file mode 100644 index 00000000..2d90c222 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2050d191642b719f651de6587571a7021eff67801eddd97c0fec397f3e9b62cd +size 3940 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/result_001.vtu b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/result_001.vtu new file mode 100644 index 00000000..317965b7 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fda0ce16348dd0d7e46f8fc3fd98dc68b71ec3d3c8363705656cfbfc7e8a1f71 +size 256530 diff --git a/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/svFSIplus.xml b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/svFSIplus.xml new file mode 100644 index 00000000..101e0f0e --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgdenModifiedAnisotropy_passive/svFSIplus.xml @@ -0,0 +1,113 @@ + + + + + 0 + 3 + 1 + 1e-2 + 0.50 + STOP_SIM + + 1 + result + + 1 + 1 + + 50 + 0 + + 1 + 1 + 1 + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/endo.vtp + + + + mesh/mesh-surfaces/epi.vtp + + + + mesh/mesh-surfaces/top.vtp + + + mesh/fibersLongCells.vtu + mesh/fibersSheetCells.vtu + + 100.0 + + + + + true + 3 + 20 + 1e-10 + + 1.0 + 1.0e6 + 0.483333 + + + 590.0 + 8.023 + 184720.0 + 16.026 + 24810.0 + 11.12 + 2160.0 + 11.436 + 100.0 + + + ST91 + 1.0e7 + + + true + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-16 + 1000 + 50 + + + + Dirichlet + 0.0 + + + + Neu + Unsteady + endo_pressure.dat + false + true + + + + + + + diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/endo_pressure.dat b/tests/cases/ustruct/LV_HolzapfelOgden_passive/endo_pressure.dat new file mode 100755 index 00000000..9762d51b --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/endo_pressure.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3cf8ccb1b43920269cf72e84d8a5332c053a1a24865ce845284987aa1a3872ad +size 20865 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/fibersLongCells.vtu b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/fibersLongCells.vtu new file mode 100644 index 00000000..61b9d063 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/fibersLongCells.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:de695cbe805a262926972c8d2e59be6d82f43bbb0e10e6e0bb4780832b32cc0d +size 93997 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/fibersSheetCells.vtu b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/fibersSheetCells.vtu new file mode 100644 index 00000000..00c92f85 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/fibersSheetCells.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:00ce8af58380bc6f682a5fbab0473818eeeb9f8024b419fbe928d38f5f8e34ec +size 99744 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-complete.exterior.vtp b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-complete.exterior.vtp new file mode 100644 index 00000000..37ab77b3 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:516ca5f67f51c432410c65a1f98d73baca653b0378664607a20eba3fd5e94db4 +size 19832 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-complete.mesh.vtu b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..cf145a1d --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9b7c72acb6c365fae7c33091f8d932c95863d95ffee7fed934a2f832aeba976d +size 27192 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/endo.vtp b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/endo.vtp new file mode 100644 index 00000000..e284a4c6 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/endo.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4028808552b03a2854ce33bf5016243c6b994e62a61928625691b72938f1480d +size 10038 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/epi.vtp b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/epi.vtp new file mode 100644 index 00000000..66272c04 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f125c7001fc4a5ad90313fcff8ef24749ccd054a6b045521dfc3d06169160a09 +size 11982 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/top.vtp b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/top.vtp new file mode 100644 index 00000000..2d90c222 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2050d191642b719f651de6587571a7021eff67801eddd97c0fec397f3e9b62cd +size 3940 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/result_001.vtu b/tests/cases/ustruct/LV_HolzapfelOgden_passive/result_001.vtu new file mode 100644 index 00000000..6c6c35a7 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b062841ceadf91cf9893ee57852a0ed4d1afbe38e6092787a214fe154b7ce309 +size 256538 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_passive/svFSIplus.xml b/tests/cases/ustruct/LV_HolzapfelOgden_passive/svFSIplus.xml new file mode 100644 index 00000000..a634fad5 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_passive/svFSIplus.xml @@ -0,0 +1,113 @@ + + + + + 0 + 3 + 1 + 1e-2 + 0.50 + STOP_SIM + + 1 + result + + 1 + 1 + + 50 + 0 + + 1 + 1 + 1 + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/endo.vtp + + + + mesh/mesh-surfaces/epi.vtp + + + + mesh/mesh-surfaces/top.vtp + + + mesh/fibersLongCells.vtu + mesh/fibersSheetCells.vtu + + 100.0 + + + + + true + 3 + 20 + 1e-10 + + 1.0 + 1.0e6 + 0.483333 + + + 590.0 + 8.023 + 184720.0 + 16.026 + 24810.0 + 11.12 + 2160.0 + 11.436 + 100.0 + + + ST91 + 1.0e7 + + + true + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-16 + 1000 + 50 + + + + Dirichlet + 0.0 + + + + Neu + Unsteady + endo_pressure.dat + false + true + + + + + + + diff --git a/tests/test_struct.py b/tests/test_struct.py index 9b02081a..01cf57d9 100644 --- a/tests/test_struct.py +++ b/tests/test_struct.py @@ -24,10 +24,13 @@ def test_LV_Guccione_passive(n_proc): run_with_reference(base_folder, test_folder, fields, n_proc) -def test_LV_Holzapfel_passive(n_proc): - test_folder = "LV_Holzapfel_passive" +def test_LV_HolzapfelOgden_passive(n_proc): + test_folder = "LV_HolzapfelOgden_passive" run_with_reference(base_folder, test_folder, fields, n_proc) +def test_LV_HolzapfelOgdenModifiedAnisotropy_passive(n_proc): + test_folder = "LV_HolzapfelOgdenModifiedAnisotropy_passive" + run_with_reference(base_folder, test_folder, fields, n_proc) def test_block_compression(n_proc): test_folder = "block_compression" diff --git a/tests/test_ustruct.py b/tests/test_ustruct.py index f3fb9112..5b3365af 100644 --- a/tests/test_ustruct.py +++ b/tests/test_ustruct.py @@ -34,6 +34,14 @@ def test_LV_Guccione_active(n_proc): test_folder = "LV_Guccione_active" run_with_reference(base_folder, test_folder, fields, n_proc) +def test_LV_HolzapfelOgden_passive(n_proc): + test_folder = "LV_HolzapfelOgden_passive" + run_with_reference(base_folder, test_folder, fields, n_proc) + +def test_LV_HolzapfelOgdenModifiedAnisotropy_passive(n_proc): + test_folder = "LV_HolzapfelOgdenModifiedAnisotropy_passive" + run_with_reference(base_folder, test_folder, fields, n_proc) + def test_LV_NeoHookean_passive_genBC(n_proc): test_folder = "LV_NeoHookean_passive_genBC" diff --git a/tests/unitTests/test.cpp b/tests/unitTests/test.cpp index 661a4680..1c0e1248 100644 --- a/tests/unitTests/test.cpp +++ b/tests/unitTests/test.cpp @@ -33,6 +33,84 @@ using namespace mat_fun; using namespace std; +// ============================================================================ +// --------------------------- Test fixture classes --------------------------- +// ============================================================================ + +/** + * @brief Test fixture class containing common setup for all material model tests in this file + * + */ +class MaterialModelTest : public ::testing::Test { +protected: + // Variables common across tests + double F[3][3] = {}; // Deformation gradient + double deformation_perturbation_small = 0.003; // Small perturbation factor + double deformation_perturbation_medium = 0.03; // Medium perturbation factor + double deformation_perturbation_large = 0.3; // Large perturbation factor + int n_F = 50; // Number of deformation gradients F to test for each small, medium, and large perturbation + double rel_tol = 1e-3; // relative tolerance for comparing values + double abs_tol = 1e-11; // absolute tolerance for comparing values + //double delta = 1e-7; // perturbation scaling factor + double delta_max = 1e-4; // maximum perturbation scaling factor + double delta_min = 1e-6; // minimum perturbation scaling factor + int order = 1; // Order of finite difference method + double convergence_order_tol = 0.02; // Tolerance for comparing convergence order with expected value + bool verbose = false; // Show values of S, dE, SdE and dPsi + + // Type alias for a 3x3 std::array + using Array3x3 = std::array, 3>; + + // Vectors to store the 3x3 arrays + std::vector F_small_list; + std::vector F_medium_list; + std::vector F_large_list; + + // Function to convert C array to std::array + Array3x3 convertToStdArray(double F[3][3]) { + Array3x3 stdArray; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + stdArray[i][j] = F[i][j]; + } + } + return stdArray; + } + + // Function to convert std::array to C array + void convertToCArray(Array3x3 stdArray, double F[3][3]) { + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + F[i][j] = stdArray[i][j]; + } + } + } + + void SetUp() override { + double F[3][3]; // Declare the C array + + // Create random deformation gradients for small perturbations + for (int i = 0; i < n_F; i++) { + create_random_perturbed_identity_F(F, deformation_perturbation_small); + F_small_list.push_back(convertToStdArray(F)); + } + + // Create random deformation gradients for medium perturbations + for (int i = 0; i < n_F; i++) { + create_random_perturbed_identity_F(F, deformation_perturbation_medium); + F_medium_list.push_back(convertToStdArray(F)); + } + + // Create random deformation gradients for large perturbations + for (int i = 0; i < n_F; i++) { + create_random_perturbed_identity_F(F, deformation_perturbation_large); + F_large_list.push_back(convertToStdArray(F)); + } + } + + void TearDown() override {} +}; + // ---------------------------------------------------------------------------- // --------------------------- Neo-Hookean Material --------------------------- // ---------------------------------------------------------------------------- @@ -42,16 +120,10 @@ using namespace std; * * This class sets up the necessary parameters and objects for testing the Neo-Hookean material model. */ -class NeoHookeanTest : public ::testing::Test { +class NeoHookeanTest : public MaterialModelTest { protected: - // Variables common across tests + // Material parameters object NeoHookeanParams 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 values - double abs_tol = 1e-11; // 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 TestNeoHookean* TestNH; @@ -59,6 +131,8 @@ class NeoHookeanTest : public ::testing::Test { // Setup method to initialize variables before each test void SetUp() override { + MaterialModelTest::SetUp(); + // Set random values for the Neo-Hookean parameters between 1000 and 10000 params.C10 = getRandomDouble(1000.0, 10000.0); @@ -74,11 +148,8 @@ 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: @@ -90,54 +161,8 @@ class STRUCT_NeoHookeanTest : public NeoHookeanTest { } }; -// Test PK2 stress zero for F = I -TEST_F(STRUCT_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(STRUCT_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(STRUCT_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(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: @@ -149,50 +174,6 @@ class USTRUCT_NeoHookeanTest : public NeoHookeanTest { } }; -// 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 - create_random_F(F); - TestNH->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); -} - - // ---------------------------------------------------------------------------- // --------------------------- Mooney-Rivlin Material ------------------------- // ---------------------------------------------------------------------------- @@ -202,17 +183,10 @@ TEST_F(USTRUCT_NeoHookeanTest, TestMaterialElasticityConsistentRandomF) { * * This class sets up the necessary parameters and objects for testing the Mooney-Rivlin material model. */ -class MooneyRivlinTest : public ::testing::Test { +class MooneyRivlinTest : public MaterialModelTest { protected: - // Variables common across tests + // Material parameters object MooneyRivlinParams 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-11; // 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 TestMooneyRivlin* TestMR; @@ -220,6 +194,8 @@ class MooneyRivlinTest : public ::testing::Test { // Setup method to initialize variables before each test void SetUp() override { + MaterialModelTest::SetUp(); + // Set random values for the Mooney-Rivlin parameters between 1000 and 10000 params.C01 = getRandomDouble(1000.0, 10000.0); params.C10 = getRandomDouble(1000.0, 10000.0); @@ -236,11 +212,8 @@ 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: @@ -252,54 +225,8 @@ class STRUCT_MooneyRivlinTest : public MooneyRivlinTest { } }; -// Test PK2 stress zero for F = I -TEST_F(STRUCT_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(STRUCT_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(STRUCT_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(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: @@ -311,54 +238,6 @@ class USTRUCT_MooneyRivlinTest : public MooneyRivlinTest { } }; -// 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 - create_random_F(F); - TestMR->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); -} - - -// ---------------------------------------------------------------------------- -// ------------------- Holzapfel-Gasser-Ogden Material ------------------------- -// ---------------------------------------------------------------------------- - // ---------------------------------------------------------------------------- // ----------------------- Holzapfel-Ogden Material --------------------------- @@ -369,16 +248,10 @@ TEST_F(USTRUCT_MooneyRivlinTest, TestMaterialElasticityConsistentRandomF) { * * This class sets up the necessary parameters and objects for testing the Holzapfel-Ogden material model. */ -class HolzapfelOgdenTest : public ::testing::Test { +class HolzapfelOgdenTest : public :: MaterialModelTest { protected: - // Variables common across tests + // Material parameters object HolzapfelOgdenParams 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-11; // 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 TestHolzapfelOgden* TestHO; @@ -386,16 +259,18 @@ class HolzapfelOgdenTest : public ::testing::Test { // Setup method to initialize variables before each test void SetUp() override { + MaterialModelTest::SetUp(); + // Set Holzapfel-Ogden parameters from cardiac benchmark paper params.a = 59.0; // Pa params.a_f = 18472.0; // Pa params.a_s = 2481.0; // Pa params.a_fs = 216.0; // Pa - params.b = 8.023; // Pa - params.b_f = 16.026; // Pa - params.b_s = 11.12; // Pa - params.b_fs = 11.436; // Pa - params.k = 100000.0; // Pa + params.b = 8.023; // no units + params.b_f = 16.026; // no units + params.b_s = 11.12; // no units + params.b_fs = 11.436; // no units + params.k = 100.0; // no units // Set random values for f between 0 and 1 and normalize params.f[0] = getRandomDouble(0.0, 1.0); @@ -428,9 +303,6 @@ class HolzapfelOgdenTest : public ::testing::Test { throw runtime_error("f and s are not orthogonal"); } - // Flag to use full anisotropic invariants for strain energy computation - params.full_anisotropic_invariants = false; - // Initialize the test object TestHO = new TestHolzapfelOgden(params); @@ -444,11 +316,8 @@ 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: @@ -460,240 +329,1468 @@ class STRUCT_HolzapfelOgdenTest : public HolzapfelOgdenTest { } }; -// Test PK2 stress zero for F = I -TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressIdentityF) { - //verbose = true; // Show values of S and S_ref +/** + * @brief Test fixture class for USTRUCT Holzapfel-Ogden material model. + */ +class USTRUCT_HolzapfelOgdenTest : public HolzapfelOgdenTest { +protected: + void SetUp() override { + HolzapfelOgdenTest::SetUp(); - // 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); -} + // Use ustruct + TestHO->ustruct = true; + } +}; -// Test PK2 stress for triaxial stretch -TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressTriaxialStretch) { - //verbose = true; // Show values of S and S_ref +// ---------------------------------------------------------------------------- +// ------------- Holzapfel-Ogden (Modified Anisotropy) Material -------------- +// ---------------------------------------------------------------------------- - // 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); -} +/** + * @brief Test fixture class for the Holzapfel-Ogden (Modified Anisotropy) material model. + * + * This class sets up the necessary parameters and objects for testing the Holzapfel-Ogden (Modified Anisotropy) material model. +*/ +class HolzapfelOgdenMATest : public MaterialModelTest { +protected: + // Material parameters object + HolzapfelOgdenMAParams params; -// Test PK2 stress for triaxial compression -TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressTriaxialCompression) { - //verbose = true; // Show values of S and S_ref + // Add the test object + TestHolzapfelOgdenMA* TestHO_ma; - // 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); -} + // Setup method to initialize variables before each test + void SetUp() override { -// Test PK2 stress with finite difference for random F -TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressRandomF) { - //verbose = true; // Show values of S, dE, SdE and dPsi + MaterialModelTest::SetUp(); - // 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); + // Set Holzapfel-Ogden parameters from cardiac benchmark paper + params.a = 59.0; // Pa + params.a_f = 18472.0; // Pa + params.a_s = 2481.0; // Pa + params.a_fs = 216.0; // Pa + params.b = 8.023; // no units + params.b_f = 16.026; // no units + params.b_s = 11.12; // no units + params.b_fs = 11.436; // no units + params.k = 100.0; // no units - // Check PK2 stress against reference value - TestHO->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); -} + // 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 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]; + params.s[2] = -params.f[1]; + } else { // If f[0] is the dominant component, use another approach + params.s[0] = -params.f[2]; + params.s[1] = 0; + params.s[2] = params.f[0]; + } -// Test PK2 stress consistent with strain energy for random F -TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressConsistentRandomF) { - //verbose = true; // Show values of S, dE, SdE and dPsi + // Normalize s + double norm_s = sqrt(params.s[0]*params.s[0] + params.s[1]*params.s[1] + params.s[2]*params.s[2]); + params.s[0] /= norm_s; params.s[1] /= norm_s; params.s[2] /= norm_s; - // 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); -} + // Check f.s = 0 + double dot_fs = params.f[0]*params.s[0] + params.f[1]*params.s[1] + params.f[2]*params.s[2]; + if (fabs(dot_fs) > 1e-6) { + cout << "f.s = " << dot_fs << endl; + cout << "f = [" << params.f[0] << ", " << params.f[1] << ", " << params.f[2] << "]" << endl; + cout << "s = [" << params.s[0] << ", " << params.s[1] << ", " << params.s[2] << "]" << endl; + throw runtime_error("f and s are not orthogonal"); + } -// Test material elasticity consistent with PK2 stress -TEST_F(STRUCT_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); -} + // Initialize the test object + TestHO_ma = new TestHolzapfelOgdenMA(params); + } + + // TearDown method to clean up after each test, if needed + void TearDown() override { + // Clean up the test object + delete TestHO_ma; + TestHO_ma = nullptr; + } +}; + +/** + * @brief Test fixture class for STRUCT Holzapfel-Ogden material model. + */ +class STRUCT_HolzapfelOgdenMATest : public HolzapfelOgdenMATest { +protected: + void SetUp() override { + HolzapfelOgdenMATest::SetUp(); + + // Use struct + TestHO_ma->ustruct = false; + } +}; -// ------------------------------ 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 { +class USTRUCT_HolzapfelOgdenMATest : public HolzapfelOgdenMATest { protected: void SetUp() override { - HolzapfelOgdenTest::SetUp(); + HolzapfelOgdenMATest::SetUp(); // Use ustruct - TestHO->ustruct = true; + TestHO_ma->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); + + + +// ---------------------------------------------------------------------------- +// ---------------- 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 MaterialModelTest { +protected: + // Material parameters object + VolumetricPenaltyParams params; + + // Add the test object + TestQuadraticVolumetricPenalty* TestQVP; + + // Setup method to initialize variables before each test + void SetUp() override { + + MaterialModelTest::SetUp(); + + // 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; + } +}; + +/** + * @brief Test fixture class for STRUCT Quadratic penalty model. + */ +class STRUCT_QuadraticVolumetricPenaltyTest : public QuadraticVolumetricPenaltyTest { +protected: + void SetUp() override { + QuadraticVolumetricPenaltyTest::SetUp(); + + // Use struct + //TestQVP->ustruct = false; + } +}; + +/** + * @brief Test fixture class for USTRUCT Quadratic penalty model. + */ +class USTRUCT_QuadraticVolumetricPenaltyTest : public QuadraticVolumetricPenaltyTest { +protected: + void SetUp() override { + QuadraticVolumetricPenaltyTest::SetUp(); + + // Use ustruct + //TestQVP->ustruct = true; + } +}; + + +// ---------------------------------------------------------------------------- +// --------------------------- 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 MaterialModelTest { +protected: + // Material parameters object + VolumetricPenaltyParams params; + + // Add the test object + TestSimoTaylor91VolumetricPenalty* TestST91; + + // Setup method to initialize variables before each test + void SetUp() override { + + MaterialModelTest::SetUp(); + + // 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; + } +}; + +/** + * @brief Test fixture class for STRUCT Simo-Taylor91 penalty model. + */ +class STRUCT_SimoTaylor91VolumetricPenaltyTest : public SimoTaylor91VolumetricPenaltyTest { +protected: + void SetUp() override { + SimoTaylor91VolumetricPenaltyTest::SetUp(); + + // Use struct + //TestST91->ustruct = false; + } +}; + +/** + * @brief Test fixture class for USTRUCT Simo-Taylor91 penalty model. + */ +class USTRUCT_SimoTaylor91VolumetricPenaltyTest : public SimoTaylor91VolumetricPenaltyTest { +protected: + void SetUp() override { + SimoTaylor91VolumetricPenaltyTest::SetUp(); + + // Use ustruct + //TestST91->ustruct = true; + } +}; + +// ---------------------------------------------------------------------------- +// ---------------------- 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 MaterialModelTest { +protected: + // Material parameters object + VolumetricPenaltyParams params; + + // Add the test object + TestMiehe94VolumetricPenalty* TestM94; + + // Setup method to initialize variables before each test + void SetUp() override { + + MaterialModelTest::SetUp(); + + // 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; + } +}; + +/** + * @brief Test fixture class for STRUCT Miehe94 penalty model. + */ +class STRUCT_Miehe94VolumetricPenaltyTest : public Miehe94VolumetricPenaltyTest { +protected: + void SetUp() override { + Miehe94VolumetricPenaltyTest::SetUp(); + + // Use struct + //TestM94->ustruct = false; + } +}; + +/** + * @brief Test fixture class for USTRUCT Miehe94 penalty model. + */ +class USTRUCT_Miehe94VolumetricPenaltyTest : public Miehe94VolumetricPenaltyTest { +protected: + void SetUp() override { + Miehe94VolumetricPenaltyTest::SetUp(); + + // Use ustruct + //TestM94->ustruct = true; + } +}; + + + + + + + + + +// ============================================================================ +// ------------------------------- TESTS -------------------------------------- +// ============================================================================ + + + +// ---------------------------------------------------------------------------- +// --------------------------- Neo-Hookean Material --------------------------- +// ---------------------------------------------------------------------------- + +// ------------------------------ STRUCT Tests -------------------------------- + +// Test PK2 stress zero for F = I +TEST_F(STRUCT_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 order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(STRUCT_NeoHookeanTest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestNH->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(STRUCT_NeoHookeanTest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestNH->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(STRUCT_NeoHookeanTest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestNH->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(STRUCT_NeoHookeanTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestNH->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(STRUCT_NeoHookeanTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestNH->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(STRUCT_NeoHookeanTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestNH->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// ------------------------------ USTRUCT Tests -------------------------------- + +// 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 order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(USTRUCT_NeoHookeanTest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestNH->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(USTRUCT_NeoHookeanTest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestNH->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(USTRUCT_NeoHookeanTest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestNH->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(USTRUCT_NeoHookeanTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestNH->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(USTRUCT_NeoHookeanTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestNH->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(USTRUCT_NeoHookeanTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestNH->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + + + +// ---------------------------------------------------------------------------- +// --------------------------- Mooney-Rivlin Material ------------------------- +// ---------------------------------------------------------------------------- + +// ------------------------------ STRUCT Tests -------------------------------- + +// Test PK2 stress zero for F = I +TEST_F(STRUCT_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 order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(STRUCT_MooneyRivlinTest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestMR->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(STRUCT_MooneyRivlinTest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestMR->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(STRUCT_MooneyRivlinTest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestMR->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(STRUCT_MooneyRivlinTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestMR->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(STRUCT_MooneyRivlinTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestMR->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(STRUCT_MooneyRivlinTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestMR->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + + +// ------------------------------ USTRUCT Tests -------------------------------- + +// 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 order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(USTRUCT_MooneyRivlinTest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestMR->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(USTRUCT_MooneyRivlinTest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestMR->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(USTRUCT_MooneyRivlinTest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestMR->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(USTRUCT_MooneyRivlinTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestMR->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(USTRUCT_MooneyRivlinTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestMR->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(USTRUCT_MooneyRivlinTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestMR->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + + + +// ---------------------------------------------------------------------------- +// ----------------------- Holzapfel-Ogden Material --------------------------- +// ---------------------------------------------------------------------------- + +// ------------------------------ STRUCT Tests -------------------------------- + +// Test PK2 stress zero for F = I +TEST_F(STRUCT_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 order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial stretch +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderTriaxialStretch) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial stretch + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.3}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial compression +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderTriaxialCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial compression + double F[3][3] = {{0.9, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 0.7}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for biaxial stretch/compression +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderBiaxialStretchCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for biaxial stretch/compression + double F[3][3] = {{1.2, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 1.0}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(STRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} +//triaxial extension, compression and biaxial extension - 3 new tests; struct and ustruct HO and HO-ma. +// Test order of convergence of consistency of material elasticity for triaxial stretch +TEST_F(STRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderTriaxialStretch) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial stretch + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.3}}; + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial compression +TEST_F(STRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderTriaxialCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial compression + double F[3][3] = {{0.9, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 0.7}}; + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for biaxial stretch/compression +TEST_F(STRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderBiaxialStretchCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for biaxial stretch/compression + double F[3][3] = {{1.2, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 1.0}}; + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(STRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(STRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(STRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// ------------------------------ USTRUCT Tests -------------------------------- + +// 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 order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial stretch +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderTriaxialStretch) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial stretch + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.3}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial compression +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderTriaxialCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial compression + double F[3][3] = {{0.9, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 0.7}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for biaxial stretch/compression +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderBiaxialStretchCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for biaxial stretch/compression + double F[3][3] = {{1.2, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 1.0}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for triaxial stretch +TEST_F(USTRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderTriaxialStretch) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial stretch + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.3}}; + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial compression +TEST_F(USTRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderTriaxialCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial compression + double F[3][3] = {{0.9, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 0.7}}; + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for biaxial stretch/compression +TEST_F(USTRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderBiaxialStretchCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for biaxial stretch/compression + double F[3][3] = {{1.2, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 1.0}}; + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(USTRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(USTRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(USTRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestHO->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + + + +// ---------------------------------------------------------------------------- +// --------------- Holzapfel-Ogden (Modified Anisotropy) Material ------------- +// ---------------------------------------------------------------------------- + +// ------------------------------ STRUCT Tests -------------------------------- + +// Test PK2 stress zero for F = I +TEST_F(STRUCT_HolzapfelOgdenMATest, 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_ma->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial stretch +TEST_F(STRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderTriaxialStretch) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial stretch + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.3}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial compression +TEST_F(STRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderTriaxialCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial compression + double F[3][3] = {{0.9, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 0.7}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for biaxial stretch/compression +TEST_F(STRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderBiaxialStretchCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for biaxial stretch/compression + double F[3][3] = {{1.2, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 1.0}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(STRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(STRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(STRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for triaxial stretch +TEST_F(STRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderTriaxialStretch) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial stretch + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.3}}; + + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial compression +TEST_F(STRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderTriaxialCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial compression + double F[3][3] = {{0.9, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 0.7}}; + + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for biaxial stretch/compression +TEST_F(STRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderBiaxialStretchCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for biaxial stretch/compression + double F[3][3] = {{1.2, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 1.0}}; + + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(STRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(STRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Set convergence order tolerances larger and special delta_max/delta_min to get this test to pass + convergence_order_tol = 0.15; + delta_max = 8e-6; + delta_min = 2e-6; + + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(STRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Set convergence order tolerances larger and special delta_max/delta_min to get this test to pass + //convergence_order_tol = 0.02; + delta_max = 4e-5; + delta_min = 1e-5; + + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + + +// ------------------------------ USTRUCT Tests -------------------------------- + +// Test PK2 stress zero for F = I +TEST_F(USTRUCT_HolzapfelOgdenMATest, 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_ma->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial stretch +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderTriaxialStretch) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial stretch + double F[3][3] = {{1.1, 0.0, 0.0}, + {0.0, 1.2, 0.0}, + {0.0, 0.0, 1.3}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial compression +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderTriaxialCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for triaxial compression + double F[3][3] = {{0.9, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 0.7}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for biaxial stretch/compression +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderBiaxialStretchCompression) { + //verbose = true; // Show order of convergence, errors, F, S + + // Create a deformation gradient F for biaxial stretch/compression + double F[3][3] = {{1.2, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 1.0}}; + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); +} + + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} + +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestHO_ma->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } } -// Test PK2 stress for triaxial stretch -TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressTriaxialStretch) { - //verbose = true; // Show values of S and S_ref +// Test order of convergence of consistency of material elasticity for triaxial stretch +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderTriaxialStretch) { + //verbose = true; // Show order of convergence, errors, F, S - // Check triaxial stretch produces PK2 stress consistent with svFSI + // Create a deformation gradient F for triaxial stretch 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); + + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); } -// Test PK2 stress for triaxial compression -TEST_F(USTRUCT_HolzapfelOgdenTest, TestPK2StressTriaxialCompression) { - //verbose = true; // Show values of S and S_ref +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for triaxial compression +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderTriaxialCompression) { + //verbose = true; // Show order of convergence, errors, F, S - // Check triaxial compression produces PK2 stress consistent with svFSI + // Create a deformation gradient F for triaxial compression 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); + {0.0, 0.8, 0.0}, + {0.0, 0.0, 0.7}}; + + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_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 +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for biaxial stretch/compression +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderBiaxialStretchCompression) { + //verbose = true; // Show order of convergence, errors, F, S - // 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); + // Create a deformation gradient F for biaxial stretch/compression + double F[3][3] = {{1.2, 0.0, 0.0}, + {0.0, 0.8, 0.0}, + {0.0, 0.0, 1.0}}; - // Check PK2 stress against reference value - TestHO->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_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 order of convergence of consistency of material elasticity for random F (small) +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S -// Test material elasticity consistent with PK2 stress -TEST_F(USTRUCT_HolzapfelOgdenTest, TestMaterialElasticityConsistentRandomF) { - //verbose = true; // Show values of CC, dE, CCdE and dS + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); - // 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); + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } } +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} -// ---------------------------------------------------------------------------- -// ---------------------------------------------------------------------------- -// ----------------------- Volumetric penalty models ------------------------- -// ---------------------------------------------------------------------------- -// ---------------------------------------------------------------------------- +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(USTRUCT_HolzapfelOgdenMATest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S -// ---------------------------------------------------------------------------- -// --------------------------- 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 + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); - // Add the test object - TestQuadraticVolumetricPenalty* TestQVP; + // Check order of convergence of consistency of material elasticity + TestHO_ma->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } +} - // 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(); +// ---------------------------------------------------------------------------- +// ---------------------- Quadratic Volumetric Penalty Material ---------------- +// ---------------------------------------------------------------------------- - // Use struct - //TestQVP->ustruct = false; - } -}; +// ------------------------------ STRUCT Tests -------------------------------- // Test PK2 stress zero for F = I TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressIdentityF) { @@ -719,73 +1816,93 @@ TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDe 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); +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S - // 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; - } + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestQVP->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); } +} - // 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 order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestQVP->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_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 +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S - // 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); + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); - // Check PK2 stress against reference value - TestQVP->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestQVP->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_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 +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); - // Check random F produces consistent PK2 stress - create_random_F(F); - TestQVP->testPK2StressConsistentWithStrainEnergy(F, n_iter, rel_tol, abs_tol, delta, verbose); + // Check order of convergence of consistency of material elasticity + TestQVP->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } } -// Test material elasticity consistent with PK2 stress -TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) { - //verbose = true; // Show values of CC, dE, CCdE and dS +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); - // Check with random F - create_random_F(F); - TestQVP->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); + // Check order of convergence of consistency of material elasticity + TestQVP->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, 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(); +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S - // Use ustruct - //TestQVP->ustruct = true; + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestQVP->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); } -}; +} + + +// ------------------------------ USTRUCT Tests -------------------------------- // Test rho, beta, drho/dp and dbeta/dp for random pressure TEST_F(USTRUCT_QuadraticVolumetricPenaltyTest, TestRhoBeta) { @@ -809,62 +1926,13 @@ TEST_F(USTRUCT_QuadraticVolumetricPenaltyTest, TestRhoBeta) { TestQVP->testRhoBetaAgainstReference(p, rho0, rho_ref, beta_ref, drhodp_ref, dbetadp_ref, rel_tol, abs_tol, verbose); } + // ---------------------------------------------------------------------------- -// --------------------------- Simo-Taylor91 Volumetric Penalty Model --------- +// ---------------------- Simo-Taylor 91 Volumetric Penalty Material ------------ // ---------------------------------------------------------------------------- -/** - * @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(); +// ------------------------------ STRUCT Tests -------------------------------- - // Use struct - //TestST91->ustruct = false; - } -}; // Test PK2 stress zero for F = I TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressIdentityF) { @@ -890,73 +1958,92 @@ TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressPrescribedIsochori 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); +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S - // 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; - } + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestST91->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); } +} - // 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 order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestST91->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_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 +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S - // 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); + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); - // Check PK2 stress against reference value - TestST91->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestST91->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_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 +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); - // Check random F produces consistent PK2 stress - create_random_F(F); - TestST91->testPK2StressConsistentWithStrainEnergy(F, n_iter, rel_tol, abs_tol, delta, verbose); + // Check order of convergence of consistency of material elasticity + TestST91->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } } -// Test material elasticity consistent with PK2 stress -TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) { - //verbose = true; // Show values of CC, dE, CCdE and dS +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); - // Check with random F - create_random_F(F); - TestST91->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); + // Check order of convergence of consistency of material elasticity + TestST91->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, 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(); +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S - // Use ustruct - //TestST91->ustruct = true; + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestST91->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); } -}; +} + +// ------------------------------ USTRUCT Tests -------------------------------- // Test rho and beta values for random p TEST_F(USTRUCT_SimoTaylor91VolumetricPenaltyTest, TestRhoBeta) { @@ -980,61 +2067,12 @@ TEST_F(USTRUCT_SimoTaylor91VolumetricPenaltyTest, TestRhoBeta) { TestST91->testRhoBetaAgainstReference(p, rho0, rho_ref, beta_ref, drhodp_ref, dbetadp_ref, rel_tol, abs_tol, verbose); } + // ---------------------------------------------------------------------------- -// --------------------------- Miehe94 Volumetric Penalty Model --------------- +// ---------------------- Miehe 94 Volumetric Penalty Material ----------------- // ---------------------------------------------------------------------------- -/** - * @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; - } -}; +// ------------------------------ STRUCT Tests -------------------------------- // Test PK2 stress zero for F = I TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressIdentityF) { @@ -1060,73 +2098,91 @@ TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDefo 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 +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (small) +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S - // Create random deformation gradient tensor - create_random_F(F); + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, 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 order of convergence between finite difference and get_pk2cc() PK2 stress + TestM94->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); } +} - // 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 order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (medium) +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S + + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestM94->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_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 +// Test order of convergence between finite difference PK2 stress and get_pk2cc() PK2 stress for random F (large) +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S - // 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); + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); - // Check PK2 stress against reference value - TestM94->testPK2StressAgainstReference(F, S_ref, rel_tol, abs_tol, verbose); + // Check order of convergence between finite difference and get_pk2cc() PK2 stress + TestM94->testPK2StressConvergenceOrder(F, delta_max, delta_min, order, convergence_order_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 +// Test order of convergence of consistency of material elasticity for random F (small) +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFSmall) { + //verbose = true; // Show order of convergence, errors, F, S - // Check random F produces consistent PK2 stress - create_random_F(F); - TestM94->testPK2StressConsistentWithStrainEnergy(F, n_iter, rel_tol, abs_tol, delta, verbose); + // Loop over F in F_small_list + for (auto F_std : F_small_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestM94->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); + } } -// Test material elasticity consistent with PK2 stress -TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) { - //verbose = true; // Show values of CC, dE, CCdE and dS +// Test order of convergence of consistency of material elasticity for random F (medium) +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFMedium) { + //verbose = true; // Show order of convergence, errors, F, S - // Check with random F - create_random_F(F); - TestM94->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose); + // Loop over F in F_medium_list + for (auto F_std : F_medium_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestM94->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, 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(); +// Test order of convergence of consistency of material elasticity for random F (large) +TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestMaterialElasticityConsistencyConvergenceOrderRandomFLarge) { + //verbose = true; // Show order of convergence, errors, F, S - // Use ustruct - //TestM94->ustruct = true; + // Loop over F in F_large_list + for (auto F_std : F_large_list) { + // Convert to C array + convertToCArray(F_std, F); + + // Check order of convergence of consistency of material elasticity + TestM94->testMaterialElasticityConsistencyConvergenceOrder(F, delta_max, delta_min, order, convergence_order_tol, verbose); } -}; +} + +// ------------------------------ USTRUCT Tests -------------------------------- // Test rho and beta values for random p TEST_F(USTRUCT_Miehe94VolumetricPenaltyTest, TestRhoBeta) { @@ -1150,7 +2206,8 @@ TEST_F(USTRUCT_Miehe94VolumetricPenaltyTest, TestRhoBeta) { TestM94->testRhoBetaAgainstReference(p, rho0, rho_ref, beta_ref, drhodp_ref, dbetadp_ref, rel_tol, abs_tol, verbose); } -// ---------------------------------------------------------------------------- -// ---------------------------------------------------------------------------- +// ============================================================================ +// ========================== END OF TESTS ==================================== +// ============================================================================ diff --git a/tests/unitTests/test.h b/tests/unitTests/test.h index 0f95e7b2..97a6f26c 100644 --- a/tests/unitTests/test.h +++ b/tests/unitTests/test.h @@ -61,86 +61,6 @@ #include "mat_models.h" #include "mat_models_carray.h" -// Class to contain material parameters -class MatParams { -public: - virtual ~MatParams() {} // Virtual destructor for proper cleanup -}; - -// Class to contain Neo-Hookean material parameters -class NeoHookeanParams : public MatParams { -public: - double C10; - - // Default constructor - NeoHookeanParams() : C10(0.0) {} - - // Constructor with parameters - NeoHookeanParams(double c10) : C10(c10) {} - -}; - -// Class to contain Mooney-Rivlin material parameters -class MooneyRivlinParams : public MatParams { -public: - double C01; - double C10; - - // Default constructor - MooneyRivlinParams() : C01(0.0), C10(0.0) {} - - // Constructor with parameters - MooneyRivlinParams(double c01, double c10) : C01(c01), C10(c10) {} - -}; - -// Class to contain Holzapfel-Ogden material parameters -class HolzapfelOgdenParams : public MatParams { -public: - double a; - double b; - double a_f; - double b_f; - double a_s; - double b_s; - double a_fs; - double b_fs; - double f[3]; // Fiber direction - double s[3]; // Sheet direction - - double k; // Smoothed Heaviside function parameter - - 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), k(0.0) { - for (int i = 0; i < 3; i++) { - f[i] = 0.0; - s[i] = 0.0; - } - } - - // 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 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]; - } - } -}; - -// 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 ---------------------- // -------------------------------------------------------------- @@ -160,6 +80,19 @@ void create_identity_F(double F[N][N]) { } } +/** + * @brief Create a ones matrix. + * + */ +template +void create_ones_matrix(double A[N][N]) { + for (int i = 0; i < N; i++) { + for (int J = 0; J < N; J++) { + A[i][J] = 1.0; + } + } +} + /** * @brief Generates a random double value. * @@ -372,6 +305,25 @@ solidMechanicsTerms calcSolidMechanicsTerms(const double F[N][N]) { return out; } +/** + * @brief Computes a linear regression line y = mx + b for given x and y data. + * + * @param x x data points. + * @param y y data points. + * @return std::pair A pair containing the slope (m) and the y-intercept (b). + */ +std::pair computeLinearRegression(const std::vector& x, const std::vector& y) { + int n = x.size(); + double sum_x = std::accumulate(x.begin(), x.end(), 0.0); + double sum_y = std::accumulate(y.begin(), y.end(), 0.0); + double sum_xx = std::inner_product(x.begin(), x.end(), x.begin(), 0.0); + double sum_xy = std::inner_product(x.begin(), x.end(), y.begin(), 0.0); + + double m = (n * sum_xy - sum_x * sum_y) / (n * sum_xx - sum_x * sum_x); + double b = (sum_y - m * sum_x) / n; + + return std::make_pair(m, b); +} // -------------------------------------------------------------- // -------------------- Mock svFSIplus object ------------------- @@ -538,7 +490,7 @@ class TestMaterialModel { * 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: + * Pseudocode (for first order finite difference): * - Compute strain energy density Psi(F) * - For each component of F, F[i][J] * - Perturb F[i][J] by delta to get F_tilde @@ -550,37 +502,68 @@ class TestMaterialModel { * @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[in] order The order of the finite difference scheme (1 for first order, 2 for second order, etc.). * @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]) { + void calcPK2StressFiniteDifference(const double F[N][N], const double delta, const int order, double (&S)[N][N]) { // Compute strain energy density given F double Psi = computeStrainEnergy(F); // Compute 1st PK stress P_iJ = dPsi / dF[i][J] using finite difference, component by component double P[3][3] = {}; - double F_tilde[N][N]; // perturbed deformation gradient - for (int i = 0; i < N; i++) { - for (int J = 0; J < N; J++) { - // Perturb the iJ-th component of F by delta - for (int k = 0; k < N; k++) { - for (int l = 0; l < N; l++) { - F_tilde[k][l] = F[k][l]; + if (order == 1){ + double F_tilde[N][N]; // perturbed deformation gradient + for (int i = 0; i < N; i++) { + for (int J = 0; J < N; J++) { + // Perturb the iJ-th component of F by delta + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + F_tilde[k][l] = F[k][l]; + } } + F_tilde[i][J] += delta; + + // Compute Psi_MR for perturbed deformation gradient + double Psi_tilde = computeStrainEnergy(F_tilde); + + // Compute differences in Psi + double dPsi = Psi_tilde - Psi; + + // Compute P[i][J] = dPsi / dF[i][J] + P[i][J] = dPsi / delta; } - F_tilde[i][J] += delta; + } + } + else if (order == 2){ + double F_plus[N][N]; // positive perturbed deformation gradient + double F_minus[N][N]; // negative perturbed deformation gradient + for (int i = 0; i < N; i++) { + for (int J = 0; J < N; J++) { + // Perturb the iJ-th component of F by +-delta + for (int k = 0; k < N; k++) { + for (int l = 0; l < N; l++) { + F_plus[k][l] = F[k][l]; + F_minus[k][l] = F[k][l]; + } + } + F_plus[i][J] += delta; + F_minus[i][J] -= delta; - // Compute Psi_MR for perturbed deformation gradient - double Psi_tilde = computeStrainEnergy(F_tilde); + // Compute Psi_MR for perturbed deformation gradient + double Psi_plus = computeStrainEnergy(F_plus); + double Psi_minus = computeStrainEnergy(F_minus); - // Compute differences in Psi - double dPsi = Psi_tilde - Psi; + // Compute differences in Psi + double dPsi = Psi_plus - Psi_minus; - // Compute P[i][J] = dPsi / dF[i][J] - P[i][J] = dPsi / delta; + // Compute P[i][J] = dPsi / dF[i][J] + P[i][J] = dPsi / (2.0 * delta); + } } } + // Compute S_ref = F^-1 * P_ref double F_inv[N][N]; @@ -589,202 +572,302 @@ class TestMaterialModel { } /** - * @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 + * @brief Computes the PK2 stress tensor S(F) from the strain energy density Psi(F) using finite differences and checks the order of convergence. * * @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. - * + * @param[in] delta_min Minimum perturbation scaling factor. + * @param[in] delta_max Maximum perturbation scaling factor. + * @param[in] order Order of the finite difference scheme (1 for first order, 2 for second order, etc.). + * @param[in] convergence_order_tol Tolerance for comparing convergence order with expected value + * @param[in] verbose Show values error and order of convergence if true. */ - 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]; - calc_JCE(F, J, C, E); + template + void testPK2StressConvergenceOrder(const double F[N][N], const double delta_max, const double delta_min, const int order, const double convergence_order_tol, const bool verbose = false) { + // Check delta_max > delta_min + if (delta_max <= delta_min) { + std::cerr << "Error: delta_max must be greater than delta_min." << std::endl; + return; + } - // Compute Psi(F) - double Psi = computeStrainEnergy(F); + // Check order is 1 or 2 + if (order != 1 && order != 2) { + std::cerr << "Error: order must be 1 or 2." << std::endl; + return; + } + + // Create list of deltas for convergence test (delta = delta_max, delta_max/2, delta_max/4, ...) + std::vector deltas; + double delta = delta_max; + while (delta >= delta_min) { + deltas.push_back(delta); + delta /= 2.0; + } - // Compute S(F) from svFSIplus + // Compute S(F) from get_pk2cc() double S[3][3], Dm[6][6]; get_pk2cc(F, S, Dm); - // Generate many random dF and check that S:dE = dPsi - // S was obtained from get_pk2cc(), and dPsi = Psi(F + dF) - Psi(F) - double dPsi, dE[3][3], SdE; - double F_tilde[3][3], J_tilde, C_tilde[3][3], E_tilde[3][3]; - for (int i = 0; i < n_iter; i++) { - // Perturb the deformation gradient - perturb_random_F(F, delta, F_tilde); + // Compute finite difference S for each delta and store error in list + std::vector errors; + double S_fd[3][3]; + for (int i = 0; i < deltas.size(); i++) { + calcPK2StressFiniteDifference(F, deltas[i], order, S_fd); - // Compute perturbed E_tilde from F_tilde - calc_JCE(F_tilde, J_tilde, C_tilde, E_tilde); + // Compute Frobenius norm of error between S and S_fd + double error = 0.0; + for (int I = 0; I < 3; I++) { + for (int J = 0; J < 3; J++) { + error += pow(S[I][J] - S_fd[I][J], 2); + } + } + error = sqrt(error); - // Compute Psi(F_tilde) - double Psi_tilde = computeStrainEnergy(F_tilde); + // Store error in list + errors.push_back(error); + } - // Compute dPsi = Psi(F_tilde) - Psi(F) - dPsi = Psi_tilde - Psi; + // Compute order of convergence by fitting a line to log(delta) vs log(error) + std::vector log_deltas, log_errors; + for (int i = 0; i < deltas.size(); i++) { + log_deltas.push_back(log(deltas[i])); + log_errors.push_back(log(errors[i])); + } + + // Fit a line to log(delta) vs log(error) + // m is the slope (order of convergence), b is the intercept + auto [m, b] = computeLinearRegression(log_deltas, log_errors); + + // Check that order of convergence is > order - convergence_order_tol + EXPECT_GT(m, order - convergence_order_tol); - // Compute dE = E(F_tilde) - E(F) + // Print results if verbose + if (verbose) { + std::cout << "Slope (order of convergence): " << m << std::endl; + std::cout << "Intercept: " << b << std::endl; + std::cout << "Errors: "; + for (int i = 0; i < errors.size(); i++) { + std::cout << errors[i] << " "; + } + std::cout << std::endl; + std::cout << std::endl; + + std::cout << "F = " << std::endl; for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - dE[i][j] = E_tilde[i][j] - E[i][j]; + for (int J = 0; J < 3; J++) { + std::cout << F[i][J] << " "; } + std::cout << std::endl; } - // Compute S:dE - double SdE = mat_fun_carray::mat_ddot<3>(S, dE); + std::cout << "S = " << std::endl; + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + std::cout << S[i][J] << " "; + } + std::cout << std::endl; + } - // Check that S:dE = dPsi - EXPECT_NEAR(SdE, dPsi, fmax(abs_tol, rel_tol * fabs(dPsi))); - - // Print results if verbose - if (verbose) { - std::cout << "Iteration " << i << ":" << std::endl; + std::cout << std::endl; + } + } - printMaterialParameters(); + /** + * @brief Compute perturbation in strain energy density (dPsi) given perturbation in the deformation gradient (dF). + * + * @param F Deformation gradient + * @param dF Deformation gradient perturbation shape + * @param delta Deformation gradient perturbation scaling factor + * @param order Order of the finite difference scheme (1 for first order, 2 for second order, etc.) + * @param dPsi Strain energy density perturbation + */ + void calcdPsiFiniteDifference(const double F[3][3], const double dF[3][3], const double delta, const int order, double &dPsi) { - std::cout << "F =" << std::endl; - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - std::cout << F[i][j] << " "; - } - std::cout << std::endl; - } + // Compute strain energy density given F + double Psi = computeStrainEnergy(F); - std::cout << "S =" << std::endl; - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - std::cout << S[i][j] << " "; + // Compute dPsi using finite difference, given dF + if (order == 1){ + double F_tilde[3][3]; // perturbed deformation gradient + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + // Perturb the iJ-th component of F by delta * dF[i][J] + for (int k = 0; k < 3; k++) { + for (int l = 0; l < 3; l++) { + F_tilde[k][l] = F[k][l] + delta * dF[k][l]; + } } - std::cout << std::endl; } + } - std::cout << "dE =" << std::endl; - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - std::cout << dE[i][j] << " "; + // Compute Psi_tilde for perturbed deformation gradient + double Psi_tilde = computeStrainEnergy(F_tilde); + + // Compute differences in Psi + dPsi = Psi_tilde - Psi; + } + else if (order == 2){ + double F_plus[3][3]; // positive perturbed deformation gradient + double F_minus[3][3]; // negative perturbed deformation gradient + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + // Perturb the iJ-th component of F by +-delta * dF[i][J] + for (int k = 0; k < 3; k++) { + for (int l = 0; l < 3; l++) { + F_plus[k][l] = F[k][l] + delta * dF[k][l]; + F_minus[k][l] = F[k][l] - delta * dF[k][l]; + } } - std::cout << std::endl; } - - std::cout << "SdE = " << SdE << ", dPsi = " << dPsi << std::endl; - std::cout << std::endl; } + + // Compute Psi_plus and Psi_minus for perturbed deformation gradient + double Psi_plus = computeStrainEnergy(F_plus); + double Psi_minus = computeStrainEnergy(F_minus); + + // Compute differences in Psi + dPsi = Psi_plus - Psi_minus; } } /** - * @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 + * @brief Compute perturbed Green-Lagrange strain tensor (dE) given perturbed deformation gradient (dF) using finite differences * - * @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. + * @param F Deformation gradient + * @param dF Deformation gradient perturbation shape + * @param delta Deformation gradient perturbation scaling factor + * @param order Order of the finite difference scheme (1 for first order, 2 for second order, etc.) + * @param dE Green-Lagrange strain tensor perturbation */ - void testMaterialElasticityConsistentWithPK2Stress(double F[3][3], int n_iter, double rel_tol, double abs_tol, double delta, bool verbose = false) { - + void calcdEFiniteDifference(const double F[3][3], const double dF[3][3], const double delta, const int order, double (&dE)[3][3]) { // Compute E from F double J, C[3][3], E[3][3]; calc_JCE(F, J, C, E); - // Compute S_ij(F) - // Compute CC_ijkl(F). - // CC is provided in Voigt notation as Dm, and we will convert it to CC - double S[3][3], Dm[6][6]; - get_pk2cc(F, S, Dm); // S from svFSI + // Compute dE using finite difference, given dF + if (order == 1){ + double F_tilde[3][3]; // perturbed deformation gradient + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + // Perturb the iJ-th component of F by delta * dF[i][J] + for (int k = 0; k < 3; k++) { + for (int l = 0; l < 3; l++) { + F_tilde[k][l] = F[k][l] + delta * dF[k][l]; + } + } + } + } - // 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++) { - EXPECT_NEAR(Dm_check[i][j], Dm[i][j], abs_tol); - } - } - // ------------------------------- - - // 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() - double dS[3][3], dE[3][3], CCdE[3][3]; - double F_tilde[3][3], J_tilde, C_tilde[3][3], E_tilde[3][3]; - double S_tilde[3][3], Dm_tilde[6][6]; - - // Loop over many random perturbations to the deformation gradient - 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 + double J_tilde, C_tilde[3][3], E_tilde[3][3]; calc_JCE(F_tilde, J_tilde, C_tilde, E_tilde); - - // Compute dE + + // Compute differences in E for (int i = 0; i < 3; i++) { - for (int J = 0; J < 3; J++) { - dE[i][J] = E_tilde[i][J] - E[i][J]; + 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]; + } + else if (order == 2){ + double F_plus[3][3]; // positive perturbed deformation gradient + double F_minus[3][3]; // negative perturbed deformation gradient for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - dS[i][j] = S_tilde[i][j] - S[i][j]; + for (int J = 0; J < 3; J++) { + // Perturb the iJ-th component of F by +-delta * dF[i][J] + for (int k = 0; k < 3; k++) { + for (int l = 0; l < 3; l++) { + F_plus[k][l] = F[k][l] + delta * dF[k][l]; + F_minus[k][l] = F[k][l] - delta * dF[k][l]; + } + } } } - - // Check that CC_ijkl dE_kl = dS_ij - mat_fun_carray::ten_mat_ddot<3>(CC, dE, CCdE); + + // Compute perturbed E_plus and E_minus from F_plus and F_minus + double J_plus, C_plus[3][3], E_plus[3][3]; + double J_minus, C_minus[3][3], E_minus[3][3]; + calc_JCE(F_plus, J_plus, C_plus, E_plus); + calc_JCE(F_minus, J_minus, C_minus, E_minus); + + // Compute differences in E for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { - EXPECT_NEAR(CCdE[i][j], dS[i][j], fmax(abs_tol, rel_tol * fabs(dS[i][j]))); + dE[i][j] = E_plus[i][j] - E_minus[i][j]; } } - + } + } + + /** + * @brief Compute contraction of PK2 stress with perturbation in Green-Lagrange strain tensor (S:dE) given perturbation in deformation gradient (dF) using finite differences. + * + * @param F Deformation gradient + * @param dF Deformation gradient perturbation shape + * @param delta Deformation gradient perturbation scaling factor + * @param order Order of the finite difference scheme (1 for first order, 2 for second order, etc.) + * @param SdE PK2 stress tensor times the perturbation in the Green-Lagrange strain tensor + */ + void calcSdEFiniteDifference(const double F[3][3], const double dF[3][3], const double delta, const int order, double &SdE) { + // Compute S(F) from get_pk2cc() + double S[3][3], Dm[6][6]; + get_pk2cc(F, S, Dm); + + // Compute dE using finite difference, given dF + double dE[3][3]; + calcdEFiniteDifference(F, dF, delta, order, dE); + + // Compute S:dE + SdE = mat_fun_carray::mat_ddot<3>(S, dE); + } + + /** + * @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 = E(F + dF) - E(F) + * - 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) { + int order = 2; + + // Generate many random dF and check that S:dE = dPsi + // S was obtained from get_pk2cc(), and dPsi = Psi(F + dF) - Psi(F) + double dPsi, SdE; + for (int i = 0; i < n_iter; i++) { + // Generate random dF + double dF[3][3]; + create_random_F(dF, 0.0, 1.0); + + // Compute dPsi + calcdPsiFiniteDifference(F, dF, delta, order, dPsi); + + // Compute SdE + calcSdEFiniteDifference(F, dF, delta, order, SdE); + + // Check that S:dE = dPsi + EXPECT_NEAR(SdE, dPsi, fmax(abs_tol, rel_tol * fabs(dPsi))); + // 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++) { @@ -792,44 +875,445 @@ class TestMaterialModel { } std::cout << std::endl; } + + std::cout << "SdE = " << SdE << ", dPsi = " << dPsi << std::endl; + std::cout << std::endl; + } + } + } + + /** + * @brief Tests the order of convergence of the consistency between dPsi and S:dE using finite differences. + * + * Analytically, we should have S = dPsi/dE. This function determines the order of convergence of 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 each component-wise perturbation dF + * - Compute dPsi + * - Compute dE + * - Compute error S:dE - dPsi + * - Compute order of convergence by fitting a line to log(delta) vs log(error) + * + * Note that the order of convergence should be order + 1, because we are comparing differences (dPsi and S:dE) + * instead of derivatives (e.g. dPsi/dF and S:dE/dF). + * @param[in] F Deformation gradient. + * @param[in] delta_max Maximum perturbation scaling factor. + * @param[in] delta_min Minimum perturbation scaling factor. + * @param[in] order Order of the finite difference scheme (1 for first order, 2 for second order, etc.). + * @param[in] convergence_order_tol Tolerance for comparing convergence order with expected value + * @param verbose Show values of errors and order of convergence if true. + */ + void testPK2StressConsistencyConvergenceOrder(double F[3][3], double delta_max, double delta_min, int order, const double convergence_order_tol, bool verbose = false) { + // Check that delta_max > delta_min + if (delta_max <= delta_min) { + std::cerr << "Error: delta_max must be greater than delta_min." << std::endl; + return; + } + + // Check that order is 1 or 2 + if (order != 1 && order != 2) { + std::cerr << "Error: order must be 1 or 2." << std::endl; + return; + } + // Loop over perturbations to each component of F, dF + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + // Generate dF with 1.0 in (i,j) component + double dF[3][3] = {}; + dF[i][j] = 1.0; + + // Create list of deltas for convergence test (delta = delta_max, delta_max/2, delta_max/4, ...) + std::vector deltas; + double delta = delta_max; + while (delta >= delta_min) { + deltas.push_back(delta); + delta /= 2.0; + } + + // Compute dPsi and S:dE for each delta and store error in list + std::vector errors; + double dPsi, SdE; + + for (int i = 0; i < deltas.size(); i++) { + calcdPsiFiniteDifference(F, dF, deltas[i], order, dPsi); + calcSdEFiniteDifference(F, dF, deltas[i], order, SdE); + + // Compute error between dPsi and S:dE + double error = fabs(dPsi - SdE); + + // Store error in list + errors.push_back(error); + } + + // Compute order of convergence by fitting a line to log(delta) vs log(error) + std::vector log_deltas, log_errors; + for (int i = 0; i < deltas.size(); i++) { + log_deltas.push_back(log(deltas[i])); + log_errors.push_back(log(errors[i])); + } + + // Fit a line to log(delta) vs log(error) + // m is the slope (order of convergence), b is the intercept + auto [m, b] = computeLinearRegression(log_deltas, log_errors); + + // Check that order of convergence is > (order + 1) - convergence_order_tol + EXPECT_GT(m, order + 1 - convergence_order_tol); + + // Print results if verbose + if (verbose) { + std::cout << "Slope (order of convergence): " << m << std::endl; + std::cout << "Intercept: " << b << std::endl; + std::cout << "Errors: "; + for (int i = 0; i < errors.size(); i++) { + std::cout << errors[i] << " "; + } + std::cout << std::endl; + std::cout << std::endl; + + std::cout << "F = " << std::endl; + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + std::cout << F[i][J] << " "; + } + std::cout << std::endl; + } + } + } + } + } + + /** + * @brief Compute perturbation in PK2 stress (dS) given perturbation in deformation gradient (dF) using finite differences + * + * @param F Deformation gradient + * @param dF Deformation gradient perturbation shape + * @param delta Deformation gradient perturbation scaling factor + * @param order Order of the finite difference scheme (1 for first order, 2 for second order, etc.) + * @param dS PK2 stress tensor perturbation + */ + void calcdSFiniteDifference(const double F[3][3], const double dF[3][3], const double delta, const int order, double (&dS)[3][3]) { + // Compute S(F) from get_pk2cc() + double S[3][3], Dm[6][6]; + get_pk2cc(F, S, Dm); + + // Compute dS using finite difference, given dF + if (order == 1){ + double F_tilde[3][3]; // perturbed deformation gradient + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + // Perturb the iJ-th component of F by delta * dF[i][J] + for (int k = 0; k < 3; k++) { + for (int l = 0; l < 3; l++) { + F_tilde[k][l] = F[k][l] + delta * dF[k][l]; + } + } + } + } + + // Compute perturbed S_tilde from F_tilde + double S_tilde[3][3], Dm_tilde[6][6]; + get_pk2cc(F_tilde, S_tilde, Dm_tilde); + + // Compute differences in S + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + dS[i][j] = S_tilde[i][j] - S[i][j]; + } + } + } + else if (order == 2){ + double F_plus[3][3]; // positive perturbed deformation gradient + double F_minus[3][3]; // negative perturbed deformation gradient + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + // Perturb the iJ-th component of F by +-delta * dF[i][J] + for (int k = 0; k < 3; k++) { + for (int l = 0; l < 3; l++) { + F_plus[k][l] = F[k][l] + delta * dF[k][l]; + F_minus[k][l] = F[k][l] - delta * dF[k][l]; + } + } + } + } + + // Compute perturbed S_plus and S_minus from F_plus and F_minus + double S_plus[3][3], Dm_plus[6][6]; + double S_minus[3][3], Dm_minus[6][6]; + get_pk2cc(F_plus, S_plus, Dm_plus); + get_pk2cc(F_minus, S_minus, Dm_minus); + + // Compute differences in S + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + dS[i][j] = S_plus[i][j] - S_minus[i][j]; + } + } + } + } + + /** + * @brief Compute material elasticity tensor contracted with perturbation in Green-Lagrange strain tensor (CC:dE) given perturbation in deformation gradient (dF) using finite differences + * + * @param F Deformation gradient + * @param dF Deformation gradient perturbation shape + * @param delta Deformation gradient perturbation scaling factor + * @param order Order of the finite difference scheme (1 for first order, 2 for second order, etc.) + * @param CCdE Material elasticity tensor times the perturbation in the Green-Lagrange strain tensor + */ + void calcCCdEFiniteDifference(const double F[3][3], const double dF[3][3], const double delta, const int order, double (&CCdE)[3][3]) { + // Compute CC(F) from get_pk2cc() + double S[3][3], Dm[6][6]; + get_pk2cc(F, S, Dm); + double CC[3][3][3][3]; + mat_models_carray::voigt_to_cc_carray<3>(Dm, CC); + + // Compute dE using finite difference, given dF + double dE[3][3]; + calcdEFiniteDifference(F, dF, delta, order, dE); + + // Compute CC:dE + mat_fun_carray::ten_mat_ddot<3>(CC, dE, CCdE); + } + + + /** + * @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 each component-wise perturbation 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] 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], double rel_tol, double abs_tol, double delta, bool verbose = false) { + int order = 2; + + // Compute E from F + double J, C[3][3], E[3][3]; + calc_JCE(F, J, C, E); + + // Compute S_ij(F) + // Compute CC_ijkl(F). + // CC is provided in Voigt notation as Dm, and we will convert it to CC + double S[3][3], Dm[6][6]; + get_pk2cc(F, S, Dm); // S from svFSI + + // 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++) { + EXPECT_NEAR(Dm_check[i][j], Dm[i][j], abs_tol); + } + } + // ------------------------------- - std::cout << "CC =" << std::endl; + // 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() + double dS[3][3], CCdE[3][3]; + + // Loop over perturbations to each component of F, dF + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + // Generate dF with 1.0 in (i,j) component + double dF[3][3] = {}; + dF[i][j] = 1.0; + + // Compute dS + calcdSFiniteDifference(F, dF, delta, order, dS); + + // Compute CC:dE + calcCCdEFiniteDifference(F, dF, delta, order, CCdE); + + // Check that CC_ijkl dE_kl = dS_ij 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++) { - std::cout << CC[i][j][k][l] << " "; + 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++) { + std::cout << F[i][j] << " "; + } + std::cout << std::endl; + } + + std::cout << "CC =" << std::endl; + 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++) { + std::cout << CC[i][j][k][l] << " "; + } + std::cout << std::endl; } - std::cout << std::endl; + } + std::cout << std::endl; + } + + std::cout << "dS =" << std::endl; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + std::cout << dS[i][j] << " "; + } + std::cout << std::endl; + } + + std::cout << "CCdE =" << std::endl; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + std::cout << CCdE[i][j] << " "; + } + std::cout << std::endl; + } + std::cout << std::endl; + } + } + } + } + + /** + * @brief Tests the order of convergence of the consistency between CC:dE and dS using finite differences. + * + * Analytically, we should have CC:dE = dS. This function determines the order of convergence of CC:dE = dS, where dE and dS are computed using finite differences in F. + * + * Pseudocode: + * - For each component-wise perturbation dF + * - For decreasing delta + * - Compute dS + * - Compute CC:dE + * - Compute error CC:dE - dS + * - Compute order of convergence by fitting a line to log(delta) vs log(error) + * + * Note that the order of convergence should be order + 1, because we are comparing differences (dS and CC:dE) + * instead of derivatives (e.g. dS/dF and CC:dE/dF). + * @param[in] F Deformation gradient. + * @param[in] delta_max Maximum perturbation scaling factor. + * @param[in] delta_min Minimum perturbation scaling factor. + * @param[in] order Order of the finite difference scheme (1 for first order, 2 for second order, etc.). + * @param[in] convergence_order_tol Tolerance for comparing convergence order with expected value + * @param[in] verbose Show values of errors and order of convergence if true. + */ + void testMaterialElasticityConsistencyConvergenceOrder(double F[3][3], double delta_max, double delta_min, int order, const double convergence_order_tol, bool verbose = false) { + // Check that delta_max > delta_min + if (delta_max <= delta_min) { + std::cerr << "Error: delta_max must be greater than delta_min." << std::endl; + return; + } + + // Check that order is 1 or 2 + if (order != 1 && order != 2) { + std::cerr << "Error: order must be 1 or 2." << std::endl; + return; + } + + // Create list of deltas for convergence test (delta = delta_max, delta_max/2, delta_max/4, ...) + std::vector deltas; + double delta = delta_max; + while (delta >= delta_min) { + deltas.push_back(delta); + delta /= 2.0; + } + + // Loop over perturbations to each component of F, dF + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + // Generate dF with 1.0 in (i,j) component + double dF[3][3] = {}; + dF[i][j] = 1.0; + + // Compute dS and CC:dE for each delta and store error in list + std::vector errors; + double dS[3][3], CCdE[3][3]; + for (int i = 0; i < deltas.size(); i++) { + calcdSFiniteDifference(F, dF, deltas[i], order, dS); + calcCCdEFiniteDifference(F, dF, deltas[i], order, CCdE); + + // Compute Frobenius norm of error between dS and CC:dE + double error = 0.0; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + error += pow(dS[i][j] - CCdE[i][j], 2); } } - std::cout << std::endl; + error = sqrt(error); + + // Store error in list + errors.push_back(error); } - - std::cout << "dE =" << std::endl; - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - std::cout << dE[i][j] << " "; - } - std::cout << std::endl; + + // Compute order of convergence by fitting a line to log(delta) vs log(error) + std::vector log_deltas, log_errors; + for (int i = 0; i < deltas.size(); i++) { + log_deltas.push_back(log(deltas[i])); + log_errors.push_back(log(errors[i])); } - - std::cout << "dS =" << std::endl; - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - std::cout << dS[i][j] << " "; + + // Fit a line to log(delta) vs log(error) + // m is the slope (order of convergence), b is the intercept + auto [m, b] = computeLinearRegression(log_deltas, log_errors); + + // Check that order of convergence is > (order + 1) - convergence_order_tol + EXPECT_GT(m, order + 1 - convergence_order_tol); + + // Print results if verbose + if (verbose) { + std::cout << "Iteration " << i << ":" << std::endl; + std::cout << "Slope (order of convergence): " << m << std::endl; + std::cout << "Intercept: " << b << std::endl; + std::cout << "Errors: "; + for (int i = 0; i < errors.size(); i++) { + std::cout << errors[i] << " "; } std::cout << std::endl; - } - - std::cout << "CCdE =" << std::endl; - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - std::cout << CCdE[i][j] << " "; + std::cout << std::endl; + + std::cout << "F = " << std::endl; + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + std::cout << F[i][J] << " "; + } + std::cout << std::endl; + } + + std::cout << "dF = " << std::endl; + for (int i = 0; i < 3; i++) { + for (int J = 0; J < 3; J++) { + std::cout << dF[i][J] << " "; + } + std::cout << std::endl; } std::cout << std::endl; } - std::cout << std::endl; } } } @@ -1011,6 +1495,124 @@ class TestMaterialModel { } }; + + +// -------------------------------------------------------------- +// --------------------- Material Parameters Classes ------------ +// -------------------------------------------------------------- + +// Class to contain material parameters +class MatParams { +public: + virtual ~MatParams() {} // Virtual destructor for proper cleanup +}; + +// Class to contain Neo-Hookean material parameters +class NeoHookeanParams : public MatParams { +public: + double C10; + + // Default constructor + NeoHookeanParams() : C10(0.0) {} + + // Constructor with parameters + NeoHookeanParams(double c10) : C10(c10) {} + +}; + +// Class to contain Mooney-Rivlin material parameters +class MooneyRivlinParams : public MatParams { +public: + double C01; + double C10; + + // Default constructor + MooneyRivlinParams() : C01(0.0), C10(0.0) {} + + // Constructor with parameters + MooneyRivlinParams(double c01, double c10) : C01(c01), C10(c10) {} + +}; + +// Class to contain Holzapfel-Ogden material parameters +class HolzapfelOgdenParams : public MatParams { +public: + double a; + double b; + double a_f; + double b_f; + double a_s; + double b_s; + double a_fs; + double b_fs; + double f[3]; // Fiber direction + double s[3]; // Sheet direction + + double k; // Smoothed Heaviside function parameter + + // 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), k(0.0) { + for (int i = 0; i < 3; i++) { + f[i] = 0.0; + s[i] = 0.0; + } + } + + // 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 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]; + } + } +}; + +// Class to contain Holzapfel-Ogden (Modified Anisortopy) material parameters +class HolzapfelOgdenMAParams : public MatParams { +public: + double a; + double b; + double a_f; + double b_f; + double a_s; + double b_s; + double a_fs; + double b_fs; + double f[3]; // Fiber direction + double s[3]; // Sheet direction + + double k; // Smoothed Heaviside function parameter + + // Default constructor + HolzapfelOgdenMAParams() : 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; + } + } + + // Constructor with parameters + HolzapfelOgdenMAParams(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]; + } + } +}; + +// 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) {} +}; + + // -------------------------------------------------------------- // --------------------- Material Model Classes ----------------- // -------------------------------------------------------------- @@ -1132,10 +1734,6 @@ class TestMooneyRivlin : public TestMaterialModel { * * 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: @@ -1165,6 +1763,141 @@ class TestHolzapfelOgden : public TestMaterialModel { dmn.stM.bss = params.b_s; dmn.stM.afs = params.a_fs; dmn.stM.bfs = params.b_fs; + dmn.stM.khs = params.k; // Smoothed Heaviside function parameter + dmn.stM.Kpen = 0.0; // Zero volumetric penalty parameter + + // Set number of fiber directions and fiber directions + nFn = 2; + Vector f = {params.f[0], params.f[1], params.f[2]}; + Vector s = {params.s[0], params.s[1], params.s[2]}; + fN.set_col(0, f); + fN.set_col(1, s); + } + + /** + * @brief Prints the Holzapfel-Ogden material parameters. + */ + void printMaterialParameters() { + std::cout << "a = " << params.a << std::endl; + std::cout << "b = " << params.b << std::endl; + std::cout << "a_f = " << params.a_f << std::endl; + std::cout << "b_f = " << params.b_f << std::endl; + std::cout << "a_s = " << params.a_s << std::endl; + 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 << "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; + } + + /** + * @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.))); + } + + /** + * @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); + + // Material parameters + double a = params.a; + double b = params.b; + double a_f = params.a_f; + double b_f = params.b_f; + double a_s = params.a_s; + double b_s = params.b_s; + double a_fs = params.a_fs; + double b_fs = params.b_fs; + + // Smoothed Heaviside parameter + double k = params.k; + + // Fiber and sheet directions + 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 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 implementation in svFSIplus + + // Invariants + double I1_bar = smTerms.Ib1; + // I4_bar_f = f . C_bar . f + double C_bar_f[3]; mat_fun_carray::mat_mul<3>(smTerms.C_bar, f, C_bar_f); + double I4_bar_f = mat_fun_carray::norm<3>(f, C_bar_f); + // I4_bar_s = s . C_bar . s + double C_bar_s[3]; mat_fun_carray::mat_mul<3>(smTerms.C_bar, s, C_bar_s); + double I4_bar_s = mat_fun_carray::norm<3>(s, C_bar_s); + // I8_bar_fs = f . C_bar . s + double I8_bar_fs = mat_fun_carray::norm<3>(f, C_bar_s); + + // Strain energy density for Holzapfel-Ogden material model with modified anisotropic invariants (bar quantities) + double Psi = 0.0; + Psi += a / (2.0 * b) * exp(b * (I1_bar - 3.0)); // Isotropic term + 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 + + return Psi; + } +}; + + +/** + * @brief Class for testing the Holzapfel-Ogden (Modified Anisotropy) material model. + * + * This class provides methods to set up and test the Holzapfel-Ogden-ma material + * model, including computing the strain energy and printing material parameters. + * + */ +class TestHolzapfelOgdenMA : public TestMaterialModel { +public: + + /** + * @brief Parameters for the Holzapfel-Ogden ma material model. + */ + HolzapfelOgdenMAParams params; + + /** + * @brief Constructor for the TestHolzapfelOgdenMA class. + * + * Initializes the Holzapfel-Ogden material parameters for svFSIplus. + * + * @param[in] params_ Parameters for the Holzapfel-Ogden ma material model. + */ + TestHolzapfelOgdenMA(const HolzapfelOgdenMAParams ¶ms_) : TestMaterialModel( consts::ConstitutiveModelType::stIso_HO_ma, consts::ConstitutiveModelType::stVol_ST91), + params(params_) + { + // Set Holzapfel-Ogden material parameters for svFSIplus + auto &dmn = com_mod.mockEq.mockDmn; + dmn.stM.a = params.a; + dmn.stM.b = params.b; + dmn.stM.aff = params.a_f; + dmn.stM.bff = params.b_f; + dmn.stM.ass = params.a_s; + dmn.stM.bss = params.b_s; + dmn.stM.afs = params.a_fs; + dmn.stM.bfs = params.b_fs; + dmn.stM.khs = params.k; // Smoothed Heaviside function parameter dmn.stM.Kpen = 0.0; // Zero volumetric penalty parameter // Set number of fiber directions and fiber directions @@ -1190,7 +1923,6 @@ class TestHolzapfelOgden : public TestMaterialModel { 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; } /** @@ -1240,55 +1972,27 @@ class TestHolzapfelOgden : public TestMaterialModel { // + 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 - double I1_bar = smTerms.Ib1; - // I4_f = f . C . f - double C_f[3]; mat_fun_carray::mat_mul<3>(smTerms.C, f, C_f); - double I4_f = mat_fun_carray::norm<3>(f, C_f); - // I4_s = s . C . s - double C_s[3]; mat_fun_carray::mat_mul<3>(smTerms.C, s, C_s); - double I4_s = mat_fun_carray::norm<3>(s, C_s); - // I8_fs = f . C . s - double I8_fs = mat_fun_carray::norm<3>(f, C_s); - - // Strain energy density for Holzapfel-Ogden material model with full anisotropic invariants - double Psi = 0.0; - Psi += a / (2.0 * b) * exp(b * (I1_bar - 3.0)); // Isotropic term - 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 - 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 - double I1_bar = smTerms.Ib1; - // I4_bar_f = f . C_bar . f - double C_bar_f[3]; mat_fun_carray::mat_mul<3>(smTerms.C_bar, f, C_bar_f); - double I4_bar_f = mat_fun_carray::norm<3>(f, C_bar_f); - // I4_bar_s = s . C_bar . s - double C_bar_s[3]; mat_fun_carray::mat_mul<3>(smTerms.C_bar, s, C_bar_s); - double I4_bar_s = mat_fun_carray::norm<3>(s, C_bar_s); - // I8_bar_fs = f . C_bar . s - double I8_bar_fs = mat_fun_carray::norm<3>(f, C_bar_s); - - // Strain energy density for Holzapfel-Ogden material model with modified anisotropic invariants (bar quantities) - double Psi = 0.0; - Psi += a / (2.0 * b) * exp(b * (I1_bar - 3.0)); // Isotropic term - 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 + // Invariants + double I1_bar = smTerms.Ib1; + // I4_f = f . C . f + double C_f[3]; mat_fun_carray::mat_mul<3>(smTerms.C, f, C_f); + double I4_f = mat_fun_carray::norm<3>(f, C_f); + // I4_s = s . C . s + double C_s[3]; mat_fun_carray::mat_mul<3>(smTerms.C, s, C_s); + double I4_s = mat_fun_carray::norm<3>(s, C_s); + // I8_fs = f . C . s + double I8_fs = mat_fun_carray::norm<3>(f, C_s); + + // Strain energy density for Holzapfel-Ogden material model with full anisotropic invariants + double Psi = 0.0; + Psi += a / (2.0 * b) * exp(b * (I1_bar - 3.0)); // Isotropic term + 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 + + return Psi; - return Psi; - } } };