Solid material model testing framework (#225)

* MR tests. Compare S directly. Compare S and Dm indirectly. Reorg unit tests a bit

* Simplify, separate, and rename S and CC tests

* Reorg test, add StVK, nHook tests. StVK doesn't pass for some reason

* Fix apparent bug in setting Mooney-Rivlin C10 parameter

* Add classes for material model params, classes for strain energies

* Add test fixtures, function for PK2Stress finite difference calc, make test names more descriptive

* Add delta parameters to calcPK2StressFiniteDifference

* Bug fix in Mooney-Rivlin model. Idx_prod was uninitialized, leading to undefined behavior.

* Add matrix-vector multiply for carrays

* Big reorg of S and CC tests. Currently only for Mooney-Rivlin material

* Add tests to compare S and S_ref directly, with prescribed S_ref and finite difference S_ref

* Add NeoHookean tests, add abs_tol

* Add Holzapfel-Ogden model test WITH BUGS

* Fix HO tests: Input fiber directions to get_pk2cc properly.
Create s orthogonal to f. Use F = I + random. Add flag for HO with/without full anisotropic invariants

* Add instructions for adding new material modeling tests

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

Code/Source/svFSI/mat_fun_carray.h

@@ -179,6 +179,20 @@ void mat_mul(const double A[N][N], const Vector<double>& v, double result[N])
  }
 }
 
+template <size_t N>
+void mat_mul(const double A[N][N], const double v[N], double result[N])
+{
+ for (int i = 0; i < N; i++) {
+ double sum = 0.0;
+
+ for (int j = 0; j < N; j++) {
+ sum += A[i][j] * v[j];
+ }
+
+ result[i] = sum;
+ }
+}
+
 template <size_t N>
 void mat_mul6x3(const double A[2*N][2*N], const Array<double>& B, Array<double>& C)
 {
@@ -203,6 +217,24 @@ void mat_mul6x3(const double A[2*N][2*N], const Array<double>& B, Array<double>&
  }
 }
 
+template <size_t N>
+void ten_mat_ddot(const double A[N][N][N][N], const double B[N][N], double C[N][N])
+{
+ for (int i = 0; i < N; i++) {
+ for (int j = 0; j < N; j++) {
+ double sum = 0.0;
+
+ for (int k = 0; k < N; k++) {
+ for (int l = 0; l < N; l++) {
+ sum += A[i][j][k][l] * B[k][l];
+ }
+ }
+
+ C[i][j] = sum;
+ }
+ }
+}
+
 template <size_t N>
 void mat_inv(const double A[N][N], double Ainv[N][N])
 {

Code/Source/svFSI/mat_models_carray.h

@@ -99,6 +99,67 @@ void cc_to_voigt_carray(const double CC[N][N][N][N], double Dm[2*N][2*N])
  } 
 }
 
+template <size_t N>
+void voigt_to_cc_carray(const double Dm[2*N][2*N], double CC[N][N][N][N]) {
+ if (N == 3) {
+ // Initialize the CC array with zeros
+ for (size_t i = 0; i < N; ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ for (size_t k = 0; k < N; ++k) {
+ for (size_t l = 0; l < N; ++l) {
+ CC[i][j][k][l] = 0.0;
+ }
+ }
+ }
+ }
+
+ // Voigt indices mapping
+ const int index_map[6][2] = {
+ {0, 0}, {1, 1}, {2, 2}, {0, 1}, {1, 2}, {2, 0}
+ };
+
+ // Fill in the CC array based on the Dm matrix
+ for (int I = 0; I < 6; ++I) {
+ for (int J = 0; J < 6; ++J) {
+ int i = index_map[I][0], j = index_map[I][1];
+ int k = index_map[J][0], l = index_map[J][1];
+ CC[i][j][k][l] = Dm[I][J];
+ CC[j][i][k][l] = Dm[I][J];
+ CC[i][j][l][k] = Dm[I][J];
+ CC[j][i][l][k] = Dm[I][J];
+ }
+ }
+ } else if (N == 2) {
+ // Initialize the CC array with zeros
+ for (size_t i = 0; i < N; ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ for (size_t k = 0; k < N; ++k) {
+ for (size_t l = 0; l < N; ++l) {
+ CC[i][j][k][l] = 0.0;
+ }
+ }
+ }
+ }
+
+ // Voigt indices mapping for 2D
+ const int index_map[3][2] = {
+ {0, 0}, {1, 1}, {0, 1}
+ };
+
+ // Fill in the CC array based on the Dm matrix
+ for (int I = 0; I < 3; ++I) {
+ for (int J = 0; J < 3; ++J) {
+ int i = index_map[I][0], j = index_map[I][1];
+ int k = index_map[J][0], l = index_map[J][1];
+ CC[i][j][k][l] = Dm[I][J];
+ CC[i][j][l][k] = Dm[I][J];
+ CC[j][i][k][l] = Dm[I][J];
+ CC[j][i][l][k] = Dm[I][J];
+ }
+ }
+ }
+}
+
 //-----------
 // get_pk2cc
 //-----------
@@ -382,7 +443,7 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn
  g1 = 4.0 * J4d * stM.C01;
 
  double CCb[N][N][N][N];
- mat_fun_carray::ten_dyad_prod<N>(Idm, Idm, CCb);
+ mat_fun_carray::ten_dyad_prod<N>(Idm, Idm, Idm_prod);
 
  for (int i = 0; i < nsd; i++) {
  for (int j = 0; j < nsd; j++) {

Code/Source/svFSI/set_material_props.h

@@ -97,7 +97,7 @@ SeMaterialPropertiesMapType set_material_props = {
 {
  lDmn.stM.isoType = consts::ConstitutiveModelType::stIso_MR;
  lDmn.stM.C10 = domain_params->constitutive_model.mooney_rivlin.c1.value();
- lDmn.stM.C01 = lDmn.stM.C10 = domain_params->constitutive_model.mooney_rivlin.c2.value();
+ lDmn.stM.C01 = domain_params->constitutive_model.mooney_rivlin.c2.value();
 } },
 
 //---------------------------//